Identifying Axion Insulator by Quantized Magnetoelectric Effect in Antiferromagnetic ${\mathrm{MnBi}}_{2}{\mathrm{Te}}_{4}$ Tunnel Junction

Intrinsic magnetic topological insulator ${\mathrm{MnBi}}_{2}{\mathrm{Te}}_{4}$ is believed to be an axion insulator in its antiferromagnetic ground state. However, direct identification of axion insulators remains experimentally elusive because the observed vanishing Hall resistance, while indicating the onset of the axion field, is inadequate to distinguish the system from a trivial normal insulator. Using numerical Green's functions, we theoretically demonstrate the quantized magnetoelectric current in a tunnel junction of atomically thin ${\mathrm{MnBi}}_{2}{\mathrm{Te}}_{4}$ sandwiched between two contacts, which is a smoking-gun signal that unambiguously confirms antiferromagnetic ${\mathrm{MnBi}}_{2}{\mathrm{Te}}_{4}$ to be an axion insulator. Our predictions can be verified directly by experiments.

Unlike other topological phases characterized by the first Chern number [12,13], an axion insulator is in a higher order topological phase characterized by the symmetry protected axion field θ = π [14][15][16][17], which can manifest as the quantized topological magnetoelectric (TME) effect [18][19][20] and other striking transport phenomena [21][22][23][24][25]. However, because the first Chern number of an axion insulator vanishes identically, the ensuing transport effect on a Hall bar device exhibits a vanishing Hall resistance accompanied by a large longitudinal resistance, which is just similar to a normal insulator. This property makes it rather elusive to properly distinguish axion insulators from normal insulators by transport experiments [11]. Therefore, to confirm the existence of axion insulator, a viable experimental scheme without ambiguity is needed.
In this Letter, we propose an axion insulator tunnel junction consisting of a few-layer MBT sandwiched between two metallic contacts as an experimental setup to unambiguously identify axion insulators through the quantized TME. We first show that a perpendicular magnetic field can induce a surface charge polarization that is physically related to the layer-resolved Chern numbers and the quantized axion field θ = π. When the magnetic field adiabatically varies with time (i.e., with a frequency far less than the insulating gap), the surface charge polarization becomes time-dependent and will generate an AC charge current through the tunnel junction. We use timedependent non-equilibrium Green's function to quantity the detectable AC current driven by a harmonic magnetic field, which agrees remarkably well with the time derivative of the induced charge polarization, thus strengthening the validity and reliability of our proposed scheme to identify axion insulators. Since archetypal materials parameters have been adopted in the modeling, we anticipate our theory to be able to inspire and guide experiments in the perceivable future.
, σ the spin-σ orbital of Bi (Te), the low energy Hamiltonian for a layered MBT can be written as [4,26,27] Here, the first term is an effective Hamiltonian for a three dimensional topological insulator, where d 0 (k) = M 0 − B 1 k 2 z − B 2 (k 2 x + k 2 y ), d 1 (k) = A 2 k x , d 2 (k) = A 2 k y , d 3 (k) = A 1 k z with A 1 (2) , B 1 (2) , M 0 being system parameters and the lattice momentum k = k x , k y , k z . Γ 0 = σ 0 ⊗τ 3 and Γ a = σ a ⊗τ 1 (a = 1, 2, 3) where σ a and τ a are Pauli matrices acting on the spin and orbital spaces, respectively. The second term describes the exchange interaction between topological electrons and magnetic ordering, where ∆ is the exchange strength and m i is the unit magnetization vector of the i-th SL [4]. Henceforth in all numerical calculations, the materials parameters are shown in Tab. I, and temperature is set to be zero.  (2) and B 1 (2) are based in Ref. [26]. J, K and MS are chosen from Ref. [28]. µB is the Bohr magneton. The exchange gap ∆ is evaluated from Refs. [29,30].
Since the topological states of MBT is intertwined with the magnetic ordering, we first need to determine its magnetic configuration. In the macro-spin approximation (spatially uniform magnetization within a particular SL), the magnetic property of an N -SL MBT can be characterized by the free energy [31] where J is the antiferromagnetic interlayer exchange interaction, K is the easy-axis anisotropy, B is external magnetic field, and M s is the saturation magnetization of each SL. The magnetization vector is parameterized as Without losing generality, we assume that B is applied along z direction and m i rotates only in the xz plane. We obtain the equilibrium magnetic configuration by minimizing the free energy U using the steepest descent method [32], which is detailed in the supplementary materials (SM) [33]. Figure. 1(a) shows the total magnetization as a function of the applied magnetic field for a six-SL MBT, where we identify the spin-flop critical points at around B ± c ≈ ±3T, beyond which the Zeeman energy overcomes the exchange and anisotropy interactions and induces non-collinear spin configurations until the system is fully polarized into a ferromagnetic state at above 10T (see Fig. S1 in the SI). Such distinct magnetic evolution is in quantitative agreement with experiments [2,11]. The complicated spin configurations in the intermediate spinflop phases are discussed in the SM [33].
In-plane transport properties on a Hall bar.-To study the electronic transport, we first discretize the continuum Hamiltonian Eq. (1) on a cubic lattice (a 0 = 5nm) invoking the k · p perturbation. Then, under a Hall bar device geometry as illustrated in the inset of Fig. 1(a), we calculate the Hall resistivity ρ xy and the longitudinal resistivity ρ xx using the Landauer-Büttiker formula [33,34]. To incorporate fluctuations, we also add a disorder potential H D = V (r)s z ⊗ σ 0 to the lattice Hamiltonian, where V (r) is uniformly distributed within [−D/2, D/2] with D being the disorder strength. The Fermi level is zero as we do not consider doping or gating.  For a six-SL MBT device reflecting real experimental setup [11], we obtain ρ xx and ρ xy by averaging 160 repeated calculations, which are plotted as functions of magnetic field B (along z) in Fig. 1(b). The results show a topological phase transition from a normal insulator (indistinguishable from an axion insulator) with a vanishing Chern number C = 0 at low magnetic fields into a quantum anomalous Hall insulator with C = ±1 at high magnetic fields. When |B| < B + c , the magnetic ground state remains antiferromagnetic with antiparallel spins on adjacent SL, and the system preserves the PT symmetry. Because the spin flips its sign under PT operation [PT : H(k, ↑) → H(k, ↓)], the bands must be doubly degenerate with a band gap of δ ≈ 2∆ at k x = 0, as shown in the left inset in Fig. 1(b). Consequently, we obtain C = 0, hence a vanishing Hall resistivity and a large longitudinal resistivity akin to a normal insulator. While ARPES experiments showed controversial results on the band gap in MBT [29,35], transport measurements strongly support the existence of large gaps in both the antiferromagnetic and ferromagnetic states of MBT [2,10,11] by confirming the insulating behavior in longitudinal transport, even though this insulating gap cannot tell axion insulators from normal insulators.
When |B| exceeds B + c , however, the magnetic moments undergo a spin-flop transition which breaks the time reversal symmetry for electrons. Correspondingly, the topological Chern number becomes C = −sgn(B), leading to a quantized Hall resistivity ρ xy = h/(Ce 2 ) and a vanishing longitudinal resistivity ρ xx = 0 [12]. The deviations of ρ xy around integer values are ascribed to the finite size effect, which can be suppressed by enlarging the system size. The in-plane resistivities shown in Fig. 1 agree quantitatively with experimental observations [11? ] widely regarded as evidences for axion insulator. Nevertheless, the topological phase transition taking place here is inadequate to determine an axion insulator because the C = 0 phase appearing at small fields by itself is indistinguishable from a normal insulator.
Surface charge polarization and layer-resolved Chern numbers.-A defining feature of axion insulator is the topological TME enabled by the quantized θ-field, which, unlike the Chern number C, can uniquely characterize the axion insulator phase. On the one hand, a magnetic field B below the spin-flop threshold will induce a quantized charge polarization P = e 2 θB/(2πh) [14], which is intimately related to the layer-resolved Chern numbers. If the applied B field is time dependent, a charge current proportional to dP /dt will be generated, enabling a directly detectable signal to be discussed later. On the other hand, the TME also manifests as the magnetization induced by an electric field [37]. However, the TME coefficient quantized by θ is typically two orders of magnitude smaller than that of ordinary magnetoelectric materials [38]. Therefore, the TME is more amenable to transport measurement as the sensitivity of detecting current is extremely high. Nonetheless, as a consistency check, we also calculated the tiny magnetization induced by an electric field, which indeed turns out to be quantized by the θ field (see the SM [33]).
To calculate P , we consider a slab of thickness L z and widths L x = L y with open boundary conditions and assume that a static magnetic field B = (0, 0, B) is applied along the z-direction, which amounts to a magnetic flux of Φ 0 = Ba 2 0 per unit cell. Using the equilibrium Green's function method [33], we obtain the charge distribution Q(r) = −e n(r) where −e is the electron charge andn(r) is the electron density operator. Figure 2 (blue dots) plots the charge distribution among each SL, Q z = x,y Q(r), with respect to an averaged background charge Q(r) = x,y,z Q(r)/L z which compensates the positive ions in the lattice. Since Q(r) is an odd function of z, as shown in Fig. 2, there is indeed a finite charge polarization P = dV rQ(r). As will be shown later, only surface charges contribute to the detectable current, thus only the surface charge polariza- is relevant to our discussion. Ideally, the surface charge polarization P z should be very close the total polarization P , but finite-size effects can bring about deviations. Fortunately, we find that the finite-size effects are well suppressed by increasing the thickness L z . It turns out that P z (L z = 6, 8, 10) = 0.91, 0.97, 0.99 (e 2 Φ/2h) with Φ = L x L y Φ 0 being the total magnetic flux penetrating the slab, rapidly approaching the quantized value determined by the axion field θ = π.
We now turn to the layer-resolved Chern numbers C z which reflect the relative contribution to the system topology by different SLs. To this end, we adopt periodic boundary conditions in the lateral dimensions under the same slab geometry used above. While the layerresolved Chern numbers can be straightforwardly obtained by projecting the wavefunctions onto each SL [33] in a clean system, here we resort to the non-commutative approach which is able to incorporate disorders [39]: C z = −2πiTr P [x,P ], [ŷ,P ] P z , wherex (ŷ) is the position operator,P is the projector onto the occupied bands,P z ≡ |ψ z ψ z | is the projector onto the z-th SL, [· · · ] is the commutator and Tr denotes the trace. In the presence of PT symmetry,P z flips sign on opposite surfaces becauseP −z = |ψ −z ψ −z | = PT |ψ z ψ z |PT = (PT ) 2 |ψ z ψ z | = −|ψ z ψ z | = −P z , ensuring that the layer-resolved Chern numbers are odd in z/L z . Figure 2 shows the layer-resolved Chern numbers C z with three different thicknesses L z = 6, 8, 10 (red squares), which agrees remarkably well with the charge distribution Q z . Even in the presence of disorders, we find that C z and Q z are very robust (See SM [33]), suggesting that they are topologically protected properties intrinsic to the axion insulator. Correspondingly, the surface Chern number C Lz surf = 0 z=−Lz/2 C z is almost half quantized: C 6 surf = 0.46, C 8 surf = 0.48, and C 10 surf = 0.49, indicating a distinct bulk axion field θ = π [40].
Charge current in MBT tunnel junction.-To detect the quantized TME in MBT using transport experiment, we need to consider a time-dependent magnetic field such that the induced surface charge polarization becomes dynamical and produces a charge current in the z direction. This approach has been utilized to characterize multiferroic materials exhibiting non-quantized magnetoelectric effects [38]. To this end, we conceive an axion insulator tunnel junction (AITJ) consisting of an even-SL MBT sandwiched between two metallic contacts [41] as illustrated in Fig. 3(a). In the adiabatic limit, the charge polarization follows the magnetic field at any instant of time, which can be detected directly as an AC output signal from the AITJ. Since the metallic contacts are connected only to the top and bottom layers, only the surface polarization P z is relevant to the transport measurement.
Using the lattice Hamiltonian, we resort to the timedependent non-equilibrium Green's function to compute the output current in the AITJ [33]. A harmonic magnetic field B(t) =ẑB sin ωt applied to the AITJ converts to a phase Φ 0 (t) = Φ 0 sin ωt for electrons, where Φ 0 = Ba 2 0 is the magnetic flux per unit cell. As a result, the effective Hamiltonian acquires a time-dependent perturbation that drives the electron motion, forming a charge current. Since the system is now periodic in time, the induced charge current can be expanded into a Fourier series as [33,42] where I n is n-th harmonic component satisfing I n = −I * −n , ensuring a real current. The total current I(t) includes a DC component I 0 and a series of AC components I n>0 . Truncating the Green's function at order n = 4 suffices to yield a converging result [33]. As an independent confirmation, we use the same harmonic field B(t) =ẑB sin ωt in the surface charge polarization P z and calculate the resulting charge current I z (t) = dP z (t)/dt [44], assuming an adiabatic condition that P undergoes a quasi-static variation without interband transitions induced by the oscillating B(t) [45]. Figure 3(b) (solid red curve) plots I z (t) within one period of oscillation for a same MBT slab, which agrees remarkably well with the harmonic signal I(t) obtained by the non-equilibrium Green's function method. To benchmark the accuracy of our numerical results, we also plot the ideal case for an infinite system, where I(t) = dP (t)/dt = θe 2 /(2πh)ωΦ cos ωt with a strictly quantized axion field θ = π (dotted black curve in Fig. 3(b)). We see that our numerical results obtained both from the non-equilibrium Green's functions and from P z (t) only slightly deviate from the ideal case, which demonstrates the validity and reliability of our proposal. We mention in passing that if the Fermi level is tuned into the conducting band (e.g., by gating the device [46]), the MBT will become metallic and the induced AC current will vanish.
For an MBT of size L x = L y = 10µm, a harmonic magnetic field of strength B ∼ 100Gs and frequency ω/2π = 1GHz induces an output AC current I ∼ 121.54nA, which is a conservative estimation. Since I scales as ωBL x L y , the output current can be amplified by increasing the driving frequency ω, the magnetic field B, or the system size in the lateral dimensions. In the ideal case, the induced surface charge polarization P z should scale linearly with the magnetic flux per unit cell Φ 0 . To evaluate potential deviations due to finite-size effects, we plot P z as a function of Φ 0 for different thicknesses against the ideal scaling in the inset of Fig. 3(b), where the finite-size effects turn out to be negligible, further confirming the validity of our calculations.
In summary, we have theoretically proposed an experimental setup to unambiguously identify antiferromagnetic MBT as an axion insulator by detecting the AC current induced by a harmonic magnetic field under the adiabatic condition. Comparing to the vanishing Hall resistance measured in previous experiments, which is inadequate to confirm the axion insulator phase, our proposed scheme provides a smoking-gun signal to identify MBT as an axion insulator.

I. LATTICE HAMILTONIAN
We use the k · p theory to discretize the effective Hamiltonian for MnBi 2 Te 4 (MBT) on a cubic lattice the annihilating operator for an electron with orbital s (p) spin σ on site i = x, y, z . Substituting the momentum k α=x,y, with a 0 the lattice constant, we can write the lattice Hamiltonian as where Here, T i is the onsite term. T α=x,y,z is the hoping term along the α direction. T zM is the exchange coupling between topological electrons and magnetic moments. ∆ is the exchange gap. m z = sin θ z , 0, cos θ z with θ z the polar angle of the magnetic moment in z-th septuple layer, which is assumed to be localized inside the xz plane for simplicity. T D represents magnetic disorders, where D(i) is a random number within [−D/2, D/2] with D the disorder strength.
Other parameters are defined in the main text.
In the presence of a static magnetic field B = 0, 0, B , the lattice Hamiltonian acquires a phase factor φ ij = 2π j i A · dr/φ 0 where φ 0 = h/(2e) is the magnetic flux quanta and A is the vector potential (∇×A = B) Here, we adopt the Landau gauge A = yB, 0, 0 . Using the Peierls substitution, we can recasts the lattice Hamiltonian into where f x (y) = e iφ ii+x = e iyBa 2 0 /φ0 . Under a time dependent magnetic field B = 0, 0, B(t) , the vector potential A(t) becomes time dependent, so is the f x (y, t) term. In particular, if the magnetic field is harmonic, i.e., B(t) = B 0 sin ωt, the lattice Hamiltonian is also periodic H(t) = H(t + T ) with T = 2π/ω, which can be written as H(t) = H 0 + H(t) on a slab of size L x × L y × L z , where H 0 is the time independent component while H(t) = i (ψ † i T x e ipB0a 2 0 sin ωt ψ i+x + h.c.). Using the expansion exp [iz sin (α)] = ∞ k=−∞ J k (z) exp (ikα) [1], we can finally rewrite the lattice Hamiltonian as where J k (±yΦ 0 ) is Bessel function of the first kind at order k, and Φ 0 = B 0 a 2 0 is the magnetic flux per unit cell. When setting Φ 0 = 0 (removing the magnetic field), Equation (S3) reduces to Eq. (S2) because J k (0) = δ k and k J k (x = 0) = 1.  Table I in the main text.

II. MAGNETIC CONFIGURATION
The magnetic properties of a N-septuple layer MBT can be captured by the free energy shown in Eq. 2 in the main text. Minimizing the free energy gives the equilibrium magnetic configurations. Figure S1 (a) shows the magnetic configuration and the total magnetization as functions of external magnetic field B for a six-septuple layers= MBT. In general, the system has seven distinct magnetic states as the magnetic field is tuned continuously from −11 Tesla to 11 Tesla. When |B| < B + c1 , the Zeeman energy is not large enough to overcome the anisotropy. Consequently, the system stays in its antiferromagnetic ground state with a vanishing magnetization. Namely, the spins in adjacent layers are antiparallel to each other, preserving the combined parity and time reversal (PT ) symmetry. For B + c1 < |B| < B + c2 , the system turns into a canted state, in which the canting angles among different layers are somehow random. Meanwhile, the magnetization exhibits a sudden change at the critical point B ± c1 where the Zeeman energy exceeds the magnetic anisotropy. For B + c2 < |B| < B + c3 , the system enters a spin flop state with the canting angle of layer z and layer −z satisfying θ z + θ −z = 2π. The system finally becomes ferromagnetic with all spins collinear with the magnetic field when |B| > B + c3 . Figure. S1(b) schematically illustrates the seven different spin configurations taking place at different magnetic fields.

III. LANDAUER BÜTTIKER FORMULISM
With the discretized Hamiltonian, we can calculate the Hall resistivity and longitudinal resistivity in a Hall-bar device geometry. The Landauer Büttiker formula defines the transmission coefficient from lead n to lead m as where G r,a = (−H − 6 m=1 Σ r,a m ) −1 are the retarded and advance Green's functions. Γ m = i(Σ r m − Σ a m ), with Σ r,a m the self energy, is the line width function. The potential V m of lead m can then be determined by solving  where T ii = − j =i T ij , and the current vector is set to beÎ = I, 0, 0, −I, 0, 0 . The Hall resistivity ρ xy = (V 2 − V 6 )/I and the longitudinal resistivity ρ xx = (V 2 − V 3 )/I in the main text are obtained by averaging 160 repeated numerical calculations.

IV. LAYER-RESOLVED CHERN NUMBERS USING THE TKNN FORMULA
The layer-resolved Chern numbers can be alternatively obtained using the TKNN formula [2] where F = 0 as the Fermi energy is inside the surface gap, m(n) (k) and |m (|n ) are the eigenenergy and eigenstate of H,P z ≡ |ψ z ψ z | is the projector onto z-th layer. Figure. S2 shows the layer-resolved Chern numbers calculated using Eq. S6 and the non-commutative formula for three different thicknesses L z = 6, 8, 10 of width L x = L y = 40. The two approaches agree remarkably well with each other.

V. CHARGE DENSITY
In the presence of a static magnetic field B = (0, 0, B), the charge density can be expressed as where −e is the electron charge,n(r) = ψ † (r, t)ψ(r, t) is the number operator, G < (r, t; r , t) is the distribution Green's function and Tr represents the trace of a matrix. After a Fourier transform, we have  Here, the blue (red) dots represent positive (negative) charges, the size of which reflects the local charge density.

VI. DISORDER EFFECT
When fabricating the experimental setup proposed in the main text, strong magnetic disorders may be present. It is thus imperative to study to what extend the magnetic disorders can affect the axion insulator state of MBT. To this end, we calculate the layer-resolved Chern number C n and layer-resolved polarized charge P z under magnetic disorders D = 0.1eV, which is strong enough to collapse the surface band gap. As a comparison, we also consider the case without disorders, and show the results together in Fig. S4. Overall, both the layer-resolved Chern numbers and the polarization remain almost unchanged in the presence of strong disorders. Specifically, the surface Chern numbers are C 6 surf (D = 0) = 0.4524 and C 6 surf (D = 0.1) = 0.4547, which are very close to each other and almost invisible in Fig. S4(a). The surface polarization is P z (D = 0) = 0.9362e 2 Φ/2h and P z (D = 0.1) = 0.9713e 2 Φ/2h, where the disorder-induced change is only 3% − 4%. These results clearly demonstrate that the axion insulator phase is robust against magnetic disorders, which is consistent with previous studies [3,4].

VII. ELECTRIC FIELD INDUCED MAGNETIZATION
The magnetoelectric response manifests in two complementary aspects. While in the main text we focus on the charge polarization induced by a magnetic field, here, for consistency, we check the magnetization induced by an electric field. For E = (0, 0, E), the Fermi level experiences a layer-dependent perturbation. The Fermi energy of layer z is thus µ F (z) = (z − 1)/L z δµ F where δµ F = −eEL z is the total Fermi energy drop. Since this electric field does not break the in-plane transnational symmetry, we can thus resort to the following equation to calculate the magnetization in a unit cell [5,6] M = e 2h n Im whereĤ k is the lattice Hamiltonian in the momentum space, E nk is the energy of the n-th band, and u nk is the corresponding Bloch wave function. The integral is performed within the first Brillouin zone and the summation is taken over the occupied bands. In Eq. S9, the magnetization is scaled into the unit of e/2h · meV. Figure. S5 shows the electric field induced magnetization for MnBi 2 Te 4 under exactly the same parameters used in the main text. In general, the results here independently confirm the quantized topological magnetoelectric effect in axion insulators, which, within the computational accuracy, equals to e/2h when the thickness approaches infinite.

VIII. NON-EQUILIBRIUM GREEN'S FUNCTION METHOD FOR THE TIME DEPENDENT CURRENT IN THE AITJ JUNCTION
The full Hamiltonian for the axion insulator tunnel junction shown in Fig. 3 in the main text is whereN u (t) is the particle number in lead u.T uC (=T † Cu ) is the coupling matrix between contact u and the central region. G < Cu (t, t) and G < uC (t, t) are lesser Green's functions. We have used the Heisenberg equation of motion in deriving Eq. S12. Since the Hamiltonian is time dependent and periodic, we can thus perform the double Fourier transform C < uC(Cu) (t, t ) = 1 2π d d e −i t/ e −i t / C < uC(Cu) ( , ), which recasts the current as As will be shown later, the Green's function G < Cu(uC) ( , ) is nonzero only when − = n ω where ω is the magnetic frequency and n is an integer, i.e., G < Cu(uC) ( , ) = n G < Cu(uC) ( , + n ω)δ( − − n ω). Consequently, the current can be further written as Equivalently, For convenience, we define G < uC(Cu),mn ( ) = G < uC(Cu) ( + m ω, + n ω). In this form, we finally obtain Equation. S16 is the central formula to calculate the time dependent current induced by a harmonic magnetic field in an axion insulator, from which we notice that the current I(t) = n I n e inωt where I n = i n + i * −n with i n = where θ(x) is the step function, Ψ i (t) [Ψ † j (t )] is the annihilation (creation) operator for electron in the system, {· · · } denotes the anticommutator and · · · is the state average on the basis of H f ull (t). Perform the equation of motion, we obtain After a double Fourier transform 1 2π dt e −i t / dte i t/ i∂ t ]G r (t, t ), we then have It is obvious that G r ( , ) is nonzero only when − = n ω with n being an integer. Eq. S19 is thus equivalent to If Φ 0 = 0 or ω = 0 (see Eqs. S2 and S10), the system returns to the equilibrium case and the Green's function is only a function of time difference. When Φ 0 = 0, J k (0) = δ k , in this case Eq. S20 can be simplified into ( + m ω)G r mn ( ) = δ mn + { α H α + H T + H 0 + i [T x + T † x ]}G r mn ( ) or G r ( ) = I + H f ull G r ( ). Similarly when ω = 0, k J k (x = 0) = 1. Since the Green's functions now are only a function of time difference, Eq. S20 can be written as G r ( ) = I + { α H α + H T + H 0 + ki [T x J k (yΦ 0 ) + T † x J k (−yΦ 0 )]}G r ( ), which is also G r ( ) = I + H f ull G r ( ). Generally, the retarted Green's function G r mn ( ) can then be solved iteratively or self-consistently [7]. In our numerical calculation, we use the wide band approximation by assuming that the Green's function of the metallic leads does not dependent on the energy.
Alternatively, the retarted Green's function can also be obtained from the Dyson equation. We separate the full Hamiltonian into two part H f ull (t) = H 0 + H 1 (t) where H 0 = α H α + H T + H 0 while H 1 (t) = ki [T x J k (yΦ 0 ) + T † x J k (−yΦ 0 )] exp (ikωt). It is clear that H 0 is time independent and solvable. The free Green's function on the basis of H 0 can then be defined as where · · · 0 refers to the state average on the basis of H 0 instead of the full Hamiltonian. According to the definition, the Dyson equation of the Green's function is G r (t, t ) = g r (t, t ) + dt 1 g r (t, t 1 ) Because H 0 is time independent, its free Green's function g r (t, t ) is only a function of the time difference t − t . Therefore, we obtain the Dyson equation after a Fourier transform G r ( , ) = g r ( )δ( − ) + g r ( ) or G r mn ( ) = g r mm ( )δ(m − n) + g r mm ( )
On the other hand, the free Green's function g r mm ( ) = 1/[ + m ω − H 0 ], Eq. S24 is exactly the same with Eq. S20 after multiplying + m ω − H 0 on both sides.
The lesser Green's function G < mn ( ) can then be obtained using the Keldysh equation where the Fermi distribution function f k ( ) = 1 when ≤ 0 while f k ( ) = 0 otherwise. G a mn ( ) = [G r nm ( )] † is the advanced Green's function. Once we obtain the Green's function G < mn ( ), the current can then be calculated straightforwardly from Eq. S16.