Dynamical axion state with hidden pseudospin Chern numbers in MnBi$_{2}$Te$_{4}$-based heterostructures

TL;DR

This paper predicts and characterizes a class of dynamical axion states in MnBi$_{2}$Te$_{4}$-based heterostructures, distinguished by hidden pseudospin Chern numbers, advancing the understanding of axion electrodynamics in topological materials.

Contribution

It introduces a method to identify dynamical axion states via hidden pseudospin Chern numbers and proposes their realization in MnBi$_{2}$Te$_{4}$ heterostructures.

Findings

Dynamical axion states can be distinguished by hidden pseudospin Chern numbers.

MnBi$_{2}$Te$_{4}$-based heterostructures can host these dynamical axion states.

The study facilitates experimental exploration of axion electrodynamics in condensed matter.

Abstract

Axion is a hypothetical elementary particle which was initially postulated to solve the charge conjugation-parity problem in particle physics. Interestingly, the axion state has emerged in effective theory of topological insulators and has attracted extensive attention in condensed matter physics. Time-reversal or inversion symmetry constrains the axion field $\theta$ to be quantized. When both the time-reversal and inversion symmetries are broken by, say, an antiferromagnetic order, the axion field $\theta$ could become unquantized and dynamical along with magnetic fluctuations, which is termed the dynamical axion field. Here, we reveal that a wide class of topological-insulator-based dynamical axion states could be distinguished from the normal-insulator-based ones by a hidden quantity derived from the pseudospin Chern number. Motivated by recent research on MnBi$_{2}$Te$_{4}$-family materials, we further show that such topological-insulator-based dynamical axion states can be hopefully achieved in MnBi$_{2}$Te$_{4}$-based heterostructures, which should greatly facilitate the study of axion electrodynamics in condensed matter physics.

Equations

$$\begin{split}\mathcal{H}_{0}(\mathbf{k})=&\epsilon_{0}(\mathbf{k})\mathbb{I}_{4\times 4}+\sum^{5}_{i=1}d_{i}(k)\Gamma^{i},\\ d_{1,2,...,5}(\mathbf{k})=&\Big{(}A_{2}k_{y},-A_{2}k_{x},A_{1}k_{z},M(\mathbf{k}),0\Big{)},\end{split}$$

$$C_{\pm }^S=\frac{1}{2\pi }\int d^{2}k{\Omega }_{\pm }^{\mu \nu },$$

$$C^{S}(k_{z})=\frac{1}{2}\Big{[}-\mathrm{sgn}(B_{2})+\mathrm{sgn}(M+B_{1}k_{z}^{2})\Big{]}.$$

$$\nu =(-1{)}^{C^S(k_z=0)+C^S(k_z=\infty )},$$

Full text

Introduction. Recently, many analogs of elementary particles, which were originally discovered in particle physics, have been found as quasi-particle excitations in condensed matter physics, such as, for example, the massless Weyl and Dirac fermions in topological semimetals Armitage et al. (2018) and mysterious Majorana fermions in topological superconductorsQi and Zhang (2011); Elliott and Franz (2015). Among them, the search for an axion state has attracted intensive interest in condensed matter physics Qi et al. (2008); Essin et al. (2009); Rosenberg and Franz (2010); Karch (2009); Bermudez et al. (2010); Nomura and Nagaosa (2011); Wan et al. (2012); Sekine and Nomura (2014); Wang et al. (2015); Morimoto et al. (2015); Turner et al. (2012); Wang and Zhang (2013); Goswami and Roy (2014); Lee et al. (2015); You et al. (2016); Wieder and Bernevig ; Yue et al. (2019); Xu et al. (2019); Hou and Wu ; Otrokov et al. (2019a); Zhang et al. (2019); Li et al. (2019); Chowdhury et al. (2019); Wu et al. (2016); Grauer et al. (2017); Mogi et al. (2017a); Okada et al. (2016); Mogi et al. (2017b); Dziom et al. (2017); Xiao et al. (2018); Varnava and Vanderbilt (2018); Gong et al. (2019); Chen et al. (2019a); Otrokov et al. (2019b); Gui et al. (2019); Allen et al. (2019); Chen et al. (2019b); Gooth et al. (2019). As a hypothetical elementary particle, the axion is a possible component of dark matter of the Universe Peccei and Quinn (1977). Axions gained new vigor from the effective theory of topological insulators (TIs) Qi et al. (2008); Qi and Zhang (2011), which is described by an axion Lagrangian density $\mathcal{L}_{\theta}=(\theta/2\pi)(\alpha/2\pi)\mathbf{E}\cdot\mathbf{B}$. Here, $\alpha=e^{2}/\hbar c$ is the fine structure constant, $\mathbf{E}$ and $\mathbf{B}$ represent electric and magnetic fields, and the dimensionless pseudoscalar $\theta$ is a parameter describing the insulator, which is referred to as the axion field in axion electrodynamics Wilczek (1987). Generically, $\theta$ is odd under time-reversal symmetry (TRS) and inversion symmetry (IS). With periodic boundary conditions in space-time, $\theta$ has a $2\pi$-periodicity and it is constrained to $\pi$ for time-reversal-invariant TIs and [math] for normal insulators (NIs). Such quantized $\theta$ in TIs will lead to the topological magnetoelectric effect Qi et al. (2008); Essin et al. (2009); Rosenberg and Franz (2010); Karch (2009); Bermudez et al. (2010); Vazifeh and Franz (2010); Nomura and Nagaosa (2011); Sekine and Nomura (2014); Wang et al. (2015); Morimoto et al. (2015) and unique magneto-optical Faraday and Kerr effects Tse and MacDonald (2010); Maciejko et al. (2010); Ochiai (2012); Mal’shukov et al. (2013); Wu et al. (2016).

Once both TRS and IS are broken, $\theta$ deviates from the quantized value of [math] or $\pi$, as has been shown for TIs in the presence of antiferromagnetic (AFM) order Li et al. (2010a); Wang et al. (2011) and antiferromagnet chromia Hehl et al. (2009). Intriguingly, an inherently fluctuating AFM order will induce fluctuations of the axion field $\theta$ as a dynamical variable with spatial and temporal dependence Li et al. (2010a); Wang et al. (2011); Sekine and Nomura (2014); Wang et al. (2016, ). Such a dynamical axion field could give rise to rich, interesting effects, such as, the dynamical chiral magnetic effect and anomalous Hall effect Wilczek (1987); Sekine and Nomura (2016); Zhang et al. , instability in an external electric field Ooguri and Oshikawa (2012); Imaeda et al. (2019), and axionic polaritons Li et al. (2010a).

Up to now, two important issues on the dynamical axion state still remain open. First, since the broken TRS invalidates the $Z_{2}$ topological invariant in TIs Fu and Kane (2007), it remains unknown whether it is available to distinguish the TI-based dynamical axion state  Li et al. (2010a) with a large static $\theta\sim\pi$ from the NI-based one with a small static $\theta\sim 0$ Coh et al. (2011). Second, the experimental demonstration of a TI-based dynamical axion state is significantly hindered by the lack of realistic materials. Here, we explicitly address the above two issues in the affirmative. First, based on an effective model analysis, we uncover a hidden quantity derived from pseudospin Chern numbers Sheng et al. (2006); Prodan (2009); Yang et al. (2011) to differentiate a wide class of TI- and NI-based dynamical axion states. Further, inspired by recent studies on the MnBi2Te4 family of intrinsic magnetic insulators Zhang et al. (2019); Li et al. (2019); Gong et al. (2019); Otrokov et al. (2019b); Vidal et al. (2019a); Chen et al. (2019a); Deng et al. ; Liu et al. ; Zhang et al. (2020); Sun et al. (2019); Vidal et al. (2019b), we find that TI-based dynamical axion states can be hopefully achieved in MnBi2Te4$/$Bi2Te3-based magnetic van der Waals (vdW) heterostructures through breaking both TRS and IS by specific structural designs. Interestingly, some natural vdW heterostructures of (MnBi2Te4)m(Bi2Te3)n have been experimentally fabricated Hu et al. (2020); Wu et al. (2019); Vidal et al. (2019b), which may provide an ideal platform to study the axion electrodynamics and to detect the axion dark matter Marsh et al. (2019).

Low-energy effective model. Without loss of generality, we start from the typical low-energy effective model of Bi2Te3-family TIs Zhang et al. (2009a) and then treat the AFM order as a perturbation Li et al. (2010a). In the nonmagnetic state with the preserved IS, the four low-lying states are the bonding states $|P1^{+}_{z},\uparrow(\downarrow)\rangle$ from Bi atoms and antibonding states $|P2^{-}_{z},\uparrow(\downarrow)\rangle$ from Te atoms, where the superscript “$+$” (“$-$”) represents positive (negative) parity. In the ordered basis of ($|P1^{+}_{z},\uparrow\rangle$, $|P1^{+}_{z},\downarrow\rangle$, $|P2^{-}_{z},\uparrow\rangle$, $|P2^{-}_{z},\downarrow\rangle$), the symmetry operations of the $D_{3d}$ group at $\Gamma$ are represented, for example, as TRS $\mathcal{T}=\tau_{0}\otimes i\sigma_{y}K$, IS $\mathcal{P}=\tau_{z}\otimes\sigma_{0}$, rotation symmetries $C_{3z}=\mathrm{exp}[\tau_{0}\otimes(-i\pi/3)\sigma_{z}]$ and $C_{2x}=\mathrm{exp}[\tau_{0}\otimes(-i\pi/2)\sigma_{x}]$, where $K$ denotes the complex conjugate, $\tau_{0}$ and $\sigma_{0}$ are $2\times 2$ identity matrices, and $\tau_{x,y,z}$ and $\sigma_{x,y,z}$ are Pauli matrices acting in the orbital and spin spaces, respectively. Keeping the most relevant terms up to quadratic order of $k$, the effective four-band Hamiltonian is given by

[TABLE]

where $\epsilon_{0}(\mathbf{k})=C+D_{1}k_{z}^{2}+D_{2}(k_{x}^{2}+k_{y}^{2})$, $M(\mathbf{k})=M+B_{1}k_{z}^{2}+B_{2}(k_{x}^{2}+k_{y}^{2})$, and the five Dirac matrices are represented as $\Gamma^{1,2,...,5}=(\tau_{x}\otimes\sigma_{x},\tau_{x}\otimes\sigma_{y},\tau_{y}\otimes\sigma_{0},\tau_{z}\otimes\sigma_{0},\tau_{x}\otimes\sigma_{z})$, which satisfy the Clifford algebra $\{\Gamma^{i},\Gamma^{j}\}=2\delta_{i,j}$. The TI phase requires the condition of $M/B_{i=1,2}<0$. Hereafter, we assume $B_{i=1,2}$ to be positive, and $M<0$ ($M>0$) corresponds to the TI (NI).

The AFM order breaks both $\mathcal{P}$ and $\mathcal{T}$ simultaneously, which is essential for the emergence of dynamical axion states, and to the leading order, it will give rise to the $\mathcal{P}$-, and $\mathcal{T}$- breaking mass term $\delta\mathcal{\mathcal{H}}_{\mathrm{AFM}}=m_{5}\Gamma_{5}$ Li et al. (2010a). Furthermore, if the AFM order preserves the combined $\mathcal{PT}$ symmetry, only the identity matrix $\mathbb{I}_{4\times 4}$ and the above five linearly independent anticommuting Dirac matrices $\Gamma^{1,2,...,5}$ are allowed Sup . Consequently, the total Hamiltonian is written as $\mathcal{H}(\mathbf{k})=\mathcal{H}_{0}(\mathbf{k})+m_{5}\Gamma_{5}$, and the antiunitary $\mathcal{PT}$ symmetry with $\mathcal{(PT)}^{2}=-1$ ensures the double degeneracy of each band. It should be emphasized that due to the nonvanishing AFM mass term, when tuning $M$ from negative (TI side) to positive (NI side), there will be no topological phase transition, but a crossover with an avoided gap closing. Nevertheless, the two magnetic states in the TI and NI sides, respectively, can still be differentiated by a hidden $Z_{2}$-like quantity derived from the pseudospin Chern numbers as we show below.

Pseudospin Chern number. The (pseudo)spin Chern number was originally proposed in the quantum spin Hall effect Sheng et al. (2006), which remains valid and robust even when TRS is broken or the (pseudo)spin is not conserved Sheng et al. (2006); Prodan (2009); Yang et al. (2011); Li et al. (2010b). The definition of the pseudospin Chern number relies on a smooth decomposition of the occupied valence bands into two pseudospin sectors through projecting a pseudospin operator $\hat{S}$ into valence bands as $h_{S}=P_{v}\hat{S}P_{v}$ with $P_{v}$ as the valence-band projecting operator Sheng et al. (2006); Prodan (2009); Yang et al. (2011); Sup . Diagonalization of $h_{S}$ then gives rise to two pseudospin branches with eigenvalues $\pm\epsilon_{S}(\mathbf{k})$, and corresponding eigenvectors $|\varphi_{\pm}\rangle$. For each pseudospin branch, the usual Chern number can be calculated in two-dimensional (2D) momentum space as

[TABLE]

where $\Omega_{\pm}^{\mu\nu}=-2\mathrm{Im}\langle\frac{\partial\varphi_{\pm}}{\partial k_{\mu}}|\frac{\partial\varphi_{\pm}}{\partial k_{\nu}}\rangle$ is the Berry curvature in the 2D $k_{\mu}$-$k_{\nu}$ slice. The pseudospin Chern number is then defined as $C^{S}=(C_{+}^{S}-C_{-}^{S})/2$.

By treating one momentum as a parameter, we can calculate $C^{S}$ for each 2D slice in the 3D Brillouin zone. Specifically, we treat $k_{z}$ as the parameter and investigate the 2D slices ($k_{x}$-$k_{y}$). We construct the pseudospin operator as $\hat{S}=\tau_{z}\sigma_{z}$ in the current basis of $\mathcal{H}(\mathbf{k})$ Sup . The $k_{z}$-dependent pseudospin Chern number can then be analytically derived as Sup

[TABLE]

Interestingly, it is independent of the AFM mass $m_{5}$. Since $B_{1}$ and $B_{2}$ are positive, $C^{S}(k_{z})$ always equals zero in the NI-based regime with $M>0$, whereas one can obtain $C^{S}(k_{z})=-1$ for $|k_{z}|<k_{zc}=\sqrt{-M/B_{1}}$ and $C^{S}(k_{z})=0$ for $|k_{z}|>k_{zc}$ in the TI-based regime with $M<0$, as shown in Fig. 1(a). It is noteworthy that the change of $C^{S}(k_{z})$ here is accompanied by the closing of the pseudospin-spectrum gap instead of the energy gap which is always kept open. This can be seen by the $k_{z}$-dependent pseudospin spectrum $\epsilon_{S}$ at (0,0,$k_{z}$) in Fig. 1(a) with a closing point at $|k_{z}|=k_{zc}$. Representative pseudospin spectra in the $k_{x}$-$k_{y}$ plane are also shown in Fig. 1(c) for $|k_{z}|<k_{zc}$ (left column), $|k_{z}|=k_{zc}$ (middle column), and $|k_{z}|>k_{zc}$ (right column), respectively. Moreover, we stress that the pseudospin Chern number remains quantized even though the pseudospin is not conserved for $k_{z}\neq 0$ ($\epsilon_{S}$ deviates from $\pm 1$) Prodan (2009). This can alternatively be reflected by the $U(1)$ Wilson loop spectra Yu et al. (2011) for each pseudospin branch in the 2D slice ($k_{x}$-$k_{y}$) of a lattice model Sup . As shown in Fig. 1(b), when $k_{x}$ goes from [math] to $2\pi$, the Wilson loop of the lower pseudospin branch exhibits one (no) winding for $C_{-}^{S}(k_{z})=+1(0)$.

In the TI-based magnetic state, such a change of $C^{S}$ from $k_{z}=0$ to $|k_{z}|>k_{zc}$ essentially results from the inverted band structure ($M<0$) at the $\Gamma$ point of the parent TI. This enables us to define a $Z_{2}$-like quantity $\nu$ as

[TABLE]

where $k_{z}=\infty$ should be replaced by $k_{z}=\pi$ in Bloch bands, and $\nu=-1$ (or $+1$) corresponds to the TI-based (NI-based) state. It is worth mentioning that when treating $k_{x}$ (or $k_{y}$) as a parameter, through similar procedures of calculating $k_{x}$- (or $k_{y}$-) dependent pseudospin Chern numbers Sup , identical $\nu$ can be obtained. Strictly speaking, $\nu$ is not a topological invariant since the above two magnetic states are topologically the same with an avoided gap-closing crossover. Nonetheless, a “nontrivial” $\nu=-1$ still acts as a useful quantity which captures the band-inversion information of the parent TI state. Such a quantity can be applied to a wide class of 3D magnetic insulators, as long as the time-reversal-invariant 2D plane could be effectively reduced to a Bernevig-Hughes-Zhang-like model Bernevig et al. (2006); Bernevig and Hughes (2013) by some unitary transformations without the TRS-breaking magnetic order.

Massive Dirac surface states. For a 3D TI, gapless surface states, with a single massless Dirac cone, emerge on the surface. However, due to the AFM mass term in TI-based states, surface states are gapped with a single massive Dirac cone on all surfaces, which can be clearly seen by transforming the above continuum model to a tight-binding one on a cubic lattice. The massive Dirac-cone surface states along the $z$ direction can be derived as h^{z}_{\mathrm{surf}}=\pm\big{[}A_{2}(k_{y}\sigma_{x}-k_{x}\sigma_{y})+m_{5}\sigma_{z}\big{]} Sup , where “$+$” and “$-$” refer to the bottom and top surfaces, respectively, indicating an energy gap of $2m_{5}$ for the surface bands, as shown in Fig. 1(d). Similarly, for the surface states along $x$ and $y$ directions, the effective Hamiltonians are given by $h^{x}_{\mathrm{surf}}=\mp(A_{2}k_{y}\sigma_{z}+A_{1}k_{z}\sigma_{y}+m_{5}\sigma_{x})$ and $h^{y}_{\mathrm{surf}}=\pm(-A_{2}k_{x}\sigma_{z}+A_{1}k_{z}\sigma_{x}-m_{5}\sigma_{y})$, respectively Sup . Note that all surfaces are gapped for our model and that no one-dimensional gapless hinge states are found, as shown in Fig. 1(e) by the energy spectrum with open-boundary conditions in both the $x$ and $y$ directions. This is different from the recently proposed (static) axion insulators with chiral hinge states between neighboring gapped surface states Yue et al. (2019); Xu et al. (2019), due to the lack of IS and higher-order topology in our system.

Magnetic van der Waals heterostructures. The recently discovered MnBi2Te4-family insulators host the static axion state with $\theta=\pi$, because the IS is preserved Zhang et al. (2019). It is promising to design (MnBi2Te4)m-(Bi2Te3)n superlattices, where both $\mathcal{T}$ and $\mathcal{P}$ are broken with suitably chosen $m$ and $n$, to realize the TI-based dynamical axion state. Here, we take (MnBi2Te4)2(Bi2Te3) as a concrete example, also written as Mn2Bi6Te11. Its crystal structure is schematically shown in Fig. 2. The basic stacking block consists of one Bi2Te3 quintuple layer (QL) sandwiched between two MnBi2Te4 septuple layers (SLs). Based on first-principles calculations, the magnetic structure is the $A$-type AFM order between nearest SLs and the out-of-plane ferromagnetic order in each SL. The exchange interaction in Mn2Bi6Te11 is expected to be weaker than that in MnBi2Te4 due to the intercalation layer Bi2Te3.

If ignoring the magnetism, Mn2Bi6Te11 takes the rhombohedral crystal structure with the same space group $D_{3d}^{5}$ ($R\bar{3}m$) as Bi2Te3. Each atomic layer has a triangular lattice structure, which is stacked in the ABC order Zhang et al. (2019). The $\mathcal{P}$ is preserved, with the inversion center at the middle Te site of each QL. Once the AFM order is considered, both $\mathcal{T}$ and $\mathcal{P}$ are broken, but their combination $\mathcal{PT}$ is preserved.

The total energies of different magnetic configurations are calculated, shown in Fig. S1 of the Suppplemental Material (SM) Sup . The (001)AFM state is found to be the magnetic ground state of Mn2Bi6Te11. The band structures of the (001)AFM state without and with spin-orbit coupling (SOC) are presented in Figs. 3(a) and 3(b), respectively. The fat band structures demonstrate that the highest valence bands and the lowest conduction bands are dominated by $p$ orbitals of Bi and Te atoms. The SOC induces the “band inversion” between the valence and conduction bands. The inverted band gap is about 0.18 eV, as seen from the density of states shown in Fig. 3(c). Furthermore, in Figs. 3(d) and 3(e), we plot the evolution of the energy gap with gradually tuning the SOC by $\lambda\ast\lambda_{0}$, where $\lambda_{0}$ denotes the realistic SOC strength. The energy gap first decreases to a minimum at about $\lambda=0.68$ [see the inset in Fig. 3(e)] where the “band inversion” occurs, and then increases again. The minimum gap $\sim 7.6$ meV results from the mass term induced by the AFM order.

The electronic structure of Mn2Bi6Te11 can be captured by the above low-energy effective model, and the detailed parameters are provided in the SM Sup , which can be obtained by fitting with band structures from first-principles calculations. The negative value of $M$ describes the inverted band structure. To investigate the surface states of Mn2Bi6Te11, we employed maximally localized Wannier functions Marzari and Vanderbilt (1997); Zhang et al. (2009b), shown in Fig. S3(a) Sup , where the MnBi2Te4 SL is chosen as the surface termination. The expected gapped Dirac surface state is clearly seen. We also calculate the Fermi surface in Fig. S3(b) Sup , which is consistent with the effective model analysis. When the energy level is increased, $k$-cubic terms must be taken into account, leading to a triangular Fermi surface, as shown in Fig. S3(c) Sup , in contrast to the hexagonal Fermi surface in Bi2Te3, since the $\mathcal{T}$ is broken.

Mn2Bi6Te11 has both inverted bulk bands and gapped surface states, and expects to realize the TI-based dynamical axion state with a hidden pseudospin Chern quantity $\nu=-1$. Moreover, based on the low-energy model, we have explicitly calculated the static part $\theta_{0}$ Li et al. (2010a); Sup of the axion field as a function of the SOC strength from $\lambda=0$ to $\lambda=1$, as shown by the blue line in Fig. 4(a). Due to the presence of the $m_{5}$ term, $\theta_{0}$ deviates from [math] and $\pi$. With increasing the SOC strength, $\theta_{0}$ continuously changes from near [math] with $M\gg|m_{5}|$ (NI-based phase) to a large value $\sim\pi$ with $M\ll-|m_{5}|$ (TI-based phase). A rapid increase of $\theta_{0}$ is found around the (crossover) transition point near $\lambda=0.68$ with a minimum gap, which is consistent with the first-principles calculations in Fig. 3(e). Intriguingly, in the presence of AFM fluctuations, the axion field will become dynamical with $\theta=\theta_{0}+\delta\theta(\mathbf{x},t)$, and to linear order, $\delta\theta(\mathbf{x},t)=\delta m_{5}/g$, where $\delta m_{5}$ is proportional to the amplitude fluctuation of AFM order along the $z$ direction, and $g$ is the coefficient Li et al. (2010a).

Discussion. The axion electrodynamics may lead to various exotic physical effects Wilczek (1987); Sekine and Nomura (2016); Zhang et al. ; Ooguri and Oshikawa (2012); Imaeda et al. (2019); Li et al. (2010a), which in turn makes it possible to experimentally detect the dynamical axion field. For example, the axionic polariton, as a collective mode emerging from the coupling between the light and axionic mode, can be detected by the attenuated total reflection measurement Li et al. (2010a). To facilitate the experimental detection of the dynamical axion state, a large proportionality coefficient $g^{-1}$ between $\delta\theta$ and the AFM fluctuation $\delta m_{5}$ is highly preferred. In Fig. 4(a), we plot $g^{-1}$ as a function of the SOC strength for Mn2Bi6Te11 (see the red line). It clearly shows that the value of $g^{-1}$ in the TI-based side ($M<0$) is generically one order of magnitude larger than that in the NI-based side ($M>0$). Notably, $g^{-1}$ becomes significantly large when approaching the transition point from the TI-based side ($M\rightarrow 0^{-}$)Zhang et al. . This can be achieved via a chemical replacement by lighter elements with weaker SOC (e.g., Bi by Sb, Te by Se). We present the detailed band structures (see Fig. S4 Sup ) and gap evolution with tuning the SOC strength in Fig. 4(b) (insets around minimum energy gaps are shown in Fig. S5 Sup ) for the family superlattices $X_{2}Y_{6}Z_{11}$ ($X$=Eu/Mn, $Y$=Bi/Sb, and $Z$=Se/Te), some of which are also promising candidates to realize the TI-based dynamical axion state.

Acknowledgement. H.Z. was supported by the National Natural Science Foundation of China (Grants No. 11674165 and 11834006), the Fok Ying-Tong Education Foundation of China (Grant No. 161006), the Fundamental Research Funds for the Central Universities (Grant No. 021314380147) and High-Performance Computing Center of Collaborative Innovation Center of Advanced Microstructures. J.W. is supported by the Natural Science Foundation of China through Grant No. 11774065, the National Key Research Program of China under Grant No. 2016YFA0300703, and the Natural Science Foundation of Shanghai under Grant Nos. 17ZR1442500, and 19ZR1471400.

H.W. and D.W. contributed equally to this work.