Electric Quantum Oscillation in Weyl Semimetals

Electronic transport in Weyl semimetals is quite extraordinary due to the topological property of the chiral anomaly generating the charge pumping between two distant Weyl nodes with opposite chiralities under parallel electric and magnetic fields. Here, we develop a full nonequilibrium quantum transport theory of the chiral anomaly, based on the fact that the chiral charge pumping is essentially nothing but the Bloch oscillation. Specifically, by using the Keldysh nonequilibrium Green function method, it is shown that there is a rich structure in the chiral anomaly transport, including the negative magnetoresistance, the non-Ohmic behavior, the Esaki-Tsu peak, and finally the resonant oscillation of the DC electric current as a function of electric field, called the electric quantum oscillation. We argue that, going beyond the usual behavior of linear response, the non-Ohmic behavior observed in BiSb alloys can be regarded as a precursor to the occurrence of electric quantum oscillation, which is both topologically and energetically protected in Weyl semimetals.

Among the greatest mysteries in physics is the asymmetry between matter and antimatter in our known universe. Baryogenesis is the hypothesized physical process producing such an asymmetry, whose precise mechanism still remains elusive. The chiral anomaly, also known as the Adler-Bell-Jackiw anomaly [1][2][3], is generally regarded as one of the most crucial elements in the mechanism of baryogenesis.
Here, we propose a full quantum signature of the chiral anomaly, which is fundamentally due to the quantization of the chiral charge pumping under strong electric fields. A main breakthrough in this work is the realization that the chiral charge pumping is essentially nothing but the Bloch oscillation in the zeroth, or chiral Landau level (LL), which is quantized to generate robust Wannier-Stark ladder (WSL) eigenstates [21][22][23] topologically protected by the chiral anomaly. Albeit somewhat less, robust WSL eigenstates can be also formed in nonchiral LLs due to the energetic protection of the Bloch oscillation in Weyl semimetals.
The formation of WSL eigenstates reveals an intriguing similarity between electricity and magnetism. The quantized cyclotron motion of electrons under strong magnetic fields gives rise to well-known magnetic quantum oscillations [24] such as de Haas-van Alphen, Shubnikov-de Haas, and eventually the quantum Hall effects. Similarly, the quantized Bloch oscillation of electrons under strong electric fields can give rise to an electric-field-induced oscillation of the DC electric current, which we call the electric quantum oscillation (EQO).
Actually, the EQO brings out one of the most fundamental differences between electricity and magnetism. That is, electric fields inevitably cause nonequilibrium, while magnetic fields do not, no matter how strong. This difference raises a pressing question. What is the nonequilibrium steady state induced by strong electric fields?
In this work, we develop a full nonequilibrium quantum transport theory of the chiral anomaly by using the Keldysh nonequilibrium Green function formalism [25] in conjunction with the Lindblad quantum master equation [26]. As a result, it is shown that there is a rich structure in the chiral anomaly transport, including the negative MR, the non-Ohmic behavior, the Esaki-Tsu peak, and finally the EQO. Being the incipient nonlinear behavior characterizing the chiral anomaly transport, the non-Ohmic behavior observed in BiSb alloys [18] can be regarded as a precursor to the occurrence of EQO, which can serve as the unmistakable quantum signature of the chiral anomaly in Weyl semimetals. We emphasize that the chiral anomaly provides a unique environment for the realization of WSL eigenstates in natural materials, which has been so far impossible except for synthetic systems such as semiconductor superlattices [27,28] and optical lattices [29,30].
In the perspective of application, this work lays a groundwork to expand the frontier of nonequilibrium quantum transport and realize novel nonlinear electronic devices by combining strong-field phenomena [31] with topological matter. It is interesting to mention that strong-field phenomena have been also investigated in combination with various many-body correlation effects such as Mott transition [26,[32][33][34][35][36][37][38] and many-body localization [39,40].

Landau-Stark quantization
To perform a concrete analysis of the chiral anomaly transport in the full quantum level, we consider a minimal tight-binding model for Weyl semimetals [41,42]: where h x (k) = 2t 1 sin k x , h y (k) = −2t 1 sin k y , and h z (k) = 2t 2 [2 − cos(k x − k y ) − cos(k x + k y )] + 2t 3 cos k z with t 1 , t 2 , t 3 being hopping amplitudes and σ x , σ y , σ z being Pauli matrices. With the time-reversal symmetry broken, this model hosts a single pair of Weyl nodes at k = (0, 0, ±π/2) with zero energy. See Fig. 1 a for the illustration of the energy spectrum at t 1 = 4, t 2 = 2, and t 3 = 1, which are to be used as hopping amplitudes throughout this work. Note that all momenta are denoted in units of corresponding inverse lattice constants. Also, unless specified otherwise, we set = c = 1 for simplicity. As elaborated later, in this work, we focus on half filling by setting the chemical potential appropriately. Let us first investigate what happens to the energy spectrum of the model Hamiltonian with the application of magnetic fields in the z direction (B = Bẑ with B > 0). Actually, the energy spectrum would develop a highly complicated fractal structure known as Hofstadter's butterfly, if magnetic fields are directly applied to the lattice model. To avoid this complication, we take the continuum limit of the model Hamiltonian within the x-y plane by replacing sin k i by k i and cos k i by 1 − k 2 i /2 for i = x, y in Eq. (1), while maintaining the full k z dispersion.
The application of magnetic fields can be implemented via minimal coupling, i.e., k → −i∇ + eA Landau with −e being the charge of electron and A Landau being the Landau-gauge vector potential. Consequently, the model Hamiltonian generates the following energy eigenvalues under magnetic fields: where n is a nonzero integer, and ω c = 4t 2 a x a y /l 2 B is the cyclotron frequency with a x and a y being the lattice constants in the x and y directions, respectively, and l B = 1/ √ eB being the magnetic length. For simplicity, the zero-point energy ω c /2 is subtracted from the energy eigenvalues to define n (k z ). Note that the energy eigenmodes are composed of the usual LL eigenstates, which are entirely dispersionless within the x-y plane, while dispersive in the z direction. Now, an interesting thing happens if one tries to set n = 0 in Eq. (2). With the sign of zero undefined, there could be two distinct energy eigenmodes corresponding to ±2t 3 | cos k z |. In reality, however, there exists only a single energy eigenmode called the chiral LL with the energy eigenvalue of 0 (k z ) = 2t 3 cos k z . Note that this singleness of the chiral LL is a unique topological property of Weyl semimetals. See Methods for details. Also, see Fig. 1 b for the illustration of chiral versus nonchiral LLs.
With the application of electric fields, each LL can be further quantized into a series of WSL eigenstates. Usually, the formation of WSL eigenstates requires wellseparated energy bands so that the Bloch oscillation can complete one full cycle without being interrupted by the Landau-Zener transition [23], which is unfortunately difficult to achieve in natural materials. Fortunately, in Weyl semimetals, there is a nice protection of the Bloch oscillation due to the aforementioned singleness of the chiral LL. Specifically, when electric fields are applied in the z direction parallel to magnetic fields (E = Eẑ), there is absolutely no Landau-Zener transition between different LLs, unless they have the same |n|. This means that the Bloch oscillation in the chiral LL is completely immune from the Landau-Zener transition. While allowed, the Landau-Zener transition is also energetically suppressed between nonchiral LLs with the same |n|, whose energies can be well-separated across the zero-point energy. Consequently, under strong magnetic fields, it is safe to assume that each LL is quantized into its own individual series of WSL eigenstates.
Also, being so-called extended states, LL eigenstates are generally known to be rather robust against disorder [43]. This means that the k z dispersion of LL eigenstates and consequently the formation of WSL eigenstates can be also robust against disorder to certain extents.
Technically, the application of electric fields can be implemented in terms of either static scalar or temporal vector potential gauge. In the static scalar potential gauge, the model Hamiltonian can be written in terms of the Stark Hamiltonian for each individual LL: where Ω = eEa z is the Bloch oscillation frequency with a z being the lattice constant in the z direction. Here, A n (k z ) is the Berry connection of the n-th LL, which turns out to be zero regardless of n in our minimal model. The Stark Hamiltonian can be diagonalized via WSL eigenstates in each LL, called Landau-Stark eigenstates, i.e., H Stark,n φ nl (k z ) = nl φ nl (k z ) with and nl =¯ n + lΩ, where¯ n = π −π dkz 2π n (k z ) is the mean energy of the n-th LL and l is the WSL index. It is important to note that WSL eigenstates are full quantum solutions of the Stark Hamiltonian, while also obtained as semiclassical solutions via the Bohr-Sommerfeld quantization [44]. See Fig. 1 c for the illustration of Landau-Stark energy levels, accompanied by the local density of states (DOS), whose details are given in Methods. Nonequilibrium quantum transport Being standing waves, WSL eigenstates cannot generate any nonzero net DC electric currents, unless there is impurity scattering, which causes the transition between different WSL eigenstates. Here, we develop a full nonequilibrium quantum transport theory of the chiral anomaly by treating the process of impurity scattering via the Keldysh nonequilibrium Green function formalism [25]. Specifically, our nonequilibrium quantum transport theory is composed of three steps. Temporal vector potential gauge. The first step is to change the gauge and implement the application of Our nonequilibrium quantum transport theory is based on the Keldysh-Dyson self-consistency loop comprising three parts; (i) the full Green function, G, is obtained by solving the Keldysh-Dyson equations with Σ being the yet-to-be-determined self-energy, (ii) Σ is then related with G via self-consistent Born approximation for impurity scattering, and (iii) the self-consistency loop is completed once the noninteracting Green function, g, is fixed in terms of Landau-Stark eigenstates. Crucially, the noninteracting lesser Green function, g < , is constructed so that Landau-Stark eigenstates are appropriately thermalized according to the WSL-wise thermalization scheme. electric fields via the temporal vector potential A Stark = −Etẑ with t being time, in which case the total vector potential is given as A = A Landau + A Stark . This particular choice of gauge is made to preserve the spatial translation symmetry so that impurity scattering can be treated via the usual method of self-consistent Born approximation (SCBA).
In the temporal vector potential gauge, the model Hamiltonian can be written as follows: which is periodic in time with the period of 2π/Ω. Such a time-dependent Hamiltonian can be analyzed by using the Keldysh nonequilibrium Green function method with nonequilibrium Green functions conveniently represented in the Floquet matrix form [33]. It is worthwhile to mention that WSL eigenstates in the static scalar potential gauge are manifested as Floquet modes in the temporal vector potential gauge.
Keldysh-Dyson self-consistency loop. The second step is to set up the Keldysh-Dyson self-consistency loop to capture the process of impurity scattering via SCBA. See Fig. 2 for the schematic diagram. Technically, the full Green functions can be obtained by self-consistently solving the Keldysh-Dyson equations [25]: where G r (g r ) and G < (g < ) are the full (noninteracting) retarded and lesser Green functions, which contain the information about the DOS and occupation, respectively. The advanced Green functions, G a and g a , are related with the retarded counterparts via complex conjugation. Meanwhile, Σ r and Σ < are the retarded and lesser selfenergies, respectively, induced by impurity scattering. In the above expressions, we drop all the subscripts (LL and Floquet indices) and arguments (k z and ω) for simplicity.
Importantly, the self-energies are related to the full Green functions via SCBA [45]: where V imp is the strength of the on-site Coulomb interaction between electron and impurity. The factor D = ω c /8πt 2 comes from the degeneracy of each LL. Now, the Keldysh-Dyson self-consistency loop is completed once the noninteracting Green functions, g r and g < , are fixed. The noninteracting retarded Green function, g r , is given in the Floquet matrix form as follows [26,33]: is the Fourier transform of Landau-Stark eigenstates in Eq. (4), and G r n (ε = ω + jΩ) = 1/(ε −¯ n + iη) is the reduced retarded Green function of Landau-Stark eigenstates. Note that the spectral information of Landau-Stark eigenstates is encoded via ϕ nl and G r n in g r n . See Methods for details. At this point, it is important to incorporate the broadening of Landau-Stark energy levels, which can be caused by any additional processes of inelastic scattering beyond SCBA. We implement such inelastic level broadening by setting η = Γ/2 in G r n with Γ being small, but finite. Specifically, we set Γ/t 3 = 0.05 throughout this work. Note that the DC electric current would be net zero in the presence of Γ alone [26]. Nonzero net DC electric currents can be only generated by the intricate interplay of both elastic and inelastic scattering.
The noninteracting lesser Green function, g < , is given in the Floquet matrix form as follows: where the reduced lesser Green function of Landau-Stark eigenstates, G < n (ε = ω + jΩ), is obtained via the WSLwise thermalization scheme, which is in turn derived as a solution of the Lindblad quantum master equation [26]. Specifically, each WSL eigenstate is individually thermalized according to the standard fluctuation-dissipation re-lation: where G a n (ε) = G r * n (ε) and f FD (ε) = 1/(e (ε−µ)/k B T + 1) is the usual Fermi-Dirac distribution function with the chemical potential µ set to be zero for half filling. See Methods for details.
DC electric current density. The third and final step is to compute the DC electric current density, J DC , from the full lesser Green function obtained as a converged solution of the Keldysh-Dyson self-consistency loop [26]: It is important to note that Eq. (12) itself is an exact expression of J DC , which means that J DC is accurate at arbitrary strengths of electric and magnetic fields if the full lesser Green function is so. See Methods for details.
In the following sections, we present the numerical results of J DC as a function of various parameters. Unless specified otherwise, all parameters with the energy unit (such as Ω, ω c , V imp , Γ, and so on) are denoted in units of t 3 throughout this work. Particularly, we set k B T /t 3 = 0.001 in this work. Also, LL indices are summed up to |n| = 6, which is necessary for the range of magnetic fields studied in this work, except for the ultra-quantum limit of strong magnetic fields, where it is sufficient to consider only the chiral LL. Meanwhile, the number of summed Floquet indices is chosen adaptively to ensure that J DC is well converged at each given Ω. See Methods for details.

Results
Electric quantum oscillation via the general Landau-Stark resonance. The DC electric current can oscillate via two different mechanisms. In this section, we first discuss the resonance between various Landau-Stark eigenstates with different LL indices, called the general Landau-Stark resonance. Fig. 3 a shows that, in a general regime of electric and magnetic fields, J DC oscillates as a function of both Ω and ω c , exhibiting a complicated, yet highly organized series of resonant peaks. Physically, the resonant behavior of J DC can be well understood in terms of the tunneling formula between adjacent sites [26]: where ρ loc and f loc are the local DOS and distribution function, respectively. Specifically, ρ loc is given as the sum of individual contributions from various LLs, i.e., Note that JDC exhibits a complicated, yet highly organized series of resonant peaks, whose trajectories are accurately described by the general Landau-Stark resonance condition,¯ n −¯ 0 = ∆l · Ω, for various cases of (n, ∆l). The strongest resonant peaks are obtained along the trajectories of (n = 1, ∆l = 1, 2, 3, · · · ) plotted in dark blue, followed by progressively weaker resonant peaks along those of n > 1 plotted in light blue. (b) Illustrated mechanism of the general Landau-Stark resonance. Here, we take the case of (n = 1, ∆l = 3) as an example, marked by the grey dot in a, where the Landau-Stark energy levels coming from the n = 1 LL (dark blue) are perfectly aligned with those from the n = 0, or chiral LL (red). ρ loc = D n ρ loc,n with where ε can cover the entire range of frequency by changing the Floquet index p while ω ∈ [−Ω/2, Ω/2]. Meanwhile, f loc can be computed via f loc = N loc /ρ loc with the local occupation number N loc given by N loc = D n N loc,n , where Note that Eq. (13) can be formally derived from Eq. (12) in the limit of strong electric fields, where WSL eigenstates form well-localized wave packets [26].
According to the tunneling formula, the DC electric current can be maximized if there is a large overlap between ρ loc (ε) and ρ loc (ε + Ω). Considering that ρ loc (ε) is composed of periodic peaks due to the Landau-Stark quantization, this means that the DC electric current can be maximized along the trajectories in the Ω-vs-ω c parameter space, satisfying the following condition of the general Landau-Stark resonance: where ∆l is an integer. See Fig. 3 b for the illustrated mechanism of the general Landau-Stark resonance. As seen from Fig. 3 a, the general Landau-Stark resonance condition describes the trajectories of maximized J DC quite accurately. Note that a similar resonance phenomenon has been observed in the transport experiment of semiconductor superlattices under parallel electric and magnetic fields [46].
Finally, to clearly show the periodicity of EQO, it is beneficial to plot J DC as a function of 1/Ω for various given ω c . Fig. 3 c shows that the resonant peaks are equally spaced as a function of 1/Ω with the period of 1/(¯ n −¯ 0 ), which is strongly reminiscent of the similar behavior in magnetic quantum oscillation. Actually, the low-electric-field data at ω c = 0.8 and 1 (red curves) reveals that, under strong magnetic fields, there is a new type of the EQO with different periodicity, which is shown below to be induced by a form of the self-resonance entirely within the chiral LL, called the chiral resonance. Electric quantum oscillation via the chiral resonance. The general Landau-Stark resonance condition can be trivially satisfied with n = 0 and ∆l = 0. If so, naïvely, J DC could be always enhanced in the ultra-quantum limit of strong magnetic fields, where the chiral LL becomes the only transport channel with all other nonchiral LLs pushed far away from the Fermi level. This naïve expectation, however, does not hold since the chiral LL alone cannot induce any actual electronic transport, at least via elastic impurity scattering alone. In this case, nonzero net DC electric current can be generated with help of the broadening of Landau-Stark energy levels due to inelastic scattering processes. Fig. 4 a shows the behavior of J DC as a function of Ω ranging from weak to strong electric fields in the ultraquantum limit of strong magnetic fields, say, at ω c = 1, where it is sufficient to consider only the chiral LL so long as Ω 0.6 (i.e., before the general Landau-Stark resonance comes into play). Particularly, in this limit, there are four distinct regimes of the chiral anomaly transport; (i) negative MR, (ii) non-Ohmic behavior, (iii) Esaki-Tsu peak, and (iv) EQO at weak (Ω 0.01), weak-to-intermediate (0.01 Ω 0.05), intermediate (Ω = Ω ET 0.05), and strong (Ω 0.05) electric fields, respectively.
First, at strong electric fields, the EQO occurs via the chiral resonance, which is distinguished from the previously described, general Landau-Stark resonance. In the case of the chiral resonance, the DC electric current oscillates as a function of 1/Ω with a constant period entirely independent of ω c , which is simply π/4 in units of 1/t 3 in our minimal model for Weyl semimetals. Fundamentally, the mechanism of the chiral resonance can be understood in terms of the wave function overlap between adjacent WSL eigenstates in the chiral LL, which oscillates asymptotically as a function of electric field. See Methods for details.
As Ω decreases, the EQO becomes less and less pronounced, finally merging into the Esaki-Tsu peak around Ω = Ω ET . Fig. 4 b shows that both Esaki-Tsu peak and subsequent EQO are closely correlated with the formation of well-separated WSL eigenstates. Note that, marking the onset of negative differential conductivity, the Esaki-Tsu peak [47] has been routinely observed in semiconductor superlattices [48]. Fig. 4 c shows that, at Ω Ω ET , J DC increases as a monotonic, but in general nonlinear function of Ω, i.e., J DC = σΩ + σ Ω 2 + σ Ω 3 + · · · , where σ denotes the usual linear Drude conductivity in the limit of weak electric fields, while σ and σ are the two lowest-order coefficients of the non-Ohmic behavior. It is important to note that the non-Ohmic behavior is an inevitable crossover phenomenon connecting between the linear Drude conductivity and Esaki-Tsu peak. Considering that both Esaki-Tsu peak and subsequent EQO are closely correlated with the formation of well-separated WSL eigenstates, the non-Ohmic behavior can be regarded as a precursor to the EQO. In this context, the non-Ohmic behavior observed in BiSb alloys [18] suggests that the observation of EQO might actually be within the reach of experiments since, in our results, the strength of electric field necessary for the occurrence of EQO is only about 10-20 times larger than that necessary for the non-Ohmic behavior. Now, we would like to confirm if the linear Drude conductivity, σ, exhibits the expected behavior of negative MR. Specifically, in the ultra-quantum limit of strong magnetic fields, σ is expected to increase as a linear function of magnetic field [3,6]. Fig. 4 d confirms that this is indeed exactly the case. Actually, Fig. 4 e shows that J DC increases as a whole with stronger magnetic fields. Finally, Fig. 4 f shows the behavior of J DC as a function of Ω for various given V imp , confirming that σ decreases with stronger impurity scattering, as expected from the Drude behavior.

Discussion
In this work, it is shown that the chiral charge pumping is essentially nothing but the Bloch oscillation. Both topologically and energetically protected in Weyl semimetals, the Bloch oscillation can be quantized to generate robust Landau-Stark eigenstates, eventually giving rise to the resonant oscillation of the DC electric current as a function of electric field.
Called the EQO, this resonant oscillation of the DC electric current can occur in Weyl semimetals via two different mechanisms. First, the EQO can occur via the resonance between various Landau-Stark eigenstates with different LL indices. Second, in the ultra-quantum limit of strong magnetic fields, the EQO can also occur via a form of the self-resonance within the chiral LL. Particularly, in this limit, there are four distinct regimes of the chiral anomaly transport; (i) negative MR, (ii) non-Ohmic behavior, (iii) Esaki-Tsu peak, and (iv) EQO at weak, weak-to-intermediate, intermediate, and strong electric fields, respectively. It is important to note that both negative MR and non-Ohmic behavior [18] have been already observed in Weyl semimetals, providing experimental support for the occurrence of EQO in natural materials.
In broad perspective, understanding nonequilibrium steady states of matter is among the foremost frontiers in physics. Induced by strong electric fields, the EQO would be one of the most salient features of nonequilibrium steady states realized in condensed matter. Usually achieved in synthetic systems such as semiconductor superlattices and optical lattices, a prerequisite for the occurrence of EQO is the formation of robust WSL eigenstates. As emphasized in this work, the chiral anomaly can provide a unique environment for the formation of robust WSL eigenstates via the combination of strongfield phenomena with topological matter. Interestingly, Weyl semimetals can be also synthetically generated by fabricating a layered structure of alternating topological and magnetic insulators [5].
Finally, there is a close analogy between the EQO studied in this work and the radiation-induced quantum oscillation observed in quantum Hall systems [49,50]. It is interesting to mention that the radiation-induced quantum oscillation has been analyzed via both Keldysh nonequilibrium Green function method and tunneling formula [51][52][53], which are also two main theoretical tools in this work. Methods Landau quantization in Weyl semimetals. We begin by writing the continuum limit of the model Hamiltonian in Eq. (1) within the x-y plane, which can be obtained by replacing sin k i by k i and cos k i by 1 − k 2 i /2 for i = x, y, while maintaining the full k z dispersion. Specifically, the model Hamiltonian can be written in the continuum limit as follows: where all momenta are denoted in units of corresponding inverse lattice constants.
With the application of magnetic fields in the z direction, the model Hamiltonian is modified via minimal coupling, i.e., k → Π = −i∇ + eA Landau with A Landau = B(0, x, 0) being the Landau-gauge vector potential. At this moment, let us assume that B > 0. The case of B < 0 is to be considered separately below. For B > 0, the model Hamiltonian can be written as where the LL raising and lowering operators, b † and b, are defined, respectively, as follows: with l B = 1/ √ eB being the magnetic length. Similarly, the pseudospin raising and lowering operators, σ + and σ − , are defined, respectively, as follows: Note that the cyclotron frequency is given by ω c = 4t 2 /l 2 B , and Π z is replaced back to its eigenvalue, k z . The Hamiltonian in Eq. (18) can be block-diagonalized by using the convenient set of basis states, {|ν ⊗ |σ }, which are composed of number eigenstates |ν (i.e., b † b|ν = ν|ν ) and the pseudospin up/down state |σ (i.e., | ↑ or | ↓ ). Now, by noting that one can obtain the block-diagonalized matrix form of the Hamiltonian as follows: which is defined in the Hilbert space spanned by two basis states, |ν ⊗ | ↑ and |ν − 1 ⊗ | ↓ with ν ≥ 1. Diagonalizing H ν (k z ) generates the energy eigenvalues of nonchiral LLs as follows: which becomes identical to n (k z ) in Eq. (2) after the LL index is defined as n = ±ν, and the zero-point energy ω c /2 is subtracted.
Meanwhile, the Hamiltonian is already fully diagonalized for ν = 0: which is defined in the Hilbert space spanned by the single basis state, |0 ⊗ | ↑ . Being diagonal, H 0 (k z ) itself is the energy eigenvalue of the chiral LL, which equals to 0 (k z ) after the subtraction of the zero-point energy.
It is important to note that the singleness of the chiral LL is a unique topological property of Weyl semimetals. To appreciate the origin of this topological property, it is beneficial to consider what happens in the case of B < 0. Actually, the model Hamiltonian can be written for the general sign of B as follows: where the LL raising and lowering operators are now generalized as follows: with l B = 1/ e|B|. After some algebra, one can show that the energy eigenvalues of nonchiral LLs are exactly the same as before regardless of the sign of B except that the zero-point energy is now generalized as sgn(B)ω c /2.
The situation is quite different for the chiral LL. That is, unlike those of nochiral LLs, the energy eigenvalue of the chiral LL depends on the sign of B: 0 (k z ) = 2t 3 sgn(B) cos k z . This sign dependence of the chiral LL is fundamentally due to the specific topological property of Weyl semimetals in our minimal model. Namely, the 2D k z slices of the Brillouin zone form Chern or trivial insulators depending on whether k z is inside or outside the region between two Weyl nodes with opposite chiralities.
Noninteracting Green functions in the Floquet matrix form. Here, we discuss how to construct the noninteracting retarded and lesser Green functions in the Floquet matrix form. We begin by writing the noninteracting Hamiltonian in the temporal vector potential gauge as follows: where c † n,kz and c n,kz are the creation and annihilation operators, respectively, for the n-th LL with k z .
The noninteracting retarded Green function is defined as follows: where c n,kz (t) = U n,kz (t, t 0 )c n,kz (t 0 ) with the unitary evolution operator, U n,kz (t, t 0 ), given by where t 0 is some arbitrary reference time.
Now, noting that the unitary evolution operator in Eq. (29) is essentially identical to the wave function of Landau-Stark eigenstates in Eq. (4), U n,kz (t, t 0 ) can be expressed in terms of φ nl as follows: where l can be chosen arbitrarily. Then, plugging Eq. (30) into the anticommutation part in Eq. (28) leads to the following result: where it is used that {c n,kz (t 0 ), c † n,kz (t 0 )} = 1 and φ nl (k z − Ωt 0 )φ * nl (k z − Ωt 0 ) = 1. Next, by using the integral representation of the Heaviside step function, one can express g r n,kz (t, t ) as follows: which is obtained after an appropriate redefinition of the integration valuable.
The mathematical form of Eq. (35) indicates that [g r n (k z , ω)] pq is nothing but the Fourier transform of g r n,kz (t, t ). Specifically, [g r n (k z , ω)] pq is the (p, q)-th element of the noninteracting retarded Green function in the Floquet matrix form [33].
Based on this realization, it is instructive to compute the noninteracting local DOS, ρ which shows that the local DOS is composed of discrete peaks at ε = nl with their weights given by the corresponding Landau-Stark eigenstates at a given site, say, origin, |ϕ nl (0)| 2 . Note that the broadening of Landau-Stark energy levels can be implemented by setting η = Γ/2 in G r n (ε) with Γ being small, but finite, in which case the delta function is replaced by the Lorentzian: with Γ quantifying the broadening width. Now, let us switch gears to the noninteracting lesser Green function. Actually, Eq. (36) suggests a very natural mathematical expression for the noninteracting lesser Green function: [g < n (k z , ω)] pq = e ikz(p−q) j ϕ np (j)G < n (ω + jΩ)ϕ * nq (j) (39) where G < n (ε = ω + jΩ) is the reduced lesser Green func-tion of Landau-Stark eigenstates. As mentioned in the main text, G < n (ε) is obtained via the WSL-wise thermalization scheme [26]. Specifically, each WSL eigenstate is individually thermalized according to the standard fluctuation-dissipation relation: where G a n (ε) = G r * n (ε) and f FD (ε) = 1/(e (ε−µ)/k B T + 1) is the usual Fermi-Dirac distribution function with the chemical potential µ set to be zero for half filling.
As done before, it is also instructive to compute the noninteracting local occupation number, N which shows that each WSL eigenstate is indeed individually thermalized with its own effective chemical potential, µ eff (l) = µ + lΩ. DC electric current density from the full lesser Green function. The DC electric current can be computed from the full lesser Green function obtained as a converged solution of the Keldysh-Dyson self-consistency loop.
To begin with, whether DC or not, the electric current density can be exactly expressed in terms of the full lesser Green function as follows: which can be understood as the sum of all contributions from each conduction mode, whose individual contribution is in turn given by the product between its group velocity and occupation number specified by the LL index n and momentum k z . It is important to note that the above expression is in principle exact at arbitrary strengths of electric and magnetic fields. As shown below, eventually, the electric current becomes strictly DC in our situation.
First, by definition, the occupation number is equal to the equal-time full lesser Green function, which can be related to its Fourier transform as follows: where [G < n (k z , ω)] pq is the (p, q)-th element of the full lesser Green function in the Floquet matrix form.
Next, the k z integration in Eq. (42) can be explicitly performed by using the constraint that [G < n (k z , ω)] pq = e ikz(p−q) [G < n (k z = 0, ω)] pq , ensuring that k z always appears as the particular form of k z − Ωt. This constraint is a manifestation of the gauge invariance in our situation. With help of Eq. (44), one can then obtain the final expression for the electric current density: −Ω/2 dω 2π n,p,q p¯ n (p) G < n (k z = 0, ω) p+q,q , where¯ n (p) = π −π dkz 2π e ikzp n (k z ). Note that the nominal time dependence in the initial expression disappears in the final expression after the k z integration. Consequently, as mentioned before, the electric current becomes strictly DC. Truncation of Floquet matrices. In the Floquet representation, Green functions are represented as infinitedimensional Floquet matrices. For practical calculations, the dimension of Floquet matrices should be truncated with an appropriate cutoff limiting the range of Floquet indices. In other words, we would like to represent retarded and lesser Green functions as finite-dimensional Floquet matrices, [G r (ω)] pq and [G < (ω)] pq , respectively, with p, q ∈ (0, ±1, · · · , ±L). The cutoff L is determined via the following procedure.
To begin with, we first estimate the cutoff by requiring that the noninteracting local DOS is properly normalized for each individual LL. Specifically, it can be said that the noninteracting local DOS for the n-th LL is properly where ρ (0) loc,n is given in Eq. (37), L n is the cutoff for the n-th LL, and δ tol is a sufficiently small tolerance. In this work, we set δ tol to be 10 −8 .
Finally, the overall cutoff L is chosen as the maximum of L n : L = max{L n }. As a general rule, the lower Ω becomes, the higher L is required. Roughly speaking, L is of the order of 1, 000 for Ω 0.01 while typically less than 100 otherwise. Mechanism of the chiral resonance. To understand the mechanism of the chiral resonance, we begin by rewriting the tunneling formula as follows: dε [ρ loc (ε + Ω)N loc (ε) − ρ loc (ε)N loc (ε + Ω)] , (47) where ρ loc and N loc are the local DOS and occupation number, respectively. Note that Eq. (47) is precisely identical to Eq. (13) since N loc = ρ loc f loc by definition. Now, assuming that WSL eigenstates are well separated in the chiral LL, the local DOS can be accurately approximated as where A l (Ω) = |J l (2t 3 /Ω)| 2 and δ Γ (ε) is the Lorentzian in Eq. (38). Similarly, the local occupation number can be also accurately approximated as which is obtained via the WSL-wise thermalization scheme as explained in Eq. (41).
After some rearrangements, Eq. (47) can be rewritten as follows: where F m (Ω) = l A l+1 (Ω)A l+m (Ω) and which is a monotonic function of Ω without any oscillatory behaviors. This means that, if any, oscillatory behaviors should come from F m (Ω), which depends on the wave function form of WSL eigenstates in the chiral LL.
In the chiral LL, WSL eigenstates are described by the Bessel function, which can be approximated as for x |l 2 − 1/4|. After some algebra making use of this asymptotic behavior of the Bessel function, one can show that F m (Ω) is an oscillatory function of 1/Ω with the period of π/4 in units of 1/t 3 . In conclusion, the chiral resonance is due to the asymptotic, oscillatory behavior of WSL eigenstate wave functions in the chiral LL.