Semiclassical Dynamics and Nonlinear Charge Current

TL;DR

This paper reviews the semiclassical framework for understanding nonlinear charge currents in electron systems, which are key to probing material properties and aiding band-structure engineering.

Contribution

It provides a universal second-order semiclassical approach to analyze nonlinear charge currents based on unperturbed states.

Findings

Unified semiclassical treatment of nonlinear charge currents

Clarification of symmetry requirements for different currents

Facilitates band-structure engineering of nonlinear responses

Abstract

Electron conductivity is an important material property that can provide a wealth of information about the underlying system. Especially, the response of the conductivity with respect to electromagnetic fields corresponds to various nonlinear charge currents, which have distinct symmetry requirements and hence can be used as efficient probes of different systems. To help the band-structure engineering of such nonlinear currents, a universal treatment of electron dynamics up to second order expressed in the basis of the unperturbed states are highly useful. In this work, we review the general semiclassical framework of the nonlinear charge currents.

Tables

Table 1. Table 1: Nonzero components of conductivity in crystals with different space inversion and time reversal symmetries. The requirement of the lower right entry should be understood as in crystals that simultaneously break time reversal and inversion symmetry but preserve the combined symmetry.

	Even	Odd
Even	${\sigma }_{(ij,k),0}$ ${\sigma }_{(ij,k),2}$ ${\sigma }_{(ij,kk),1}$  ${\sigma }_{(ij,kk),3}$	${\sigma }_{ijk,1}$
Odd	${\sigma }_{(ij,k),1}$ ${\sigma }_{(ij,kk),0}$ ${\sigma }_{(ij,kk),2}$	${\sigma }_{ijk,0}$ ${\sigma }_{ijk,2}$

Equations

$$J_i={\sigma }_{ij}{E}_j+{\sigma }_{ijk}{E}_j{E}_k+{\sigma }_{(ij,k)}{E}_j{B}_k+{\sigma }_{(ij,kk)}{E}_j{B}_k{B}_k+\cdots .$$

$${\sigma }_{ijk}$$

$${\sigma }_{(ij,k)}$$

$${\sigma }_{(ij,kk)}$$

$${\hat{H}}_{\text{f}}={\hat{H}}_0[\hat{\mathbf{p}}+e\mathbf{A}(\mathbf{r});\mathbf{r}]+\frac{g{\mu }_B}{\hslash }\mathbf{B}\cdot \hat{\mathbf{s}}+e\mathbf{E}\cdot \mathbf{r},$$

$${\hat{H}}_c={\hat{H}}_0[\hat{\mathbf{p}}+e\mathbf{A}({\mathbf{r}}_c);\mathbf{r}]+\frac{g{\mu }_B}{\hslash }\mathbf{B}\cdot \hat{\mathbf{s}}+e\mathbf{E}\cdot {\mathbf{r}}_c,$$

$$|W⟩=\int d\mathbf{p}{C}_0(\mathbf{p})e^{i\mathbf{p}\cdot \mathbf{r}}|u_0[\mathbf{p}+e\mathbf{A}({\mathbf{r}}_c)/\hslash ]⟩.$$

$${\mathbf{p}}_c=⟨W|\mathbf{p}|W⟩.$$

$${\mathbf{r}}_c=⟨W|\mathbf{r}|W⟩={\frac{\partial \gamma }{\partial \mathbf{p}}|}_{\mathbf{p}={\mathbf{p}}_c}+{\mathcal{A}}_0({\mathbf{k}}_c),$$

$$L$$

$$=-({\mathbf{r}}_c-{\mathcal{A}}_0)\cdot \hslash {\dot{\mathbf{k}}}_c-\frac{1}{2}e\mathbf{B}\times {\mathbf{r}}_c\cdot {\dot{\mathbf{r}}}_c-{\tilde{\epsilon }}_0,$$

$$L=\left(\hslash {\mathbf{k}}_c-\frac{1}{2}e\mathbf{B}\times {\mathbf{r}}_c\right)\cdot {\dot{\mathbf{r}}}_c-{\tilde{\epsilon }}_0.$$

$${\mathbf{r}}_c$$

$${\mathbf{k}}_c$$

$$L=\hslash \mathbf{p}\cdot \dot{\mathbf{q}}-{\tilde{\epsilon }}_0.$$

$$\mathcal{D}=1+\frac{e}{\hslash }\mathbf{B}\cdot {\mathbf{\Omega }}_0({\mathbf{k}}_c).$$

$${\dot{\mathbf{r}}}_c$$

$$\hslash {\dot{\mathbf{k}}}_c$$

$$\mathbf{P}=-e\int \frac{d{\mathbf{k}}_c}{(2\pi {)}^3}{\mathcal{A}}_0({\mathbf{k}}_c).$$

$${\mathbf{m}}_{\text{orb}}$$

$$=-\frac{e}{2}\sum _{n\ne 0}{\mathcal{A}}_{0n}\times {\mathbf{v}}_{n0},$$

$$F=\int \mathcal{D}\frac{d\mathbf{k}}{(2\pi {)}^3}(-k_{B}T)\ln \left[1+\exp \left(\frac{{\tilde{\epsilon }}_0-\mu }{k_{B}T}\right)\right].$$

$$\mathbf{M}=-\frac{\partial F}{\partial B}=\int \frac{d\mathbf{k}}{(2\pi {)}^3}\left(f_0\mathbf{m}-\frac{e}{\hslash }{g}_e{\mathbf{\Omega }}_0\right),$$

$$|W⟩=\int d\mathbf{p}{e}^{i\mathbf{p}\cdot \mathbf{r}}\left(C_0(\mathbf{p})|u_0⟩+\sum _{n\ne 0}{C}_n(\mathbf{p})|u_n⟩\right).$$

$$⟨u_n|e^{-i\mathbf{p}\cdot \mathbf{r}}(i\hslash {\partial }_t-{\hat{H}}_{\text{f}})|W⟩=0,\forall n\ne 0.$$

$$C_n=\frac{G_{n0}}{{\epsilon }_0-{\epsilon }_n}{C}_0-\frac{i}{2}\frac{e}{\hslash }[\mathbf{B}\times (i{\partial }_{\mathbf{p}}+{\mathcal{A}}_0-{\mathbf{r}}_c)C_0]\cdot {\mathcal{A}}_{n0},$$

$$\int d\mathbf{p}{e}^{i\mathbf{p}\cdot \mathbf{r}}\frac{-i}{2}\frac{e}{\hslash }[\mathbf{B}\times (i{\partial }_{\mathbf{p}}+{\mathcal{A}}_0-{\mathbf{r}}_c)C_0]\cdot \sum _{n\ne 0}{\mathcal{A}}_{n0}|u_n⟩$$

$$=$$

$$L=-({\mathbf{r}}_c-{\mathcal{A}}_0^{\text{t}})\cdot \hslash {\dot{\mathbf{k}}}_c-\frac{1}{2}e\mathbf{B}\times {\mathbf{r}}_c\cdot {\dot{\mathbf{r}}}_c-{\tilde{\epsilon }}_0,$$

$${\mathbf{r}}_c$$

$${\mathbf{k}}_c$$

$$\mathcal{D}=1+\frac{e}{\hslash }\mathbf{B}\cdot {\mathbf{\Omega }}_0^{\text{t}}.$$

$${\dot{\mathbf{r}}}_c$$

$$\hslash {\dot{\mathbf{k}}}_c$$

$${\mathbf{r}}_c=⟨W|\mathbf{r}|W⟩={\frac{\partial \gamma }{\partial \mathbf{p}}|}_{\mathbf{p}={\mathbf{p}}_c}+{\mathcal{A}}_0({\mathbf{k}}_c)+{\mathcal{A}}_0^{\prime }({\mathbf{k}}_c).$$

Full text

I Introduction

Electron transport has long been a central topic in solid state physics, not only due to its fascinating complexity but also because it can provide rich information about the underlying electronic structure and dynamics. In the most simple Drude theory of metals, the conductivity is determined by the density, effective mass of carriers, and the phenomenological relaxation time Ashcroft and Mermin (1976). Such simple picture can be generalized using the dynamics of Bloch electrons, to take into account the real energy spectrum and the Fermi surface geometry, as well as the microscopic scattering process Ashcroft and Mermin (1976).

In recent years, this classic paradigm of electron transport receives profound modifications, due to a deeper understanding of the geometric structure in the momentum space induced by the wave functions instead of the energy spectrum. The most well-known examples are various Hall effects Nagaosa et al. (2010); Xiao et al. (2010). Central to the Hall effect is the concept of the Berry curvature. Its integration yields the Chern number, responsible for a topological classification of matter Qi and Zhang (2011). Moreover, it directly modifies the electron motion and hence is indispensable in describing the static and optical transport phenomenon Xiao et al. (2010); Morimoto and Nagaosa (2016); Zhong et al. (2016); Ma and Pesin (2015). Berry curvature is also the analog of the magnetic field in the momentum space. As such, the source of Berry curvature is the momentum-space analog of magnetic monopoles, which can be realized by Weyl/Dirac points Son and Spivak (2013); Armitage et al. (2018). In three dimensions, such isolated monopoles in momentum space can response differently to electromagnetic fields and cause a characteristic strong negative longitudinal magnetoresistance Son and Spivak (2013); Armitage et al. (2018).

On the other hand, the charge current can be expanded with respect to the electromagnetic fields as follows

[TABLE]

Berry curvature offers a first order correction to the electron dynamics. As such, it can give an accurate description of $\sigma_{ij}$ whose antisymmetric part is the anomalous Hall effect. In addition to this linear current, there are various nonlinear currents with coefficients satisfying different symmetry requirements and hence can be used to probe distinct systems. To further demonstrate the symmetry requirement, we can expand the last three coefficients in terms of transport relaxation time $\tau$ 111the highest order of $\tau$ is the same as the highest order of the field, which is the property of the asymptotic solution to the Boltzmann equation, as implied in Eq. (64).

[TABLE]

The symmetry of different coefficients is summarized in Table. 1. It is thus tempting to find how different coefficients depend on the geometric properties of Bloch electrons in the momentum space. This can provide experimental guide and help the band structure engineering of larger nonlinear charge currents. For this purpose, the electron dynamics should be extended to second order.

In this work, we will focus on the semiclassical treatment of the electron dynamics and review the effort of generalizing such theory up to second order and deriving the nonlinear currents in the Boltzmann transport framework. Our paper is organized as follows. In Sect. II, we discuss the generalization of the semiclassical equations of motion up to second order through an appropriately constructed wave packet up to first order. We then discuss the validity of the semiclassical theory though the magnetoelectric coefficient and magnetic susceptibility. With the correct semiclassical theory up to second order, in Sect. III, we extend the Boltzmann transport theory and derive its solution in steady-state up to third order. We then use this universal treatment to derive various nonlinear currents, including the electric-field-correction to the nonlinear Hall effect ($\sigma_{jii}$), linear magnetoresistance ($\sigma_{(ij,k)}$), and quadratic magnetoresistance ($\sigma_{(ij,kk)}$).

II Semiclassical Theory of Bloch electrons in electromagnetic fields

To derive various nonlinear charge currents, one will need the correct electron dynamics beyond linear order under electromagnetic fields. For this purpose, certain perturbation technique should be employed. However, the periodic nature of crystals causes difficulties for such process. We note that the perturbation from electromagnetic fields reads $e\bm{E}\cdot\bm{r}+\frac{1}{2}e[\bm{A}(\bm{r})\cdot\hat{\bm{v}}+\hat{\bm{v}}\cdot\bm{A}(\bm{r})]$. It contains the position operator, which is unbounded in extended systems. Therefore, direct application of the time-independent perturbation is inappropriate.

To circumvent such difficulty, for magnetoconductivty one can first obtain the eigenstate of Bloch electrons under magnetic field, which are Landau levels for effective continuum models or Hofstadter spectrum for tight-binding models. Then the static magnetoconductivity is obtained using the standard Kubo formula under an appropriate limit. The dependence of the conductivity on magnetic field is encoded in the quantization of original Bloch states under magnetic field. Although this method is standard and has wide applicability, it cannot give universal understanding of the physics behind the magnetoconductivity, as the quantization into Landau levels or the Hofstadter spectrum varies for different systems.

Alternatively, one can construct a phase space using an appropriate position variable $\bm{R}$ and momentum variable $\bm{P}$ from the Bloch states. The electron dynamics is then contained in the equations of motion in the phase space. If $\bm{R}$ and $\bm{P}$ are canonical variables, their dynamics is solely determined by the energy. Therefore, the whole procedure amounts to deriving an effective energy in terms of $\bm{R}$ and $\bm{P}$ under external fields. If $\bm{R}$ and $\bm{P}$ are non-canonical, their Poisson brackets are also needed, which determine the structure of the equations of motion. The magnetoconductivity is then derived in combination with the Boltzmann equation for the distribution function. Although this method treats external fields as perturbations, and is only derived up to the second order currently, it has the benefit that it yields analytical results independent of the model detail. In fact, all terms in the equations of motion are expressed using the unperturbed Bloch functions, and hence can directly access both the spectrum and the geometric structure in the unperturbed Hilbert space.

In this section, we follow the second method by treating Bloch electrons as wave packets, whose center of mass position and momentum naturally span the phase space. We will first sketch the construction of the wave packet and show how its internal structure modifies the electron dynamics in phase space at first order of external fields. A detailed discussion of both aspects can be found in Ref. Xiao et al. (2010). We then extend such theory up to second order and derive the second order correction to equations of motion. The validity of the semiclassical theory can be confirmed through the magnetoelectric effect and the quantization of effective energies into Landau levels under magnetic field. Since all elements in the semiclassical dynamics are expressed using unperturbed Bloch states, they can be evaluated using the output from the first-principles codes.

II.1 Properties of the wave packet

The dynamics of Bloch electrons under uniform electromagnetic fields is governed by the following crystal Hamiltonian in non-relativistic limit

[TABLE]

where $\hat{H}_{0}(\hat{\bm{p}};\bm{r})$ is the unperturbed Hamiltonian for periodic crystals with $\hat{\bm{p}}$ being the momentum operator, $g$ is the gyromagnetic ratio, $\mu_{B}$ is the Bohr magneton, and $\hat{\bm{s}}$ is the spin operator. The magnetic field enters through both the minimal coupling and the Zeeman coupling. Generally speaking, $\hat{H}_{\text{f}}$ does not respect the translational symmetry, unless the magnetic flux through a unit cell is a rational number times the flux quantum such that the magnetic translational symmetry is recovered Zak (1964a, b); Hofstadter (1976); Chang and Niu (1995).

Although $\hat{H}_{\text{f}}$ exactly determines the dynamics, to directly diagonalize it is a difficult task. Instead, we can simplify the problem by assuming that the solution to the Schr$\ddot{\rm o}$dinger equation $(i\hbar\partial_{t}-\hat{H}_{\text{f}})\psi=0$ takes the ansartz of a wave packet, which is the essence of the semiclassical theory. To construct the wave packet, we first evaluate the electromagnetic potential in the full Hamiltonian at the center of mass position $\bm{r}_{c}$ of the wave packet, and hence recover the translational symmetry of the original unperturbed crystal. The resulting local Hamiltonian reads:

[TABLE]

Its eigenenergy forms the local Bloch bands $\varepsilon_{n}[\bm{p}+(e/\hbar)\bm{A}(\bm{r}_{c})]+e\bm{E}\cdot\bm{r}_{c}$ with $n$ being the band index and $\hbar\bm{p}$ being crystal momentum. Its eigenfunction is the Bloch function $e^{i\bm{p}\cdot\bm{r}}|u_{n}[\bm{p}+e\bm{A}(\bm{r}_{c})/\hbar]\rangle$.

In the following, we focus on a single band with index [math] and assume that it is well separated from all the other bands. When electromagnetic fields are weak, the dynamics of a Bloch electron starting from some state in this band will still be confined in the same band. Therefore, the wave packet is constructed as the superposition of Bloch states from the band [math] Sundaram and Niu (1999)

[TABLE]

$|W\rangle$ should be normalized, indicating $\int d\bm{p}|C_{0}|^{2}=1$.

There are two constraints on the coefficient $C_{0}$. First, the wave packet is assumed to be sharply localized in the momentum space. This requires that the magnitude of $C_{0}$ satisfies $|C_{0}|^{2}\approx\delta(\bm{p}-\bm{p}_{c})$. Here $\hbar\bm{p}_{c}$ is the center of mass momentum of the wave packet, i.e.

[TABLE]

Secondly, the construction of the wave packet requires a local Hamiltonian around the center of mass position $\bm{r}_{c}$. This $\bm{r}_{c}$ has to be determined in a self-consistent manner, which exerts a constraint on the phase of $C_{0}$ Sundaram and Niu (1999):

[TABLE]

where $\gamma=-arg(C_{0})$, $\bm{\mathcal{A}}_{0}(\bm{p})=\langle u_{0}(\bm{p})|i\bm{\partial}_{\bm{p}}|u_{0}(\bm{p})\rangle$ is the intraband Berry connection, and $\hbar\bm{k}_{c}=\hbar\bm{p}_{c}+e\bm{A}(\bm{r}_{c})$ is the gauge-independent physical momentum. These two constraints complete the construction of the wave packet. Although the wave packet is for Bloch states from a single band, it can be generalized to account for the general multiband case Culcer et al. (2005); Shindou and Imura (2005).

II.2 Semiclassical dynamics up to first order

With the knowledge of the wave packet, we now derive its evolution. The coefficient $C_{0}$ in the wave packet can be determined through the variational principle, i.e. the least action principle. Since $C_{0}$ is a complex function, with independent magnitude and phase, determined by $\bm{k}_{c}$ (or equivalently, $\bm{p}_{c}$) and $\bm{r}_{c}$ respectively, applying the variational principle to the Lagrangian will yield the dynamics of $\bm{r}_{c}$ and $\bm{k}_{c}$, i.e. the phase space equations of motion.

By evaluating the Lagrangian under the wave packet up to first order, we obtain Sundaram and Niu (1999)

[TABLE]

where $\hat{H}_{1}=\frac{e}{4}\bm{B}\cdot[(\bm{r}-\bm{r}_{c})\times\hat{\bm{v}}-\hat{\bm{v}}\times(\bm{r}-\bm{r}_{c})]+e\bm{E}\cdot(\bm{r}-\bm{r}_{c})$ is the first order correction to the local Hamiltonian $\hat{H}_{c}$ and $\tilde{\varepsilon}=\varepsilon_{0}-\bm{B}\cdot\bm{m}+e\bm{E}\cdot\bm{r}_{c}$ with $\bm{m}=\bm{m}_{orb}-(g\mu_{B}/\hbar)\langle u_{0}|\hat{\bm{s}}|u_{0}\rangle$ is the modified Band energy after taking account of both the magnetic moment and electric dipole of the wave packet, which will be discussed later. Here we use the symmetric gauge for the vector potential.

The Berry connection $\bm{\mathcal{A}}_{0}$ is an essential ingredient in the Lagrangian. In the absence of $\bm{\mathcal{A}}_{0}$, the above Lagrangian reduces to the conventional canonical form for electrons under electromagnetic fields:

[TABLE]

The appearance of $\bm{\mathcal{A}}_{0}$ indicates that $\bm{r}_{c}$ and $\bm{k}_{c}$ are not canonical variables. To make this statement clearer, we make the following transformation Chang and Niu (2008)

[TABLE]

where $\bm{\Omega}_{0}=\bm{\nabla}_{\bm{p}}\times\bm{\mathcal{A}}_{0}$ is the Berry curvature, which can be viewed as the geometrical Berry phase per unit area Berry (1984). Here the argument of $\bm{\mathcal{A}}_{0}$ and $\bm{\Omega}_{0}$ is $\bm{p}+\frac{e}{2\hbar}\bm{B}\times\bm{q}$. Then the Lagrangian in Eq. (II.2) recovers the canonical form for $\bm{p}$ and $\bm{q}$

[TABLE]

Eq. (12) and (13) thus describe the connection between physical variables and canonical variables.

One consequence of the noncanonicality is that the phase space measure for the volume element $d\bm{r}_{c}d\bm{k}_{c}$ has to change based on the Jacobian $\partial(\bm{r}_{c},\bm{k}_{c})/\partial(\bm{q},\bm{p})$. The resulting phase space density of states reads Xiao et al. (2005)

[TABLE]

Since the physical position and momentum take the crystal volume and the Brillouin zone as their range, it is natural to evaluate the statistical average of any operator in the physical phase space spanned by $(\bm{r}_{c},\bm{k}_{c})$. Therefore, the modified density of states $\mathcal{D}$ should always be used if the magnetic field is involved.

The noncanonicality of $\bm{r}_{c}$ and $\bm{k}_{c}$ will also lead to a nontrivial symplectic form, and hence affects the structure of the equations of motion. From Eq. (II.2), the Euler-Lagrangian equations of motion yield the following phase space dynamics Sundaram and Niu (1999)

[TABLE]

The second term in the velocity equation is the anomalous velocity arising from the noncanonicality. The above equations of motion have wide applicability in the static spin, heat, and charge transport phenomena Xiao et al. (2010) and optical phenomena Moore and Orenstein (2010); Ma and Pesin (2015); Zhong et al. (2015); Morimoto and Nagaosa (2016); Zhong et al. (2016); Morimoto et al. (2016); Deyo et al. (2019).

II.3 Electromagnetic Dipoles

One direct application of the semiclassical theory is to derive the electromagnetic multipoles in crystals. This derivation confirms the validity of the semiclassical theory as the results are consistent with those derived using other methods. As the wave packet is sharply localized in the momentum space, it has a finite width in the real space. The charge and current distributions across such spread need not to be uniform. To represent such internal anisotropy of the wave packet, it is natural to use electromagnetic multipoles. For illustration purpose, we will consider the electromagnetic dipoles.

The electric dipole of the wave packet reads $\langle W|-e\bm{r}|W\rangle=-e\bm{r}_{c}$, with $\bm{r}_{c}$ given in Eq. (9). Under the periodic gauge Resta (2000); Xiao et al. (2010), the Bloch wave function is periodic in the momentum space, so is the coefficient $C_{0}$. Therefore, by integrating $\bm{r}_{c}$ in the Brillouin zone, the first term in Eq. (9) from the phase of $C_{0}$ can at most yield an integer multiples of $2\pi$ and hence can be ignored. The remaining term yields the electric polarization in crystals

[TABLE]

Originally, the electric polarization is derived by first considering the adiabatic charge pumping current and then identifying it as the displacement current from the temporal variation of the polarization King-Smith and Vanderbilt (1993); Resta (1994). Here our method yields the same electric polarization under the periodic gauge for Bloch functions.

The magnetic dipole can be evaluated as follows Chang and Niu (1996); Sundaram and Niu (1999)

[TABLE]

where $\bm{\mathcal{A}}_{0n}=\langle u_{0}|i\bm{\partial}_{\bm{k}_{c}}|u_{n}\rangle$ is the interband Berry connection and $\bm{v}_{n0}=\langle u_{n}|\hat{\bm{v}}|u_{0}\rangle$ is the interband velocity element. In this work, we use the convention that for intraband matrix elements, we only keep one band index. The magnetic moment $\bm{m}_{\text{orb}}$ represents how a wave packet couples to an external magnetic field, i.e. shifting the band energy through an effective Zeeman coupling $-\bm{B}\cdot\bm{m}_{\text{orb}}$. As such, $\bm{m}_{\text{orb}}$ contributes to the orbital magnetization as shown later. $\bm{m}_{\text{orb}}$ can usually be accessed through the optical activity of crystals. For example, in ferromagnetic materials, the cross-gap part of $\bm{m}_{\text{orb}}$ (with [math] standing for occupied bands and $n$ for the unoccupied bands in Eq. (II.3)) can be measured through the $f$-sum rule for the circular dichroism Souza and Vanderbilt (2008); Yao et al. (2008). Moreover, in noncentrosymmetric metals, the natural optical activity in the semiclassical regime is determined by the integration of the momentum space dipole of the magnetic moment (defined as $v_{i}(m_{\text{orb}})_{j}$ with $\bm{v}$ being the band velocity) on the Fermi surface Ma and Pesin (2015); Zhong et al. (2015, 2016).

Interestingly, if one starts from the relativistic Dirac Hamiltonian and construct a coherent wave packet for the upper two bands, the orbital magnetic moment of this wave packet naturally recovers the spin magnetic moment Chang and Niu (2008). Mathematically, this is equivalent to the Foldy-Wouthuysen transformation that reduces the Dirac Hamiltonian to the non-relativistic Schr$\ddot{\rm o}$dinger Hamiltonian Foldy and Wouthuysen (1950); Blount (1962a). Here the semiclassical theory offers a heuristic and clear picture for the emergence of spin.

To derive the orbital magnetization in crystals, one should evaluate the free energy under magnetic field. On one hand, the magnetic moment of the wave packet couples to the magnetic field and modifies the band energy. On the other hand, the phase space density of states is changed by the magnetic field. Taking both corrections into consideration, one obtains the following free energy

[TABLE]

Therefore, the orbital magnetization can be derived by taking the derivative of the free energy with respect to the magnetic field Xiao et al. (2006)

[TABLE]

where $g_{e}=-(k_{B}T)\ln[1+\exp((\varepsilon_{0}-\mu)/k_{B}T)]$ is the grand potential density, and $f_{0}$ is the equilibrium Fermi distribution function. Here and hereafter we will ignore the subscript $c$ in the integration over $\bm{k}_{c}$. This is the contribution of the orbital magnetization from band [math]. The total contribution can be obtained by summing over the band index. The orbital magnetization thus contains two contributions. The first one is due to the relative motion inside the wave packet, or the self-rotation of the wave packet. The second one is due to the rotation of the wave packet as a whole, or the revolution of the wave packet. This same magnetization can also be obtained through the Wannier function approach Thonhauser et al. (2005); Ceresoli et al. (2006), the exact Hofstadter spectrum Gat and Avron (2003a, b), or the linear response theory Shi et al. (2007); Qin et al. (2011).

II.4 Semiclassical dynamics up to second order

Although the semiclassical dynamics in Eq. (16) and (17) is very powerful, it cannot fully account for the second order response functions such as the magnetoresistance, magnetoelectric coefficient, and so on. For this purpose, we need to extend the above semiclassical dynamics up to second order.

The construction of the wave packet has to be modified accordingly. The external electromagnetic fields generally affect the eigenfunction. As a result, the wave packet should be the superposition of the true eigenstates instead of the unperturbed ones. We can still expand the true eigenstate in the basis of the unperturbed ones, and obtain the following wave packet Gao et al. (2014)

[TABLE]

Here for simplicity we drop the argument $\bm{p}+(e/\hbar)\bm{A}(\bm{r}_{c})$ of $|u_{0}\rangle$ and $|u_{n}\rangle$. If $|W\rangle$ satisfies the Schr$\ddot{\rm o}$dinger equation, we have

[TABLE]

Therefore, we use this as an additional constraint on the wave packet in Eq. (22), which determines the connection between $C_{n}$ and $C_{0}$. Specifically, solving this constraint at the linear order of electromagnetic fields yields Gao et al. (2014)

[TABLE]

where $G_{n0}=-\bm{B}\cdot\bm{M}_{n0}+e\bm{E}\cdot\bm{\mathcal{A}}_{n0}$ with $\bm{M}_{n0}=\frac{e}{2}\sum_{m\neq 0}(\bm{v}_{nm}+\bm{v}_{0}\delta_{mn})\times\bm{\mathcal{A}}_{m0}-g\mu_{B}\bm{s}_{n0}/\hbar$ being the interband element of the total (spin plus orbital) magnetic moment. We comment that different from the discussion here, in Ref. Gao et al. (2014) and Gao et al. (2015), the spin magnetic moment has not been explicitly added.

Equation (24) indicates that the correction to the wave packet from the magnetic field has two contributions with distinct structures and hence different origins. It is useful to first analyze the nature of these two perturbations as they prelude the form of the second order correction to the phase space dynamics. The first term has the conventional form of the perturbation correction with the energy gap as the denominator. It is called the vertical mixing as it mixes Bloch states from different bands but at the same $\bm{p}$ point in the Brillouin zone. It also contains both adiabatic and nonadiabatic perturbation, as the wave packet $|W\rangle$ is time-dependent through the argument $\bm{p}+(e/\hbar)\bm{A}(\bm{r}_{c})$ of the Bloch state.

In comparison, the second term in $C_{n}$ modifies the wave packet in the following way:

[TABLE]

where $\hat{D}=\bm{\partial}_{\bm{p}}+i\bm{\mathcal{A}}_{0}$ is the covariant derivative that ensures the gauge invariance. Notice that the momentum argument of $|u_{0}\rangle$ is $\bm{p}+(e/\hbar)\bm{A}(\bm{r}_{c})$. The first term in Eq. (II.4) has the meaning of shifting the momentum to $\bm{p}+(e/\hbar)\bm{A}(\bm{r})$. This is a unique property of the perturbation from the magnetic field, and it suggests that the correction to the wave function also obeys the Peierls substitution. The same property is also derived using the Moyal product and phase space formulation of quantum mechanics Blount (1962b). The second term in Eq. (II.4) makes the total correction normal to the original state $e^{i\bm{p}\cdot\bm{r}}|u_{0}\rangle$ to eliminate redundancy in the perturbative correction. Since the second term in $C_{n}$ effectively mixes Bloch states in the same band but at neighbouring $\bm{p}$ points, it is referred to as the horizontal mixing. Both corrections can be easily visualized in Fig. 1.

Using the modified wave packet, we obtain the following Lagrangian

[TABLE]

where $\bm{\mathcal{A}}_{0}^{\text{t}}$ is the Berry connection evaluated in the true eigen state instead of the unperturbed one. It can be put in the form of the original Berry connection plus a positional shift correction: $\bm{\mathcal{A}}_{0}^{\text{t}}=\bm{\mathcal{A}}_{0}+\bm{\mathcal{A}}_{0}^{\prime}$. Before discussing the expression of the positional shift $\bm{\mathcal{A}}_{0}^{\prime}$ and modified energy $\tilde{\varepsilon}$, we comment that the Lagrangian up to second order has the exact same form with the previous one in Eq. (II.2). Consequently, the connection between physical and canonical variables also has the same form except that the Berry connection $\bm{\mathcal{A}}_{0}^{\text{t}}$ in the exact eigenstate should be used

[TABLE]

where $\bm{\Omega}_{0}^{\text{t}}=\bm{\nabla}_{\bm{p}}\times\bm{\mathcal{A}}_{0}^{\text{t}}$ is the Berry curvature evaluated using the exact eigenstate. As a result, the phase space density of states reads

[TABLE]

More importantly, the equations of motion up to second order also keep the same form Gao et al. (2014)

[TABLE]

The positional shift is gauge-independent. It represents the additional shift of the wave packet center after applying external electromagnetic fields

[TABLE]

The positional shift at first order reads

[TABLE]

where $\epsilon_{ijk}$ is the Levi-Civita symbol and $\gamma_{ji\ell}=\frac{1}{2}(\partial_{\ell}g_{ij}+\partial_{i}g_{j\ell}-\partial_{j}g_{i\ell})$ is the Christoffel symbol. Here and hereafter, we use the Einstein summation convention for repeated indices and the partial derivative is with respect to the crystal momentum unless otherwise specified. $g_{ij}$ is a Fubini-Study metric tensor, usually called quantum metric Provost and Vallee (1980); Neupert et al. (2013); Anandan and Aharonov (1990); Resta (2011). It measures the distance between two neighbouring Bloch states in the momentum space. In fact, it is a part of a more general concept, the quantum geometrical tensor, which, for band [math], is defined as

[TABLE]

The real and imaginary part of $\mathcal{G}_{ij}$ yield the quantum metric and Berry curvature, respectively Provost and Vallee (1980); Neupert et al. (2013).

In the positional shift, the first term is due to the vertical mixing and the second term is due to the horizontal mixing. It is interesting that the horizontal mixing yields a purely geometric correction to the Berry connection. This is not coincidental, as the horizontal mixing addresses neighbouring Bloch states in the momentum space, whose difference is proportional to the Berry connection and whose distance is determined by the quantum metric. We comment that this positional shift is also envisioned in Ref. Berry (1987), using a similar technique.

The positional shift reflects the change in the electric dipole moment of the wave packet. When external fields are applied, the electric polarization has to pick up the following change. First, the density of states has to be changed to $\mathcal{D}$. Second, the Berry connection $\bm{\mathcal{A}}_{0}\rightarrow\bm{\mathcal{A}}_{0}^{\text{t}}+\frac{e}{2\hbar}(\bm{B}\times\bm{\mathcal{A}}_{0}^{\text{t}}\cdot\bm{\partial}_{\bm{p}})\bm{\mathcal{A}}_{0}^{\text{t}}+\frac{e}{2\hbar}\bm{\Omega}_{0}^{\text{t}}\times(\bm{B}\times\bm{\mathcal{A}}_{0}^{\text{t}})$ according to Eq. (27). As a result, the first order change in electric polarization reads

[TABLE]

The first term is the Abelian Chern-Simons 3-form, which yields the topological part of the orbital magnetoelectric coefficient. The second term yields the electric polarizability and the cross-gap part of the orbital magnetoelectric polarizability, consistent with the calculation using the linear response theory Essin et al. (2010). Such consistency confirms the validity of the first order correction to the Berry phase.

The modified band energy in Eq. (26) should contain the correction up to second order, i.e. it can be put in the following form: $\tilde{\varepsilon}_{0}=\varepsilon_{0}-\bm{B}\cdot\bm{m}+e\bm{E}\cdot\bm{r}_{c}+\varepsilon_{(2)}$. The second order correction describes the change in the electric dipole and magnetic moment of the wave packet in response to external fields. It reads Gao et al. (2015)

[TABLE]

where $\alpha_{k\ell}=\partial_{k}\partial_{\ell}\varepsilon_{0}/\hbar^{2}$ is the inverse effective mass tensor, $m$ is the free electron mass, $G_{0n}=(G_{n0})^{\star}$, and

[TABLE]

In the energy correction, the first two terms are purely due to the horizontal mixing and hence are geometrical corrections involving Berry curvature and quantum metric. The third term is due to the cross effect of the horizontal and vertical mixing. $\bm{P}_{\text{E}}$ is called energy polarization as it accounts for the energy dipole in the momentum space. The fourth term is purely due to the horizontal mixing and has the conventional form of the perturbation correction to the energy. The fifth term shifts the momentum argument of $\varepsilon_{0}$. The last term is the expectation of the second order correction to the local Hamiltonian.

We comment that if one starts from an effective tight-binding Hamiltonian, the last term of Eq. (II.4) will change. Instead of the free electron mass, one has to use the Hessian operator $\hat{\Gamma}_{ij}=\partial^{2}\hat{H}_{0}/\partial\hat{p}_{i}\partial\hat{p}_{j}$. Then the last term in Eq. (II.4) should be replaced by the following two terms

[TABLE]

Here $(\hat{\Gamma}_{ij})_{mn}=\langle u_{m}|\hat{\Gamma}_{ij}|u_{n}\rangle$.

II.5 Landau level quantization and magnetic susceptibility

If one starts from a continuum tight-binding Hamiltonian, electronic states will be quantized into Landau levels under magnetic field. Therefore, as long as the semiclassical theory can successfully reproduce the Landau level quantization up to second order, it is enough to account for various response functions at second order in magnetic field.

There are various quantization rules for Landau levels, such as the Onsager’s rule Onsager (1952), the Maslov canonical operator method Reijnders et al. (2013), the Gutzwiller trace formula Gutzwiller (1971), and so on. For semiclassical theory, the most relevant one is the Onsager’s rule Onsager (1952), which reads

[TABLE]

where $S(\varepsilon_{n})$ is the area in the momentum space enclosed by an equal-energy contour. Later, this quantization rule is generalized up to first order to take account of corrections from Berry phase and magnetic moment Wilkinson (1984); Rammal (1985); Mikitik and Sharlai (1999); Chang and Niu (1996); Carmier and Ullmo (2008); Xiao et al. (2010); Fuchs et al. (2010).

To further generalize the Onsager’s rule, one observes that the area $S(\varepsilon_{n})$ is proportional to the total electron density below $\varepsilon_{n}$. Since the area $S$ as a function of energy is a smooth function, such electron density should be understood as the smooth semiclassical electron density $\rho_{\text{semi}}$. Then it has been proved that the Onsager’s rule can be replaced by the following density quantization rule Gao and Niu (2017)

[TABLE]

This suggests that the smooth semiclassical electron density always intersects with the true stepwise electron density at the half-filling points, as shown in Fig. 2.

In general $\rho_{\text{semi}}$ can be expanded near $B=0$ in power series of $B$. Up to second order we have

[TABLE]

where $M$ and $\chi$ are the magnetization and magnetic susceptibility respectively, evaluated at zero magnetic field and zero temperature. It is interesting to note that Eq. (41) indicates the nonlinearity in the Landau level fan diagram when plotted as $B$ against $1/n$ due to the appearance of the susceptibility. Such nonlinearity is consistent with experiments Taskin et al. (2010); Analytis et al. (2010); Ren et al. (2010); Sacépé et al. (2011); Br ne et al. (2011); Xiu et al. (2011); Xiang et al. (2015). Based on Eq. (41), as long as the semiclassical theory can yield the correct susceptibility, it can give correct second order responses in magnetic field.

To evaluate the magnetic susceptibility, we need to calculate the free energy up to second order. In the semiclassical framework, the free energy can be expressed as follows

[TABLE]

Here the second term $g_{\text{L}}$ is the Peierls-Landau magnetic free energy: $g_{\text{L}}=-(e^{2}f_{0}^{\prime}/48\hbar^{2})B_{s}B_{t}\epsilon_{sik}\epsilon_{tj\ell}\alpha_{ij}\alpha_{k\ell}\,,$ where $f_{0}^{\prime}$ is the energy derivative of the Fermi distribution function $f_{0}$. For isotropic bands, the effective mass tensor $\alpha$ is diagonal, and $g_{\text{L}}$ will reduce to its familiar form Ashcroft and Mermin (1976). This term originates from higher order corrections to the replacement of the quantum mechanical commutator with the classical Poisson bracket Blount (1962b).

After collecting terms at second order in $B$, one can find the following free energy $F=\int g^{\prime\prime}(d\bm{k}/8\pi^{3})$ with

[TABLE]

The magnetic susceptibility can be readily obtained through its definition

[TABLE]

The first two terms in Eq. (II.5) are Peierls-Landau diamagnetic and Pauli paramagnetic contributions. The physical meaning of the other terms in Eq. (II.5) can be illustrated by re-expressing it in the Wannier function basis and then taking the atomic insulator limit. In particular, the two terms in the second line reduce to the familiar form of the Van-Vleck paramagnetic free energy and Langevin diamagnetic free energy in atomic physics, respectively. The last term is purely geometrical as it solely comes from the horizontal mixing. The energy polarization and geometrical magnetic free energy approach zero in the atomic limit, indicating that they have no analog in atomic physics and are novel terms in crystals.

The susceptibility obtained from Eq. (II.5) in the semiclassical framework is consistent with the result derived using other methods. For example, the susceptibility has been calculated in the honeycomb lattice using Hofstadter spectrum or linear response theory Gómez-Santos and Stauber (2011); Raoux et al. (2014). Eq. (II.5) yields the exactly same result as shown in Fig. 3. In Ref. Ogata and Fukuyama (2015), it is also proved that Eq. (II.5) is consistent with the susceptibility from linear response theory 222In Ref. Ogata and Fukuyama (2015), it is found that except the energy polarization contribution, all the other terms are consistent. For the energy polarization contribution, Ref. Gao et al. (2015) contains a typo. When inserting the second order energy in Eq. (42), the energy polarization in Ref. Gao et al. (2015) has an additional $1/4$ factor by mistake. After removing such factor as given in Eq. (II.5), the energy polarization has the same expression with Eq. (2.31) in Ref. Ogata and Fukuyama (2015). This confirms the validity of the semiclassical dynamics under magnetic field.

II.6 Evaluating positional shift and energy correction in tight-binding models

The positional shift and energy correction can be implemented in first-principles codes. For this purpose, besides the energy spectrum, one needs three additional matrix elements: the velocity matrix element $\bm{v}_{mn}=\langle m|\hat{\bm{v}}|n\rangle$, the Hessian matrix element $(\Gamma_{ij})_{mn}=\langle u_{m}|\hat{\Gamma}_{ij}|u_{n}\rangle$, and the spin Pauli matrix element $\bm{\sigma}_{mn}=\langle m|\hat{\bm{\sigma}}|n\rangle$, all of which can be evaluated in first-principles codes in principle Marzari and Vanderbilt (1997); Souza et al. (2001); Yates et al. (2007); Marzari et al. (2012). In the following, we sketch the process of re-expressing the positional shift and energy correction in terms of the energy spectrum and these three matrix elements.

We first examine the positional shift. The interband part of the Berry connection is related to the interband velocity element:

[TABLE]

Using this identity, we manipulate the interband mixing element $G_{n0}$

[TABLE]

Then the first term in the positional shift can be put in the desired form

[TABLE]

For the remaining term we note that

[TABLE]

The only unknown quantity has the form $\partial_{i}(v_{j})_{mn}$, which can be manipulated as follows

[TABLE]

The first three terms are in the desired form. The last term contains the intraband Berry connection which is gauge-dependent. Therefore, it cannot be effectively evaluated in the first-principles codes. However, it can be shown that it does not contribute to Eq. (II.6) and hence need not to be evaluated. Equations (46) to (II.6) are enough to transform the positional shift in the desired form.

We now examine the energy correction. Using the following identities

[TABLE]

together with Eq. (45) to (II.6), it is straightforward to put all terms except the energy polarization contribution in Eq. (II.4) in the desired form. We comment that for the tight-binding Hamiltonian, one may have to substitute the last term in Eq. (II.4) with Eq. (II.4) which has realistic Hessian matrix. In this case, we have

[TABLE]

To evaluate the last two terms, one notes that

[TABLE]

Using Eq. (55) and Eq. (II.6), one can put the second term in $\varepsilon_{\text{H}}$ in the desired form plus extra terms from intraband Berry connection. After inserting Eq. (II.6) into $\varepsilon_{\text{H}}$, the first term in the second line of Eq. (II.6) does not contribute due to the anti-symmetrization of the spatial indices. The last term in Eq. (II.6) will cancel the extra gauge-dependent terms in Eq. (55). This completes the manipulation of $\varepsilon_{\text{H}}$.

The manipulation of energy polarization contribution in Eq. (II.4) can be done as follows:

[TABLE]

where

[TABLE]

Again, terms contain intraband Berry connection can be ignored as they will cancel each other in the end.

We comment that if one starts from the Schr$\ddot{\rm o}$dinger Hamiltonian, the Hessian matrix can be directly evaluated without using wave functions in first-principles codes. For Schrodinger Hamiltonian, $\hat{\Gamma}_{ij}=1/m\delta_{ij}$. Therefore, one has $(\Gamma_{ij})_{mn}=1/m\delta_{ij}\delta_{mn}$. Then for the semiclassical dynamics, one only needs to evaluate the velocity and spin matrix elements using output from first-principles codes.

As a concrete example, we consider the following tight-binding Hamiltonian describing a 3D TI in the Bi2Se3 family on a simple cubic lattice Fu and Berg (2010); Qi and Zhang (2011); Vazifeh and Franz (2013)

[TABLE]

where $\bm{\sigma}$ and $\bm{s}$ are Pauli matrices in orbital and spin space respectively, and $M_{k}=\epsilon-2t(\cos k_{x}+\cos k_{y}+\cos k_{z})$. At $\epsilon=6t$, the above Hamiltonian is at the phase transition point between trivial phase and topological insulator phase, and supports a Dirac point at the $\Gamma$ point. The Dirac point can split into Weyl point if the following additional term is added to the Hamiltonian

[TABLE]

These two terms break inversion and time reversal symmetry, respectively. The resulting Weyl semimetal phase is displayed in Fig. 4.

Here we choose the semiconductor phase and calculate the first and second order correction to the band energy, as shown in Fig. 5. It can be seen that both first and second order correction vary drastically near the band crossing points and the small band gap. In fact, both corrections diverge at the $\Gamma$ point due to the band crossing. Across the small global band gap as shown in Fig. 4, both first and second order corrections change sign. It is also interesting to note that near the $\Gamma$ point, the lowest two bands have similar first order corrections, but are opposite to those of the remaining two bands. This property is consistent with the calculation of magnetic moment for the low-energy Dirac model. In comparison, the second order correction have different properties. The lowest and the second highest band have similar second order corrections, opposite to those of the remaining two bands. This can be most easily seen from Eq. (II.4). The Berry curvature, effective mass tensor and energy gap reverse sign for the two bands forming a Dirac cone, leading to the sign change in their second order energy correction.

III Nonlinear charge current

With the semiclassical theory up to second order at hand, we can derive the nonlinear charge current up to third order, i.e. conductivity up to second order in fields. We focus on the semiclassical regime with weak electromagnetic fields. There is one important issue that we want to discuss before the detailed derivation of nonlinear currents. In the previous section, the semiclassical theory is generally derived for a single Bloch band that is well separated from all the other bands. In reality, this condition is hardly met as band-crossing points are generally present. In this case, the previous semiclassical theory can still be used to derive the conductivity as long as the Fermi surface is not close to those band-crossing points. The reason is as follows. As a perturbation theory, the semiclassical theory may fail when the band gap is small, such as near the band crossing points. However, the conductivity is a Fermi surface property. As a result, as long as the band-crossing points are deep inside the Fermi sea, the perturbation theory still works near the Fermi surface.

In this section, we will first establish the general theory of the charge current beyond the linear order. Then we will discuss several important examples, including the response of the anomalous Hall conductivity to electromagnetic fields, linear magnetoresistance in time-reversal-broken materials, chiral anomaly in Weyl semimetals and two different mechanisms for the negative longitudinal magnetoresistance: the intrinsic quadratic magnetoresistance and the current jetting.

III.1 General theory of nonlinear charge currents

The current in the semiclassical theory reads

[TABLE]

where $f$ is the electron distribution function, which is the Fermi function in equilibrium. Here and hereafter, we drop the subscript $c$ in the center of mass position $\bm{r}_{c}$ and momentum $\bm{k}_{c}$ for simplicity. From the semiclassical dynamics in Eq. (30) and (31), one can show that

[TABLE]

where $\tilde{\bm{v}}=\partial\tilde{\varepsilon}/\hbar\partial\bm{k}$ is the modified band velocity. Here and hereafter, we ignore the band index [math] in relevant intraband quantities for simplicity and the subscript [math] only has the meaning of zeroth order.

The remaining factor in the current is the distribution function $f$, which is typically solved from the Boltzmann equation. If the sample is homogeneous and reaches a steady state under external fields, under the relaxation time approximation, the Boltzmann equation reads

[TABLE]

Here the argument of equilibrium distribution $f_{0}$ is the modified band energy $\tilde{\varepsilon}$. Solving the Boltzmann equation perturbatively, one has Pal and Maslov (2010)

[TABLE]

Equations (61), (62) and (64) are enough to derive the current except the anomalous Hall current up to third order of the external fields, i.e. the conductivity up to second order. Other than the anomalous Hall current, the equilibrium distribution $f_{0}$ alone does not contribute to the current. Therefore, at least the first order correction to $f_{0}$ is required, which automatically contains an external field due to $\dot{\bm{k}}$. As a result, one only needs to evaluate the factor $\mathcal{D}\dot{\bm{r}}$ up to second order, which is given in the semiclassical dynamics. For the anomalous Hall conductivity, the semiclassical theory can only yields its first order correction, as only the first order correction to Berry curvature is obtained in the semiclassical theory.

To explicitly write out the solution in Eq. (64) up to third order in fields, we note that

[TABLE]

Therefore, we write the distribution function as follows

[TABLE]

where the subscript $1$, $2$, and $3$ stand for different orders in $\tau$, and

[TABLE]

[TABLE]

There is another contribution to the nonlinear conductivity. The system can reach a thermodynamic equilibrium together with a uniform magnetic field. During such process, if the sample remains charge neutral, the carrier density is fixed. As the magnetic field modifies the band energy, the chemical potential has to vary accordingly.

Without loss of generality, we assume $\bm{B}$ is along the $z$ direction, i.e. $B=B_{z}$. At first order in $B$, the change in chemical potential can be fixed using the following expansion of the electron density $n$

[TABLE]

Note that $\left.\frac{\partial n}{\partial\mu}\right|_{T,B=0}=g(\mu)$, which is simply the density of states at the chemical potential for the unperturbed band. Moreover, $\left.\frac{\partial n}{\partial B}\right|_{\mu,T}$ can be connected to the magnetization through the Maxwell equation:

[TABLE]

We comment that for insulators at zero temperature, $-\frac{\partial M_{z}}{\partial\mu}\rightarrow\sigma_{xy}$, and Eq. (71) reduces to the familiar Streda formula. Combining Eq. (70) and (71), we obtain

[TABLE]

At second order in $B$, the expansion of $n$ reads

[TABLE]

$\mu_{(2)}$ can be solved by setting the above equation equal to zero. The result can be simplified using Eq. (71) and the definition of the isothermal magnetic susceptibility $\chi_{zz}=\partial M_{z}/\partial B$, yielding

[TABLE]

These change in chemical potential will also contribute to the nonlinear conductivity:

[TABLE]

where $f_{0}^{\prime\prime}$ and $f_{0}^{\prime\prime\prime}$ are second and third order energy derivatives of $f_{0}$.

Equation (61) with Eq. (62) and (67)-(69), combined with Eq. (72), (74), and (III.1) yields a complete description of the nonlinear conductivity up to second order in the framework of the semiclassical theory. This procedure gives a complete account for the drifting part in the Boltzmann equation. However, it ignores the dependence of the collision integral on electromagnetic fields. Nevertheless, the semiclassical theory can yield a qualitatively valid result and offers a fresh perspective of how the band properties beyond the spectrum affects the conductivity.

Finally, we comment that to get the corresponding change in the resistivity, one can invert the conductivity tensor. For example, if change in conductivity $\delta\sigma_{xx}$ is obtained, the corresponding change in resistivity is

[TABLE]

provided that the Hall conductivity $\sigma_{xy}$ is much smaller than the longitudinal conductivity $\sigma_{xx}$.

III.2 Nonlinear anomalous Hall conductivity

The anomalous Hall effect refers to a Hall-type current in ferromagnets solely driven by an electric field. It is a topic under extensive studies Nagaosa et al. (2010). In the past, there are three mechanisms identified: the intrinsic contribution, the skew-scattering contribution, and the side-jump contribution. Here we focus on the intrinsic contribution, which is due to the anomalous velocity in Eq. (16) Xiao et al. (2010). From the perspective of the semiclassical theory, the Berry curvature acts as the magnetic field in the momentum space and bends the electron trajectory likewise, leading to a Hall-type current. The anomalous Hall effect has also been studied in noncolinear antiferromagnets, in which the net spin magnetization vanishes but the orbital magnetization is still present Ohgushi et al. (2000); Shindou and Nagaosa (2001); Taillefumier et al. (2006); Kalitsov et al. (2009); Takatsu et al. (2010); Udagawa and Moessner (2013); Chen et al. (2014); Suzuki et al. (2017); Guo and Wang (2017).

From symmetry consideration, the anomalous Hall effect requires broken time-reversal symmetry. In other words, it requires either a net spin magnetization as in ferromagnets or orbital magnetization as in noncolinear antiferromagnets. It is forbidden if time reversal symmetry is present. However, the nonlinear Hall effect may still be present, due to the manipulation of the magnetization through electromagnetic fields.

In systems with both time reversal and inversion symmetry, a net magnetization can be induced through the magnetic susceptibility by a magnetic field. This magnetization can then lead to a Hall-type current. The analytical expression for this correction can be obtained by plugging the anomalous velocity (second term in Eq. (62)) into the current in Eq. (61). The result reads Gao et al. (2014)

[TABLE]

where $\bm{\mathcal{A}}^{\prime}(\bm{B})$ stands for the part of $\bm{\mathcal{A}}^{\prime}$ solely dependent on $\bm{B}$ and $\bm{m}$ is the spin plus orbital magnetic moment. It can be explicitly checked that compared with the susceptibility in Eq. (II.5), the second term in Eq. (77) is part of the geometric contribution to the magnetic susceptibility. Moreover, the first term is also part of the Van-Vleck and energy polarization contributions to the susceptibility. This current yields the conductivity $\sigma_{(ij,k),0}$ in Eq. (3).

This current has exactly the same dependence on electromagnetic fields with the ordinary Hall current. However, they have different origins. The ordinary Hall current is due to the Lorentz force and has the form of $\omega_{c}\tau$ with $\omega_{c}$ being the cyclotron frequency. In contrast, the current in Eq. (77) requires a nontrivial structure in the momentum space and does not involve the relaxation time. In fact, the ratio of the resistivity $\rho_{xy}^{\prime}$ from Eq. (77) and the ordinary Hall resistivity $\rho_{xy}^{\text{ord}}$ can be put in the following form

[TABLE]

where the first factor is universal and $S(\mu)$ is a model-dependent factor but independent of the scattering process. The universal scaling factor means that $\rho^{\prime}_{xy}$ will dominate the ordinary Hall effect when $\rho_{xx}$ is large, i.e. in dirty metals/semiconductors with relatively small relaxation time. In a typical Hall-bar measurement set-up, both $\rho_{xy}^{\prime}$ and $\rho_{xy}^{\text{ord}}$ will contribute to the Hall current. Therefore, the total Hall resistivity should be dependent on the relaxation time. To differentiate one contribution from the other, one should change the universal scaling factor through temperature, film thickness, or doping, and measure the scaling behaviour based on Eq. (78).

In systems that simultaneously break the time reversal and inversion symmetry but preserve the combined symmetry, the anomalous Hall current also vanishes. This symmetry is exactly the magnetoelectric symmetry which allows a magnetization to be induced by an electric field through the magnetoelectric coefficient Landau et al. (1984). This magnetization can then lead to a Hall-type current. The analytical expression can be obtained in a similar way with the magnetic-field-correction to anomalous Hall effect. The result reads Gao et al. (2014)

[TABLE]

where $\bm{\mathcal{A}}^{\prime}(\bm{E})$ stands for the part of $\bm{\mathcal{A}}^{\prime}$ solely dependent on $\bm{E}$. This current yields the conductivity $\sigma_{ijk,0}$ in Eq. (2). The connection between this Hall-type current and the magnetoelectric effect can be mostly seen for a two-band model, as shown in Fig. 6.

In noncentrosymmetric but time-reversal-invariant materials, a net magnetization can be induced in the nonequilibrium steady state through the Edelstein effect. Such magnetization can further lead to a Hall-type current. The analytical expression can be obtained by plugging Eq. (62) and (67) into Eq. (61) and keeping terms at second order in $\bm{E}$, and reads Sodemann and Fu (2015)

[TABLE]

This current yields the conductivity $\sigma_{ijk,1}$ in Eq. (2). Interestingly, the requirement of the mirror symmetry in Fig. 6(b) also works for this current. But one should keep in mind that for the current in Eq. (79), as the time reversal symmetry is broken, the inherent spin texture in the sample is also subjected to the mirror operation. The observation of this nonlinear Hall current has been reported in several recent experiments Xu et al. (2018); Kang et al. (2018); Ma et al. (2018). We comment that in Eq. (80), $\bm{v}\bm{\Omega}$ together constitutes a pseudotensor which represents the first order moment of the Berry curvature in the momentum space. Hence it is referred to as the Berry curvature dipole. This Berry curvature dipole, together with the magnetic moment dipole also play essential roles in the study of the natural optical activity Landau et al. (1984); Malashevich and Souza (2010); Orenstein and Moore (2013); Ma and Pesin (2015); Zhong et al. (2015, 2016).

Finally, using the quantum kinetic equations for the evolution of the density matrix under both electric field and disorders, it has been found that in noncentrosymmetric materials, besides the Berry curvature dipole contribution, there are also nonlinear side jump and skew scattering mechanisms responsible for the nonlinear anomalous Hall effect at second order of electric field Nandy and Sodemann (2019). It can compete with the nonlinear anomalous Hall effect due to Berry curvature dipole.

III.3 Linear Magnetoresistance

Linear magnetoresistance, referring to the linear dependence of the resistance on magnetic field, was first observed in simple metals Reitz and Overhauser (1967); Penz and Bowers (1968); Jones (1969) and under continuous studies afterwards. There is a surge of interest in recent years, as it appears in various novel materials, such as multilayer graphene Friedman et al. (2010), topological insulators Qu et al. (2010); Tang et al. (2011); Wang et al. (2012); Tian et al. (2014), and Weyl/Dirac semimetals He et al. (2014); Liang et al. (2015); Feng et al. (2015a); Novak et al. (2015); Narayanan et al. (2015); Huang et al. (2015). It has been shown that the linear magnetoresistance can arise in the ultra-quantum regime for Dirac systems when only the lowest Landau level is partially filled Abrikosov (1998, 2000); Wang and Lei (2012), as well as in inhomogeneous samples with mobility fluctuation Herring (1960); Parish and Littlewood (2003); Porter and Marrows (2012); Kozlova et al. (2012). In the semiclassical regime, linear magnetoresistance is subject to a symmetry requirement. Generally speaking Onsager relation requires that $\sigma_{xx}(\bm{B})=\sigma_{xx}(-\bm{B})$, which forbids odd power dependence of the conductivity on magnetic field. However, if the time reversal symmetry is broken and a magnetization is present, Onsager relation becomes $\sigma_{xx}(\bm{B},\bm{M})=\sigma_{xx}(-\bm{B},-\bm{M})$, which can be satisfied in principle for terms with odd powers of $B$.

To derive the linear magnetoresistance, or equivalently the linear magnetoconductivity in the semiclassical regime, one collects terms linear in $B$ in the current in Eq. (61). To simplify the result, we first consider the transverse magnetoresistance and assume that the electric field is along $\hat{x}$ direction and $\bm{B}$ is along $\hat{z}$ direction. This configuration will yield the following conductivity (here we use the definition in Eq. (3)) Chen et al. (2015)

[TABLE]

For longitudinal magnetoresistance, we assume that the electric field is parallel to the magnetic field, i.e. it is along $\hat{z}$ direction. Then the magnetoconductivity reads

[TABLE]

Both $\sigma_{(xx,z),1}$ and $\sigma_{(zz,z),1}$ are clearly violations to the Kohler’s rule Pippard (1989). In other words, the magnetic field affects the conductivity not through the Lorentz force, but through the two extra corrections in the first order semiclassical equations of motion. On one hand, it couples to the $k$-dependent magnetic moment to modifies the band energy and hence band velocity. On the other hand, it couples to the Berry curvature and changes the phase of space density of states, or equivalently, the carrier density. In the longitudinal magnetoconductivity, there is an additional term with the factor $\bm{v}\cdot\bm{\Omega}$. This actually measures the flux strength of the momentum space magnetic field $\bm{\Omega}$. This unique contribution to the longitudinal magnetoresistance is due to the coupling between the real and momentum space lorentz force, i.e. the magnetic field along $\hat{z}$ can first bend the velocity along $\hat{x}$ (or $\hat{y}$) to the direction $\hat{y}$ (or $\hat{x}$), which is further bent to the $\hat{z}$ direction due to the momentum space magnetic field $\Omega_{x}$ (or $\Omega_{y}$). Finally, we comment that both $\sigma_{(xx,z),1}$ and $\sigma_{(zz,z),1}$ are linear in magnetic moment and Berry curvature, which are odd in time reversal operation. As a result, these two conductivities only appear in materials that break time reversal symmetry.

The fact that $\sigma_{(xx,z),2}$ and $\sigma_{(zz,z),2}$ vanish identically is consistent with the Onsager’s relation, as formally they do not need to break the time reversal symmetry according to Table. 1. We comment that if a detailed consideration of scattering process is performed, there is a contribution to $\sigma_{(xx,z)}$ that is of second order in $\tau$, due to the skew scattering process Chen et al. (2015). Therefore, it requires a nontrivial antisymmetric part of the scattering probability $W_{\bm{k}\bm{k}^{\prime}}$ and hence the breaking of time reversal symmetry.

When the chemical potential approaches the hot spot of the Berry curvature and the magnetic moment, the latter can have a large contribution to the linear magnetoresistance, as shown in Fig. 7. As the chemical potential approaches the band edge, the Berry curvature and magnetic moment increases and so does the linear magnetoresistance. However, the semiclassical theory is only valid when $(\mu-\Delta)\tau/\hbar>1$. Otherwise, the scattering from different impurities will be correlated and a fully quantum mechanical treatment should be used.

III.4 Negative longitudinal magnetoresistance induced by chiral anomaly

It is well known that the chemical potential difference in real space can lead to a current. Its counterpart in momentum space, i.e. the current due to the chemical potential difference in momentum space, has completely different origins and attracts great attention recently as it manifests in topological semimetals.

To understand this effect, we note that from the semiclassical theory, the Berry curvature is the analog of the magnetic field in the momentum space. The flux of the Berry curvature, defined as $(1/2\pi)\oint d\bm{S}\cdot\bm{\Omega}$, can be nonzero through a closed surface, indicating the appearance of a monopole charge in the momentum space. This is most easily realized by a pair of Weyl points

[TABLE]

where $\lambda_{i}=\pm 1$ stands for a positive or negative charge, or equivalently, the chirality of the Weyl node. For such pair of Weyl nodes, there is a current driven by the magnetic field. This is the chiral magnetic current, derived using the last term in Eq. (62)

[TABLE]

where $\mu_{i}$ is the chemical potential for each Weyl node. This current has to vanish in equilibrium Vazifeh and Franz (2013), even when the two Weyl nodes are shifted in energy and the relative chemical potential is different.

However, away from equilibrium, the chemical potential for different valleys can respond to external fields in different ways, leading to an inhomogeneous chemical potential distribution in the momentum space and hence a net current through $J_{\text{CME}}$ Son and Spivak (2013). To see this, we first write down the Boltzmann equation:

[TABLE]

where $W_{\bm{k}\bm{k}^{\prime}}$ is the scattering probability. Using the semiclassical equations of motion in Eq. (16) and (17), one finds that

[TABLE]

Therefore, the Boltzmann equation becomes

[TABLE]

In multivalley systems such as a pair of Weyl nodes, a special solution to the above equation may emerge. To see this, we note that usually the transport current is determined by the fastest scattering rate which is the intravalley scattering. For the intravalley scattering with $f_{0}(\bm{k})$ and $f_{0}(\bm{k}^{\prime})$ from the same type of Weyl node, one can multiply Eq. (III.4) by $\mathcal{D}(\bm{k})$ and integrate it around the Weyl node. The result is

[TABLE]

However, for a Weyl node, although the integration of the first two terms vanish, the integration of the third term is proportional to the charge of the Weyl node and hence does not vanish. It indicates that charge can be pumped from the positive Weyl node to the negative Weyl node. This contradiction shows that the intervalley scattering is indispensable in solving the above Boltzmann equation.

To explicitly derive this solution, one can take $f_{0}(\bm{k})$ to be from the positive Weyl node and $f_{0}(\bm{k}^{\prime})$ from the negative Weyl node. Then under the relaxation time approximation, Eq. (90) for the positive Weyl node becomes

[TABLE]

where $\Delta N_{+}$ and $\Delta N_{-}$ is the density of particles near the Fermi surface, and $\tau_{\text{int}}$ is the intervalley scattering time. Eq. (91) means that the pumping of particles from one Weyl node to another through drifting has to be balanced by the diffusive collision through the intervalley scattering process. It leads to an inhomogeneous chemical potential distribution in the momentum space

[TABLE]

where $g(\mu)=\mu^{2}/(2\pi^{2}\hbar^{3}v^{3})$ is the density of states at the Fermi surface and $\tau_{\text{int}}$ is the intervalley scattering time. This chemical potential difference for Weyl nodes with opposite chiralities is the chiral anomaly effect Nielsen and Ninomiya (1983); Aji (2012); Son and Spivak (2013); Burkov (2015a); Armitage et al. (2018).

This imbalance in the chemical potential can lead to a current through Eq. (86)

[TABLE]

This current is only present when $\bm{E}$ is parallel to $\bm{B}$. The corresponding conductivity is quadratic in magnetic field and increases with it, indicating a negative longitudinal magnetoresistance. It corresponds to $\sigma_{(ij,kk),1}$ in Table. 1 and does not need to break time reversal or inversion symmetry.

Besides the semiclassical theory, the negative magnetoresistance has also been confirmed using chiral kinetic theory Stephanov and Yin (2012); Son and Yamamoto (2013); Spivak and Andreev (2016); Hidaka et al. (2017); Sekine et al. (2017) and exact quantized Landau levels Burkov (2014). In the semiclassical regime, the negative magnetoresistance for one pair of Weyl nodes is also extended in other situations. It has been shown that samples with Dirac points also have negative magnetoresistance, provided that the chiral charge is approximately conserved and the chiral relaxation time is large Burkov (2015b). Moreover, when the Dirac point opens a gap due to the suppression of the lattice symmetry, the chiral relaxation time is still present, so is the negative magnetoresistance Andreev and Spivak (2018). Experimentally, the negative longitudinal magnetoresistance in Weyl and Dirac semimetals has been reported and interpreted as the chiral magnetic effect Huang et al. (2015); Shekhar et al. (2015); Yang et al. (2015); Wang et al. (2016); Zhang et al. (2016); Kim et al. (2013); Xiong et al. (2015); Feng et al. (2015b); Li et al. (2015); Liang et al. (2015); Li et al. (2016); Zhang et al. (2017), but similar signals could also arise from two different mechanisms discussed in the following two sections.

Finally, we comment that the Weyl nodes is gapless and near those band-crossing points, the semiclassical theory may fail. However, we can set the Fermi surface to be well above the Weyl point. The flux through the Fermi surface does not change with the Fermi energy. As a result, the discussion using the semiclassical theory still works.

III.5 Intrinsic quadratic longitudinal magnetoresistance

The magnetoconductivity due to chiral anomaly and chiral magnetic effect is one contribution to $\sigma_{(ij,kk)}$ in Eq. (4). It is not a complete description of $\sigma_{(ij,kk)}$. To fully account for the drifting contribution to $\sigma_{(ij,kk)}$, one has to use the semiclassical theory up to second order, as presented in Sect. III (A). Specially, one can plug Eq. (67) and (62) into Eq. (61) and collect terms up to second order in magnetic field. The current responsible for $\sigma_{(ij,kk),1}$ is complicated but it can be put in several groups:

[TABLE]

The first contribution has the form

[TABLE]

This term is purely due to the anomalous velocity as both $\bm{E}\times\bm{\Omega}$ and $(\bm{v}\times\bm{B})\times\bm{\Omega}$ are parts of the anomalous velocity. This current is particularly interesting. If the electromagnetic fields are along the same direction, which can be assume to be the $\hat{z}$ direction without loss of generality, this current becomes

[TABLE]

This will always lead to a negative longitudinal magnetoresistance.

The second contribution has the form

[TABLE]

This is solely due to the Zeeman coupling between magnetic field and the magnetic moment, containing both the spin and orbital part.

The third contribution has the form

[TABLE]

This current is due to the coupling between the anomalous velocity and the Zeeman energy correction.

The fourth contribution has the form

[TABLE]

The first line is due to the second order correction to the band energy and the remaining terms are due to the first order correction to the Berry curvature, i.e. $\bm{\Omega}^{\prime}$.

From the expressions of $\bm{J}_{1}$ to $\bm{J}_{4}$, we find that $\bm{J}_{1}$ is purely due to the orbital motion of Bloch electrons, while $\bm{J}_{2}$ to $\bm{J}_{4}$ can contain both spin and orbital contribution. The orbital part has been calculated in Ref. Gao et al. (2017) and the spin part has been calculated in Ref. Dai et al. (2017). Here $\bm{J}_{1}$ to $\bm{J}_{4}$ are complete expressions that contain the spin part, the orbital part, and the spin-orbital part contribution.

The last contribution is due to the correction to the chemical potential:

[TABLE]

In nonmagnetic metals, $\sigma_{(ij,kk),3}$ is also nonzero. The corresponding current reads

[TABLE]

This current requires a peculiar Fermi surface such that the integration does not vanish Pal and Maslov (2010).

The current $\bm{J}$ and $\bm{J}^{\prime}$ complete the description of the quadratic magnetoresistance. In $\bm{J}^{\prime}$ the magnetic field always couples to the relaxation time. Therefore, the corresponding magnetoresistance follows Kohler’s rule Pippard (1989). In comparison, $\bm{J}$ violates Kohler’s rule and its corresponding magnetoresistance does not depend on the relaxation time. Moreover, $\bm{J}$ does not need a highly anisotropic Fermi surface. As long as the Berry curvature, magnetic moment, and magnetic susceptibility does not vanish at the same time, $\bm{J}$ will be nonzero. This intrinsic contribution to magnetoresistance is also confirmed by a fully quantum mechanical treatment Wang et al. (2018).

The current $\bm{J}$ can yield a negative longitudinal magnetoresistance. In fact, for the Dirac Hamiltonian in Eq. (85), the ratio of the quadratic longitudinal magnetoresistance $\delta\rho$ from $\bm{J}$ to that from chiral anomaly reads

[TABLE]

The global positive sign means that $\delta\rho$ is negative, like $\delta\rho_{\text{CA}}$. Usually $\tau\ll\tau_{\text{int}}$. Therefore, in Weyl semimetals, the negative magnetoresistance should be dominated by the contribution from the chiral anomaly in the semiclassical regime. However, the negative magnetoresistance from $\bm{J}$ can persist in metals without any Dirac/Weyl point, as shown in Fig. 8. Panel (a) suggests that the longitudinal magnetoresistance is negative while the transverse one is positive. This causes the sign change of the effective magnetoresistance as the angle between $\bm{E}$ and $\bm{B}$ continuously changes from [math] to $\pi/2$, as shown in Panel (b).

It is interesting to further check the behaviour of the magnetoconductivity in a lattice model. We use the model Hamiltonian in Eq. (59) and (60) and choose the parameters such that it is in the semiconductor phase as shown in Fig. 4 (d). We then calculate the quadratic magnetoconductivity for the transverse and longitudinal configuration, using the expression $\bm{J}_{1}$ to $\bm{J}_{5}$. Both spin and orbital part of the magnetic moment is considered. The result is presented in Fig. 9. It can be found that both magnetoconductivity show strong dependence on the chemical potential. Consistent with the result in the low-energy model, the magnetoresistance will be strongly enhanced when chemical potential is near the band edge, as both the anomalous velocity and energy correction are enhanced. Moreover, the longitudinal magnetoconductivity is generally positive, indicating a negative longitudinal magnetoresistance, while the transverse magnetoconductivity is negative. Near the band edge, the longitudinal magnetoconductivity is larger than the transverse one, indicating a sizeable transition angle $\theta$ from negative to positive magnetoresistance, as the relative angle between $\bm{E}$ and $\bm{B}$ changes from [math] to $\pi/2$. This is similar to the negative magnetoresistance induced by the chiral anomaly.

III.6 Current Jetting

Beyond the Boltzmann transport theory, the magnetoresistance is also greatly affected when the current distribution in the sample is inhomogeneous. For example, for a matchstick sample with two leads for current input and output, as shown in Fig. 10, the current in the sample can be highly inhomogeneous provided that the transverse magnetoresistance is sufficiently large Yoshida (1980); Pippard (1989). This is the current jetting problem.

The current jetting problem can be solved through the charge continuity equation: $\partial\rho/\partial t+\bm{\nabla}\cdot\bm{J}=0$. In the steady state $\partial\rho/\partial t=0$. Therefore, in the sample, the following equation for the electrostatic potential is satisfied

[TABLE]

Since the Hall conductivity is antisymmetric, it does not affect the above equation. However, it can affect the boundary condition. If the Hall conductivity can be ignored, and the longitudinal conductivity is determined by the Drude formula: $\sigma_{xx}=\sigma_{yy}=\sigma_{0}/(1+\omega_{c}^{2}\tau^{2})$ and $\sigma_{zz}=\sigma_{0}$, the above equation can be put in the form of the Poisson equation

[TABLE]

where $R=\sqrt{\sigma_{zz}/\sigma_{xx}}$. The boundary condition should also be determined by the charge continuity equation.

A large $R$ will greatly stretch the equipotential contour along $\hat{z}$ direction. It can cause a great portion of current lies in a sharp cone centered around the $\hat{z}$ axis Pippard (1989). This is similar to the behaviour of the longitudinal magnetoresistance due to chiral anomaly. In fact, if the injection of the current only localized near a small region in the surface and does not spread to the whole surface, the current injection may cause a strong negative magnetoresistance Huang et al. (2015); dos Reis et al. (2016); Yuan et al. (2016).

IV Outlook

We have shown how to derive various nonlinear currents in the framework of the semiclassical theory and Boltzmann transport theory. Band properties beyond the spectrum play essential roles in those nonlinear currents. However, the semiclassical theory ignores the electron interaction and the modification to electron-impurity scattering from electromagnetic fields. The semiclassical theory has been extended in Fermi liquid using Keldysh formalism up to first order Shindou and Balents (2006, 2008). Moreover, the Boltzmann transport theory is also extended to treat the drifting and collision in equal-footing up to first order using the Wigner distriubtion function and quantum kinetic theory Sekine et al. (2017); Culcer et al. (2017). It is tempting to extend such theories up to second order to give a complete treatment of various nonlinear currents.

For the three nonlinear anomalous Hall effects, the one due to Berry curvature dipole has attracted several experimental efforts, while no experiments are reported to measure the other two corrections. The correction due to magnetic field should be generally present in any system and competes with the ordinary Hall effect. Distinguishing the magneto nonlinear Hall effect from the ordinary Hall effect is similar to extracting the intrinsic anomalous Hall from total anomalous Hall signal Tian et al. (2009); Chun et al. (2007); Lee et al. (2004); Mathieu et al. (2004); Sales et al. (2006); Zeng et al. (2006).

Moreover, in samples that have combined time reversal and inversion symmetry (but break each one separately), the linear anomalous Hall effect as well as the nonlinear Hall effect due to Berry curvature dipole vanishes identically and only the intrinsic nonlinear Hall current contributes. This type of materials is ideal to test the intrinsic nonlinear Hall effect. If the time reversal and inversion symmetry, as well as their combined symmetry are all broken, the anomalous Hall and its two electric nonlinear version should all be present. There are still methods to distinguish them. One can first perform a scaling of the relaxation time to distinguish the intrinsic contribution from the extrinsic contribution. This is similar to extracting intrinsic, skew-scattering and side-jump contribution from the total anomalous Hall signal Tian et al. (2009); Chun et al. (2007); Lee et al. (2004); Mathieu et al. (2004); Sales et al. (2006); Zeng et al. (2006). The nonlinear part can be further extracted by reversing the direction of the electric field and finding the difference in the Hall signal.

For the negative longitudinal magnetoresistance, the difficulty is to distinguish the contribution from chiral anomaly from the intrinsic and current jetting contributions. The current jetting effect may be reduced by attaching the lead across the whole end of the sample so that the current is uniformly injected. In comparison, the main difference between the chiral anomaly and the intrinsic contribution is a hierarchy of relaxation time. One needs to control the ratio of the intervalley or chiral relaxation time to the transport relaxation time to distinguish these two contributions.

There is another issue in evaluating nonlinear currents in first-principles codes. Depending on the type of the unperturbed Hamiltonian, the second order energy has different forms for Schrodinger Hamiltonian or effective tight-binding Hamiltonian, mainly due to different Hessian matrices. Then there are two methods to evaluate the second order energies and related effects: one can start from the Schr$\ddot{\rm o}$dinger Hamiltonian and use the wave function from first-principles codes as approximation to evaluate the analytical expression derived for the Schr$\ddot{\rm o}$dinger Hamiltonian; on the other hand, one can start from the effective tight-binding Hamiltonian and uses its response to electromagnetic fields to approach the real one. The second one requires to evaluate the effective Hessian matrices while that of the first one is simply $1/m\delta_{ij}$. It is tempting to study which method is more accurate.

Acknowledgements.

We acknowledge useful discussions with Xiao Di. This work is supported by the Department of Energy, Basic Energy Sciences, Grant No. DE-SC0012509.