Nonlinear response in a non-centrosymmetric topological insulator

TL;DR

This paper investigates nonlinear optical responses in non-centrosymmetric topological insulators, explaining frequency up-conversion phenomena through a generalized Kubo response theory.

Contribution

It introduces a theoretical framework based on a generalized Kubo formula to analyze nonlinear responses in topological insulators with broken inversion symmetry.

Findings

Nonlinear response theory explains frequency up-conversion.

Breaking inversion symmetry enables second-harmonic generation.

Hexagonal warping affects nonlinear optical properties.

Abstract

Nonlinear phenomena are inherent in most systems in nature. Second or higher-order harmonic generations, three-wave and four-wave mixing are typical phenomena in nonlinear optics. To obtain a nonzero signal for second-harmonic generation in the long-wavelength limit ($q\rightarrow0$), the breaking of inversion symmetry is required. In topological materials, a hexagonal warping term which breaks the rotation symmetry of the Fermi surface is observed by angular-resolved photo-emission spectroscopy (ARPES). If a gap opens (e.g., by doping with magnetic impurities) the inversion symmetry will be broken. Here we use a nonlinear response theory based on a generalized Kubo formula to explain the frequency up-conversion in topological materials.

Equations

$${\tilde{\sigma }}_{xxx}(\omega ,\omega )=\frac{e^3}{{\omega }^2}\frac{i}{4{\pi }^2}{\int }_0^{2\pi }d\theta {\int }_0^{k_{\mathrm{cut}}}kdkT\sum _l\mathrm{Tr}⟨v_x\hat{G}(\mathbf{k},i{\omega }_l)v_x\hat{G}(\mathbf{k},i{\omega }_l+i{\omega }_n)v_x\hat{G}(\mathbf{k},i{\omega }_l-i{\omega }_n){⟩}_{i{\omega }_n\to \omega +i\delta },$$

$$\displaystyle\tilde{\sigma}^{\textrm{inter}}_{xxx}=\frac{ie^{3}}{\hbar^{3}\omega^{2}\pi^{2}}\int kdkd\theta F(k,\theta)\frac{[f(E)-f(-E)]}{E}\Big{[}\frac{1}{\hbar\omega+i\delta+2E}-\frac{1}{\hbar\omega+i\delta-2E}+\frac{2}{\hbar\omega+i\delta-E}-\frac{2}{\hbar\omega+i\delta+E}\Big{]}$$

$$H_0=v_k(k_x{\sigma }_y-k_y{\sigma }_x)+\frac{\lambda }{2}(k_{+}^3+k_{-}^3){\sigma }_z+M{\sigma }_z,$$

$$F(k,\theta )$$

$${\tilde{\sigma }}_{xx}(\omega )=\frac{e^2}{i\omega }\frac{1}{4{\pi }^2}{\int }_0^{k_{cut}}kdkd\theta \times$$

$$T\sum _l\mathrm{Tr}⟨v_x\hat{G}(\mathbf{k},i{\omega }_l)v_x\hat{G}(\mathbf{k},i{\omega }_n+i{\omega }_l){⟩}_{i{\omega }_n\to \omega +i\delta }$$

$${\tilde{\sigma }}_{xx}^{\text{inter}}(\omega )=\frac{i{e}^2}{4{\pi }^2{\hslash }^2\omega }{\int }_0^{k_{cut}}kdkd\theta H(k,\theta )$$

$$\displaystyle\frac{f(E)-f(-E)}{E}\Big{[}\frac{1}{\hbar\omega+i\delta+2E}-\frac{1}{\hbar\omega+i\delta-2E}\Big{]}$$

$$H(k,\theta )=\frac{v_k^2}{E}[9{\lambda }^2{k}^6{\cos }^2(2\theta )+(M+\lambda k^3\cos (3\theta ){)}^2+$$

$$v_k^2{k}^2{\sin }^2\theta -6\lambda k^3\cos (2\theta )\cos \theta (M+\lambda k^3\cos (3\theta ))]$$

$$\hat{G}(\mathbf{k},i{\omega }_n)=\frac{1}{2}\sum _{s=\pm }(1+s{\mathbf{F}}_k\cdot \sigma )G_0(\mathbf{k},s,i{\omega }_n),$$

$${\mathbf{F}}_k=\frac{(-v_{k}k\sin \theta ,v_{k}k\cos \theta ,\lambda k^3\cos (3\theta )+M)}{\sqrt{v_k^2{k}^2+[\lambda k^3\cos (3\theta )+M{]}^2}},$$

$$E=\sqrt{v_k^2{k}^2+[\lambda k^3\cos (3\theta )+M{]}^2}.$$

$$v_x$$

$$v_y$$

$$\mathbf{A}\cdot \mathbf{B}=a_0{b}_0+\mathbf{a}\cdot \mathbf{b},$$

$$\mathbf{A}\times \mathbf{B}=a_0\mathbf{b}+b_0\mathbf{a}+i(\mathbf{a}\times \mathbf{b}).$$

$$(a_0+\mathbf{a}\cdot \sigma )(g_{01}+{\mathbf{g}}_1\cdot \sigma )=a_0{g}_{01}+g_{01}\mathbf{a}\cdot \sigma +a_0{\mathbf{g}}_1\cdot \sigma +\mathbf{a}\cdot {\mathbf{g}}_1+i(\mathbf{a}\times {\mathbf{g}}_1)\cdot \sigma$$

$$=\mathbf{A}\cdot {\mathbf{G}}_1+(\mathbf{A}\times {\mathbf{G}}_1)\cdot \sigma .$$

$$\mathrm{Tr}⟨(a_0+\mathbf{a}\cdot \sigma )(g_{01}+{\mathbf{g}}_1\cdot \sigma )(a_0+\mathbf{a}\cdot \mathbf{p})(g_{02}+{\mathbf{g}}_2\cdot \sigma )(a_0+\mathbf{a}\cdot \sigma )(g_{03}+{\mathbf{g}}_3\cdot \sigma )⟩$$

$$=$$

$$+i(A\times G_1)\times (A\times G_2)\cdot (A\times G_3).$$

$$\hat{G}(\mathbf{k},i{\omega }_n)={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\frac{\hat{A}(\mathbf{k},\omega )}{i{\omega }_n-\omega },$$

$${\tilde{\sigma }}_{xxx}(\omega ,\omega )$$

$$\hat{A}(\mathbf{k},\omega )=A_I(\mathbf{k},\omega )+A_x(\mathbf{k},\omega ){\sigma }_x+A_y(\mathbf{k},\omega ){\sigma }_y+A_z(\mathbf{k},\omega ){\sigma }_z.$$

$$A_I(\mathbf{k},\omega )=\delta (\omega +\mu -E)+\delta (\omega +\mu +E),$$

$$A_x(\mathbf{k},\omega )=\frac{-v_{k}k\sin \theta [\delta (\omega +\mu -E)-\delta (\omega +\mu +E)]}{\sqrt{v_k^2{k}^2+[\lambda k^3\cos (3\theta )+M{]}^2}},$$

$$A_y(\mathbf{k},\omega )=\frac{v_{k}k\cos \theta [\delta (\omega +\mu -E)-\delta (\omega +\mu +E)]}{\sqrt{v_k^2{k}^2+[\lambda k^3\cos (3\theta )+M{]}^2}},$$

$$A_z(\mathbf{k},\omega )=\frac{[\lambda k^3\cos (3\theta )+M][\delta (\omega +\mu -E)-\delta (\omega +\mu +E)]}{\sqrt{v_k^2{k}^2+[\lambda k^3\cos (3\theta )+M{]}^2}}.$$

$$\mathrm{Tr}⟨v_x\hat{A}(\mathbf{k},{\omega }_1)v_x\hat{A}(\mathbf{k},{\omega }_2)v_x\hat{A}(\mathbf{k},{\omega }_3)⟩$$

$$=$$

$$+\delta ({\omega }_1+\mu -E)\delta ({\omega }_2+\mu -E)\delta ({\omega }_3+\mu +E)-\delta ({\omega }_1+\mu -E)\delta ({\omega }_2+\mu +E)\delta ({\omega }_3+\mu +E)$$

$$-\delta ({\omega }_1+\mu +E)\delta ({\omega }_2+\mu -E)\delta ({\omega }_3+\mu +E)-\delta ({\omega }_1+\mu +E)\delta ({\omega }_2+\mu +E)\delta ({\omega }_3+\mu -E)]$$

$$+$$

$$F(k,\theta )$$

Full text

I introduction

The nonlinear response to an external driving electromagnetic field $\mathbf{E}$ or $\mathbf{B}$ can be characterized by a conductivity tensor $\tilde{\sigma}$ which is not a constant, but depends on the magnitude of $\mathbf{E}$ or $\mathbf{B}$. Nonlinearity is often found to be important in optical devices, especially in the recent discovery of high-efficiency solar energy harvesting in non-centrosymmetric crystal structures such as perovskite oxides Grinberg ; Nie ; Shi ; Quile ; Cook . In 3D topological insulator (TI) and ferroelectric materials, Dirac cones Hasan ; Qi1 ; Moore ; Hsieh1 ; Chen obeying spin-momentum locking Hsieh2 ; Nori with in-plane spin component perpendicular to the momentum $\mathbf{k}$ were verified by spin-sensitive angular-resolved photo emission spectroscopy (ARPES). The quasiparticles (helical Dirac fermions) observed in topological materials possess an important feature: the Fermi contours are circular for small values of the chemical potential $\mu$, and acquire a snowflake Chen shape as $\mu$ increases. Analyzing the experiment results, Fu Fu assigned a hexagonal warping term in the Hamiltonian of such quasiparticles. This term has a strong signature in the optical conductivity, spin texture, Hall conductivity and circular dichroism of topological insulators Li1 ; Xiao3 ; Li2 .

The optical conductivity was predicted to show Li1 a large near-linear increase with photon energy above the inter-band threshold as compared to the usual flat background Carbotte1 ; Li08 ; Orlita inter-band optical conductivity in graphene. The spin texture Xu1 (specifically, out of plane spin $S_{z}$) shows a mixture of up-and-down directions; in contrast to the normal all-up or all-down hedgehog type Xu2 distribution for massive Dirac fermions (see, e.g. Fig. 5 of Ref. [Li2, ]). It is also possible to introduce a gap in the topological surface quasiparticles (massive Dirac fermions) by magnetic doping Chen1 ; Tokura in B$\textrm{i}_{2}$S$\textrm{e}_{3}$ Chen1 and recently in $\textrm{Cr}_{x}(\textrm{Bi}_{1-y}\textrm{Sb}_{y})_{2\text{\textminus}x}\textrm{Te}_{3}$ Tokura ; Yasuda . Considerable particle-hole asymmetry of the surface Dirac cone of a 3D TI usually displays, which can be modeled with a small sub-dominant Schrödinger quadratic-in-momentum term in addition to the dominant Dirac Hamiltonian. While perhaps small, the Schrödinger term has been shown to provide important modifications Li3 in the chiral nonlinear magneto optical conductivity (MOC) which is related to the absorption of left and right circularly polarised light of 3D TI. This is to be compared with what is found in graphene Carbotte4 ; Xiao1 ; Yao or the related single layer silicene Tabert .

In this work we focus on the nonlinear optical conductivity induced by an electric field $\mathbf{E}$ in contrast to the nonlinear MOC which is induced by a magnetic field $\mathbf{B}$. We consider three-wave mixing Shen (e.g., second-harmonic generation) from non-centrosymmetric topological materials. Second-harmonic generation (SHG) was first demonstrated by projecting a laser beam through crystalline quartz Franken . Later on this effect was found in other materials (e.g., silicon surfaces) Heinz with broken inversion symmetry. Theoretically, SHG was predicted to be nonzero in semiconductors Sipe , and more recently in single-layer graphene Mikhailov ; Daria with oblique incidence of radiation on the 2D electron layer. For oblique incidence, the incident radiation has a nonzero wave vector component $\mathbf{q}$ parallel to the plane of the 2D layer. In the long-wavelength limit ($\mathbf{q}\rightarrow 0$, normal incidence), the SHG vanishes because graphene is a centrosymmetric material. However, higher-order harmonics (e.g., third-harmonic generation) could be nonzero in graphene Wright or generally Dirac Fermions system Nagaosa . The nonlinear coupling of three monochromatic waves, thus called three-wave mixing, has been successfully used to generate optical frequency up-conversion or down-conversion. Nonlinear optical analogs, including SHG, have also been studied recently in various contexts, including Josephson plasma waves Savel and cavity quantum electrodynamics Kockum ; Kockum1 ; Stassi ; Gu .

In the following paragraphs we present a Green’s function formalism for calculating the nonlinear conductivity in section II. We use a two-band hexagonal warping model which can be found in surface states of 3D TI and ferroelectric materials. The inversion symmetry of the Fermi surface is broken by a magnetic doping in the hexagonal warping model. In section III we present the linear optical conductivity from the Green’s function formalism. In section IV we present our numerical results of the nonlinear and linear conductivity for different sets of parameters (e.g., chemical potential, gap parameter, temperature, etc.). In section V we summarize our results with a conclusion.

II Nonlinear optical conductivity

The linear conductivity $\tilde{\sigma}_{xx}(\omega)$ is related to the current $J_{x}(\omega)=\tilde{\sigma}_{xx}(\omega)E_{x}(\omega)$, while the nonlinear conductivity $\tilde{\sigma}_{xxx}(\omega,\omega)$ is related to the current $J_{x}(2\omega)=\tilde{\sigma}_{xxx}(\omega,\omega)E^{2}_{x}(\omega)$. In general, the nonlinear conductivity is a tensor $\tilde{\sigma}_{\alpha\beta\gamma}$. However, here for simplicity we only consider the $xxx$ component of the tensor; the other components of the conductivity tensor can be obtained in a similar way. The nonlinear conductivity has been well studied in earlier references; for example, in Ref. [Bloemb, ] the Eq. (2-48) defines the nonlinear conductivity as a product of momentum matrix elements and then in Eq. (2-49) the momentum matrix elements were connected to velocity matrix elements.

For the linear conductivity it has been shown in chapter 8 of the book [Economou, ] that the Eq. (8.53) uses a trace of momentum operators and Green’s functions and then in Eq. (8.55) this was connected to the product of velocity matrix elements. The velocity matrix element is connected to the position matrix element and the shift vector. ZhouLi For the nonlinear conductivity, instead of using velocity matrix elements Bloemb directly, we define the nonlinear conductivity as a trace of velocity operators and Green’s functions Economou ; the imaginary frequency in each Green’s function is set by using a triangle Feynman diagram.

[TABLE]

here $v_{x}$ is the velocity operator and $\widehat{G}(\mathbf{k,}i\omega_{l})$ is the matrix Green’s function, $e$ is the charge of the electron, $k$ the absolute value of the momentum $\mathbf{k}$ with direction $\theta$ and cutoff $k_{\mathrm{cut}}$, $T$ is the temperature with $\omega_{n}=2n\pi T$, $\omega_{l}=(2l+1)\pi T$ the Boson and Fermion Matsubara frequencies, $n$ and $l$ are integers and $\mathrm{Tr}$ is a trace. To obtain the nonlinear conductivity, which is a real frequency quantity, we needed to make an analytic continuation from imaginary $i\omega_{n}$ to real $\omega$ and $\delta$ is infinitesimal. This is valid for the long wavelength limit $q\rightarrow 0$;

Consider a two-band model as an example, the velocity operators and matrix Green’s functions are 2 $\times$ 2 matrices, and can be expanded onto the basis of Pauli matrices $\mathbf{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$, as $v_{x}=a_{0}+\mathbf{a}\cdot\mathbf{\sigma}$, and $\widehat{G}(\mathbf{k},i\omega_{n})=g_{0}+\mathbf{g}\cdot\mathbf{\sigma}$. We can use the algebra $(\mathbf{a}\cdot\mathbf{\sigma})(\mathbf{g}\cdot\mathbf{\sigma})=(\mathbf{a}\cdot\mathbf{g)}I_{2}+i(\mathbf{a}\times\mathbf{g)\cdot\mathbf{\sigma}}$ to evaluate the trace and the complicated results will be contained in the function $F(k,\theta)$ to be integrated further in momentum space $(k,\theta)$. We can also perform the sum over the internal Fermion Matsubara frequencies $\omega_{l}$ and the result is Fermi-Dirac distribution function defined as $f(x)=1/[\exp(x/T-\mu/T)+1]$. After tedious but straightforward algebra (details in the appendix), we finally obtained both the inter-band and intra-band nonlinear optical conductivity, the intra-band optical conductivity contributes to the frequency region of $\omega\approx 0$ and was given in the appendix, in the equations below we present the results of the inter-band optical conductivity ($\tilde{\sigma}^{\textrm{inter}}_{xxx}=\tilde{\sigma}^{\textrm{inter}}_{xxx}(\omega,\omega)$),

[TABLE]

Here $E$ is the quasiparticle energy which depends on the momentum $(k_{x},k_{y})=(k\mathrm{cos}(\theta),k\mathrm{sin}(\theta))$. Take a two-band hexagonal warping model as an example, the Hamiltonian is given by,

[TABLE]

this model has been used to describe the surface states band structure near the $\Gamma$ point in the surface Brillouin zone of a 3D TI and also recently in ferroelectric materials. The Dirac fermion velocity to second order is $v_{k}=\hbar v_{\mathrm{F}}(1+\alpha k^{2})$, with $v_{F}$ the usual Fermi velocity and $\hbar v_{\mathrm{F}}$ measured to be 2.55 eV$\cdot\mathrm{\mathring{A}}$ and $\alpha$ is a constant which is fit along with $m$ to the measured band structure in Ref. [Fu, ]. Here $m$ appears in the quadratic term $\hbar^{2}k^{2}/(2m)$ which, for simplicity, is dropped in the Hamiltonian $H_{0}$. The inclusion of the quadratic term provides particle-hole asymmetry; however the wave function is not changed Li2 , thus the Berry curvature and Berry connection (defined from the wave function) are not modified by this quadratic term. For simplicity, the quadratic correction to the velocity $\alpha$ is also discarded. The magnitude of the hexagonal warping parameter is $\lambda=$200 eV$\cdot\mathrm{\mathring{A}}^{3}$, estimated from the measured Fermi velocity. The same value was used in Ref. [Fu, ]. The $\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$ are Pauli matrices here referring to spin, while in graphene these would relate instead to pseudospin. Finally $k_{\pm}=k_{x}\pm ik_{y}$, with the $k_{x}$, $k_{y}$ momentum along the $x$ and $y$ axis, respectively. $M$ is the strength of the gap opened when the topological thin film is in proximity to magnetic impurities.

The quasiparticle energy dispersion relation is given by $E=\sqrt{v_{k}^{2}k^{2}+[\lambda k^{3}\cos(3\theta)+M]^{2}}$, and the function $F(k,\theta)$ is given by

[TABLE]

Note that if $M=0$ or $\lambda=0$ the integration $\int_{0}^{2\pi}d\theta\int_{0}^{k_{\mathrm{cut}}}kdk$ will be zero because the integrand is an odd function of $k_{x}$. So only when both $M\neq 0$ and $\lambda\neq 0$ we obtain a non-vanishing second-harmonic generation nonlinear conductivity in the long-wavelength limit $q\rightarrow 0$.

III Linear optical conductivity

It is well known Li2 that the linear optical conductivity is obtained from the standard Kubo formula in terms of the matrix Green’s function and velocity operators, the longitudinal conductivity is given by

[TABLE]

which works out to be

[TABLE]

where the function $H(k,\theta)$ is given by

[TABLE]

It is interesting to check the units of $H(k,\theta)$ and $F(k,\theta)$. We find that $F(k,\theta)\times k$ and $H(k,\theta)$ have the same unit as $v_{k}^{2}E$. So $\tilde{\sigma}_{xxx}$ has the same unit as $\tilde{\sigma}_{xx}\times(e/\hbar\omega k)$. Then the product of the nonlinear conductivity $\tilde{\sigma}_{xxx}$ and external electric field, $\tilde{\sigma}_{xxx}\times E_{x}$ has the same unit as $\tilde{\sigma}_{xx}$, as expected.

IV Numerical Results

To evaluate the nonlinear optical conductivity, we need to perform an integration in momentum space which is restricted by the Fermi-Dirac distribution function $f(x)$. At zero temperature, the restricted area is the Fermi surface shown in Fig. 1. In (a), for a small chemical potential $\mu=0.1$ eV, the Fermi surface is very close to but not a perfect circle because of the small gap $2M=20$ meV. At larger chemical potential $\mu=0.4$ eV, the Fermi surface deviates a bit from a snowflake shape. In (b) a much larger gap $2M=200$ meV is used and the Fermi surface is significantly distorted. The inversion symmetry is broken in both (a) and (b).

In Fig. 2 and Fig. 3 we plot the numerical results of the real part of the inter-band nonlinear optical conductivity $\tilde{\sigma}^{\textrm{inter}}_{xxx}(\omega,\omega)$ and linear optical conductivity $\tilde{\sigma}^{\textrm{inter}}_{xx}(\omega)$, respectively. In Fig. 2, we find that if the chemical potential $\mu$ is smaller than half the gap $M$, the onset frequency is $M$. Because the chemical potential $\mu$ lies in the gap, the minimum energy for the inter-band transition is $2M$. The energy of absorbing two photons is $2\omega$, so that the onset frequency $2\omega\geq 2M$. If the chemical potential $\mu$ is larger than $M$, the onset frequency is the chemical potential $\mu$, because in this case $2\omega\geq 2\mu$. There is a small drop in the nonlinear optical conductivity at $\omega=2\mu$, because this is the onset frequency for another inter-band transition involving one photon absorbing. Thus the number of photons in the process of frequency doubling decreases. For $\omega\geq 2\mu$, curves with different values of $\mu$ fall on top of each other.

In Fig. 3 we find that the onset frequency of the linear optical conductivity is $2\mu$, in contrast to the onset frequency $\mu$ of the nonlinear optical conductivity. When the frequency $\omega$ is larger than the onset frequency $2\mu$, the linear optical conductivity warps up, in contrast to the nonlinear optical conductivity which decreases as the frequency increases. Curves with different values of $\mu$ also fall on top of each other for the linear optical conductivity. In Fig. 4 and Fig. 5 we show the corresponding imaginary parts of the inter-band optical conductivity $\tilde{\sigma}_{xxx}^{\textrm{inter}}(\omega,\omega)$ and $\tilde{\sigma}_{xx}^{\textrm{inter}}(\omega)$: respectively. The absolute value of the imaginary part of the nonlinear conductivity decreases to zero faster than that of the linear conductivity.

V Conclusion

In conclusion, we developed a method based on the trace of the velocity operator and Green’s function to calculate the nonlinear response functions in a non-centrosymmetric topological insulator. Our method is equivalent to the velocity matrix element method Bloemb if a two-band free-electron approximation was considered. We obtained the nonlinear conductivity for frequency up-conversion in the second harmonic generation. In the model used here (the two-band hexagonal warping model), the energy scale is around 200 meV in the far-infrared region, relevant for the thermal energy. This model describes surface states of 3D TI. If the 3D TI $(\textrm{Bi}_{1-y}\textrm{Sb}_{y})_{2}\textrm{Te}_{3}$ (BST) is doped with magnetic impurities, a small gap is opened in $\textrm{Cr}_{x}(\textrm{Bi}_{1-y}\textrm{Sb}_{y})_{2\text{\textminus}x}\textrm{Te}_{3}$ (CBST) Tokura ; Yasuda thus the inversion symmetry was broken.

Acknowledgements.

The authors thank Chong Wang, Yong Xu, A. W. Frisk Kockum, Mauro Cirio and Zubin Jacob for useful discussions. Z. L. acknowledges the support of a JSPS Foreign Postdoctoral Fellowship under Grant No. PE14052 and P16027. F.N. is supported in part by the: MURI Center for Dynamic Magneto-Optics via the Air Force Office of Scientific Research (AFOSR) (FA9550-14-1-0040), Army Research Office (ARO) (Grant No. W911NF-18-1-0358), Asian Office of Aerospace Research and Development (AOARD) (Grant No. FA2386-18-1-4045), Japan Science and Technology Agency (JST) (the Q-LEAP program, the ImPACT program and CREST Grant No. JPMJCR1676), Japan Society for the Promotion of Science (JSPS) (JSPS-RFBR Grant No. 17-52-50023, JSPS-FWO Grant No. VS. 059.18N), RIKEN-AIST Challenge Research Fund, and the John Templeton Foundation.

Appendix A Derivation of the nonlinear optical conductivity from the Green’s function technique

In this appendix we present a general formula for the calculation of the nonlinear conductivity tensor. We also expand the imaginary frequency Green’s function into a sum of the real frequency spectral function which can be measured directly from ARPES experiments. From this expansion we derive the intra-band and inter-band contribution to the nonlinear conductivity. We present how to perform the sum in imaginary frequency and obtain concise results in the non-interacting electron approximation. Consider the hexagonal warping model as an example. The Green’s function can be rewritten in the basis of Pauli matrices

[TABLE]

where

[TABLE]

and $G_{0}(\mathbf{k},s,i\omega_{n})=[{i\omega_{n}+\mu-sE}]^{-1}$, where $s=\pm 1$ and the energy spectrum is given by

[TABLE]

The velocity operator can be obtained as (for simplicity we set $\hbar=1$ )

[TABLE]

In general $v_{x}=a_{0}+\mathbf{a}\cdot\mathbf{\sigma}$ , $v_{y}=b_{0}+\mathbf{b}\cdot\mathbf{\sigma}$, and $\widehat{G}(\mathbf{k},i\omega_{n})=g_{0}+\mathbf{g}\cdot\mathbf{\sigma}$. If we define $\mathbf{A}=(a_{0},\mathbf{a}),\mathbf{B}=(b_{0},\mathbf{b}),\mathbf{G}=(g_{0},\mathbf{g})$, and use the following rules for dot and cross product of two vectors

[TABLE]

[TABLE]

Then the products of $v_{x}\widehat{G}(\mathbf{k,}\omega_{1})$ can be evaluated as

[TABLE]

The trace can be carried out in general as

[TABLE]

The matrix Green’s function $\widehat{G}(\mathbf{k},i\omega_{n})$ can be conveniently written in terms of a matrix spectral function $\widehat{A}(\mathbf{k},\omega)$ with

[TABLE]

then the conductivity in the long-wavelength limit becomes

[TABLE]

For two-band models, the spectral function $\widehat{A}(\mathbf{k},\omega)$ can be expanded in the basis of Pauli matrices,

[TABLE]

In the free-electron approximation (ignoring impurity scattering and electron-phonon scattering), the spectral functions are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The trace can be carried out as

[TABLE]

where we have defined two functions,

[TABLE]

[TABLE]

These terms can be separated into inter-band and intra-band contributions to the nonlinear conductivity.

Appendix B Intra-band nonlinear conductivity

The intra-band nonlinear conductivity includes those terms proportional to $\delta(\omega_{1}+\mu-E)\delta(\omega_{2}+\mu-E)\delta(\omega_{3}+\mu-E)$ and $\delta(\omega_{1}+\mu+E)\delta(\omega_{2}+\mu+E)\delta(\omega_{3}+\mu+E)$, which will contribute to the zero-frequency DC conductivity. Performing the sum over Matsubara frequencies, we obtain

[TABLE]

And the intra-band conductivity becomes

[TABLE]

the intra-band conductivity can be numerically evaluated, by replacing the $\delta$ function with the broadened Lorentzian function. One can also evaluate the intra-band conductivity analytically; one example was given in the appendix of Ref. [Li1, ].

Appendix C Inter-band nonlinear conductivity

The other terms like $\delta(\omega_{1}+\mu-E)\delta(\omega_{2}+\mu+E)\delta(\omega_{3}+\mu-E)$ are included in the inter-band nonlinear conductivity which will contribute to the nonzero-frequency AC conductivity, written as

[TABLE]

which is further simplified as

[TABLE]

Finally we obtained

[TABLE]