Introduction

Tunneling phenomenon is one of the most remarkable and unique consequences of the wave nature of particles in quantum mechanics, where a particle can penetrate through classically forbidden regions. In solids, the quantum mechanical wavefunctions of electrons form the band structure separated by the energy gaps, and the tunneling can occur between these bands when an electric field is applied. This is called Zener tunneling through the energy gap and has been actively studied1,2,3,4,5,6,7,8,9,10,11,12,13. A concise formula, i.e., Landau–Zener formula1,2, has been obtained for a model Hamiltonian describing the two-band system as

$$H=\left(\begin{array}{ll}v\hslash k&\delta \\ \delta &-v\hslash k\end{array}\right),$$
(1)

where ±vℏk are the energy dispersions and 2δ is the energy gap. Under an external electric field E, the wavenumber k is accelerated as \(\hslash \dot{k}=-eE\) as shown in Fig. 1. The transition probability from the lower band to the upper band reads

$$P=\exp \left(-\frac{\pi {\delta }^{2}}{e\hslash Ev}\right),$$
(2)

which is essentially singular with respect to E showing the nonperturbative nature of the quantum tunneling.

Fig. 1: Schematics of Landau–Zener tunneling.
figure 1

A wave packet (characterized by the wavenumber k and band dispersion ϵ(k)) driven by an electric field can tunnel into the conduction band with transition probability P.

Full size image

At a pn-junction of semiconductors, the tunneling shows an asymmetric behavior, which is utilized as a tunneling diode for rectifying devices14. Because of the broken inversion symmetry, the tunneling probability, and hence, the I–V characteristics depend strongly on the direction of the electric field E. For the uniform bulk crystal, however, the asymmetry in the Zener tunneling probability is a highly nontrivial issue even when the crystal lacks the inversion symmetry. This can be seen in the band dispersion εn(k) (n: band index); the relation εn(k) = εn(−k) holds due to the time-reversal symmetry even in the absence of the spatial inversion symmetry. Therefore, the inversion symmetry is rather hidden in wave mechanics15. Intuitively, the extended wave state is rather insensitive to the broken inversion symmetry compared with the localized wave packet. Therefore, a fundamental question is how the nonreciprocal behavior, i.e., the asymmetry between the opposite direction of the electric field E, is realized in the tunneling processes of the bulk crystals, reflecting the wave nature of the electrons. This is also an important issue in terms of device applications; Ferroelectric random access memory utilizes nonreciprocal current response at the time of read-out operations of recorded polarization direction16, while its working mechanism has not been fully understood so far.

The nonreciprocal phenomena in noncentrosymmetric crystals have been extensively studied in these days, including both the dc transport17,18,19,20,21,22 and photo-excited current23,24,25,26,27,28,29,30,31,32,33. In particular, the no-go theorem has been proposed for the nonreciprocal transport of independent particles induced by the static electric field, in terms of a perturbative expansion with respect to E34. Thus nonreciprocal dc transport requires some sort of interaction effects in the perturbative regime. On the other hand, this theorem does not apply for the photocurrent induced by the light irradiation that induces the inter-band transitions, which is called shift current. The shift current is formulated in terms of the Berry connection of the Bloch wavefunctions, which correspond to the intracell coordinates of the electrons30,31,32,33,35. The optical transition causes the shift in the intracell coordinates, i.e., shift vector, since intracell coordinates are generally different for the valence and conduction bands in noncentrosymmetric crystals. The steady pumping of polarization of photoexcited electron–hole pairs results in the dc photocurrent. Therefore, it is concluded that the wavefunctions encode the information of the noncentrosymmetry in sharp contrast to the energy dispersion. In fact, the Berry phase becomes zero (or trivial) when the system preserves both the inversion and time-reversal symmetries.

As discussed above, the tunneling is a nonperturbative effect, and cannot be captured by the perturbative expansion with respect to E. Hence, it is possible that the nonreciprocal nature appears in the Landau–Zener tunneling even in the independent particle approximation. In this paper, we show that this is indeed the case by deriving the generalized Landau–Zener formula including the shift vector, i.e., the information of the Bloch wavefunctions. The generalized Landau–Zener formula shows that nonreciprocal tunneling generally appears in inversion broken systems, even in the presence of the time-reversal symmetry. The nonreciprocity has a geometric origin, dominated by the Berry connection difference between the valence and conduction bands. We also give a demonstration of nonreciporcal tunneling by applying the obtained formula to the Rice–Mele model, an archetypal one-dimensional model of a ferroelectrics.

Results

Tunneling formula with a shift vector

Let us consider a time evolution of a system under a slow change of parameters. In particular, here we focus on a change of momentum k under a DC electric field, k → k(t) = k − eEt/ℏ. It is well known that the solution of the time-dependent Schrödinger equation in the adiabatic limit is given by snapshot eigenstates

$$H(t)\left|n,k(t)\right\rangle ={\varepsilon }_{n}(t)\left|n,k(t)\right\rangle$$
(3)

multiplied by dynamical and Berry phase factors [see Eq. (4)]. The diabatic correction is derived from the transition dipole matrix elements. To see this, let us expand a state vector \(\left|\Psi \right\rangle\) by the adiabatic solutions as

$$\left|\Psi \right\rangle =\sum _{n}{a}_{n}(t){e}^{-i\int\nolimits_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{n}({t}_{1})+eE{A}_{nn}({t}_{1})]/\hslash }\left|n,k(t)\right\rangle ,$$
(4)

where \({A}_{\mathrm{nm}}(t)=i\left\langle n,k(t)\right|{\partial }_{k}\left|m,k(t)\right\rangle\) is the Berry connection. (We note that the “off-diagonal” Berry connections for n ≠ m correspond to transition dipole matrix elements.) With paying attention in dealing with the Berry phase factor, we can reduce the time-dependent Schrödinger equation \(i\hslash {\partial }_{t}\left|\Psi \right\rangle =H(t)\left|\Psi \right\rangle\) to

$$i{\partial }_{t}{a}_{n}(t)=\frac{eE}{\hslash }\sum _{m\ne n}| {A}_{nm}(t)| {e}^{i\int_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{n}-{\varepsilon }_{m}+{{{eER}}_{nm}}]/\hslash +i\arg {A}_{nm}({t}_{0})}{a}_{m}(t),$$
(5)

with

$${R}_{nm}={A}_{nn}-{A}_{mm}-{\partial }_{k}\arg {A}_{nm}.$$
(6)

See “Methods” for details. Here we have used ℏ∂t = −eE∂k for \(\left|n,k(t)\right\rangle\) and \(\arg {A}_{nm}(t)\). Rnm is so called shift vector which is a gauge-invariant object describing the polarization difference between two bands n, m30,31,32,33,35,36. Specifically, the shift vector can be interpreted as the spatial shift of the wave packets between the valence and conduction bands, since the Berry connections Ann correspond to the intracell coordinates of the Bloch electrons. It is known that the shift vector plays an important role in bulk photogalvanic effect in noncentrosymmetric crystals, so called shift current. The fact that Rnm appears in Eq. (5) is usually overlooked since Ann = 0 is assumed in many cases. A similar geometric contribution has also been discussed as the quantum geometric potential37,38, in the context of the adiabatic condition.

Let us focus on a tunneling process between two bands, n = ± with a−(t0) = 1, a+(t0) = 0. Our goal is to derive the tunneling probability P = ∣a+(t)∣2 after one cycle of the Bloch oscillation. For simplicity, we consider only the first-order correction w.r.t. ∣A+−∣ here. By integrating Eq. (5) and using it recursively, we obtain

$${a}_{+}(t)= \ i{e}^{i\arg {A}_{+-}({t}_{0})}\\ \times \int\nolimits_{{k}_{0}}^{k(t)}d{k}_{1}\left|{A}_{+-}\right|\exp \left[-i\mathop{\int}\nolimits_{{k}_{0}}^{{k}_{1}}d{k}_{2}\left(\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{eE}+{R}_{+-}\right)\right]$$
(7)

with k0 = k(t0), as we detail in “Methods”.

A two-band Hamiltonian can be represented as H =  d(k) ⋅ σ with σ being Pauli matrices (when we subtract a constant energy shift). The quantities necessary for the evaluation of the tunneling amplitude are given as

$${\varepsilon }_{+}-{\varepsilon }_{-}=2\sqrt{{{\boldsymbol{d}}}^{2}},$$
(8)
$$| {A}_{+-}| =\frac{\sqrt{{({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})}^{2}}}{2{{\boldsymbol{d}}}^{2}},$$
(9)
$${R}_{+-}=\frac{({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot ({\partial }_{k}^{2}{\boldsymbol{d}})}{{({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})}^{2}}\sqrt{{{\boldsymbol{d}}}^{2}}.$$
(10)

In order to evaluate the integral in an asymptotic manner, we employ the Dykhne–Davis–Pechukas (DDP) method3,4 in accordance with Ref. 3. Namely, we perform the integral by means of contour integration in the complex plane. The contour of the integral, which is originally the real axis [blue lines in Fig. 2], can be deformed within an analytic region, thanks to the Cauchy’s integral theorem4,39.

Fig. 2: Complex k plane that governs the tunneling process.
figure 2

The branching point kc is a point where the band gap vanishes in the complex plane of momentum k. Cauchy's theorem allows to deform the original integration path on the real axis (blue) into the contour in the complex plane that passes through the branching point (magenta). The absolute value of the exponential factor is constant along the black contours. The direction of the original path is different for a positive electric field E > 0 and b negative electric field E < 0, and the choice of the branching point is also different.

Full size image

This treatment is advantageous since one can utilize a (complex) branching point kc, where the energy gap vanishes \({\boldsymbol{d}}{({k}_{c})}^{2}=0\) [Such a point is indeed a branching point when the Hamiltonian is analytic, as \({\varepsilon }_{+}-{\varepsilon }_{-}\propto {(k-{k}_{c})}^{1/2}\) in the vicinity of kc]. This point essentially governs the tunneling process between the two bands: since the prefactor ∣A+−∣ diverges as we approach kc [see Eq. (9)], only this divergent part contributes to the asymptotic value of the integral, when the integration path is deformed to pass through the vicinity of the branching point kc.

We show the integration path by magenta lines in Fig. 2. The main part of the contour is one along which the absolute value of the exponential factor is constant (i.e., the imaginary part of the k2 integral in Eq. (7) is constant). This contour passes through the branching point kc, but we make a detour around it since kc itself is a singular point of the integrand. Due to the divergence mentioned above, this detoured part contributes dominantly against the main part. While the branching points appear in a pairwise manner \(({k}_{c},{k}_{c}^{* })\), we need to choose one of them such that the exponential factor becomes smaller than unity in accordance with the detailed derivation given in “Methods”. We note that the integral on the first and last vertical lines have finite but small contributions in general. They share the same absolute value (that is perturbative in E) but their phases are different. This leads to a small oscillation in the tunneling amplitude with respect to E, on top of the nonperturbative contribution from the branching point kc. Since we are interested in the nonperturbative contribution, we neglect those perturbative corrections in the following.

In a generic situation, one can assume that the leading order term of d2, \({({\partial }_{k}{\boldsymbol{d}})}^{2}\) and \(({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot ({\partial }_{k}^{2}{\boldsymbol{d}})\) in the expansion around kc is given as d2 ~ iα(k − kc), \({({\partial }_{k}{\boldsymbol{d}})}^{2} \sim \beta\), and \(({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot ({\partial }_{k}^{2}{\boldsymbol{d}}) \sim \eta\), respectively. By evaluating the detoured part of the integral [circular arc around kc in Fig. 2] with these expanded forms, we arrive at

$$P \sim \exp \left[2\,{\text{Im}}\,\int\nolimits_{{k}_{0}}^{{k}_{c}}{{d{k}}}_{2}\left(\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{{eE}}+{R}_{+-}\right)\right],$$
(11)

as we describe in “Methods”. We note that there actually is a prefactor (π/3)2 ~ 1.1, which appears because we have evaluated the probability within the first order of the adiabatic perturbation. However, we dropped this because it is known3 that the prefactor becomes exactly 1 if we take into account all order contributions in the adiabatic perturbation (for details, see “Methods”). We also note that the integral in Eq. (11) is invariant against replacing the lower limit k0 to an arbitrary k on the real axis (usually it is set to the gap minimum).

The obtained formula, Eq. (11), includes the geometric correction described by the shift vector. This contribution is absent in the original DDP formula, because they assumed that the 2 × 2 Hamiltonian is real; under this assumption d is a two-dimensional vector, where \(({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot ({\partial }_{k}^{2}{\boldsymbol{d}})=0\) always holds. On the other hand, Refs. 6,7,8 dealt with a generic 2 × 2 Hamiltonian, so that the geometric correction indeed appeared in their study. However, their calculation was done in a particular gauge, and the obtained expression is not written with gauge-invariant quantities. Comparison of the formulae is given in “Methods”. Our result is expressed in terms of gauge-invariant quantities, the shift vector, so that it provides much clearer understanding on the physical interpretation of the correction and remarkable physical consequences due to it.

The shift vector Rnm plays a crucial role in the nonlinear transport of inversion-broken systems. It satisfies Rnm(k) = Rnm(−k) when the system is time-reversal symmetric, while Rnm(k) = −Rnm(−k) when the system is inversion symmetric. Thus when the system has both symmetries, there is no correction to the tunneling probability; this is consistent with the fact that one can make the Hamiltonian real in such cases. A nontrivial result thus can appear when either symmetry is broken. In particular, a qualitatively new phenomenon appears when the inversion symmetry is broken: when we change the sign of the electric field E, we need to change the choice of the branching point from kc to  −kc in order to have a correct result P ≤ 1 with the formula (11) (see Fig. 2a, b). Under this alternation, the shift vector contribution is invariant in the time-reversal broken system, which leads to a simple correction of the probability independent of the field strength/direction. On the other hand, when the inversion symmetry is broken, the exponent of the shift vector correction is odd under this alternation, so that it leads to an exponentially-large difference in the tunneling probability when the direction of the electric field is reverted.

The physical meaning of the shift vector as an intracell coordinate provides an intuitive understanding of its role in the tunneling process, as follows. The tunneling process can be interpreted as a propagation of a wave packet thorough a classically forbidden region in real space, which has a thickness of 2δ/e|E| as drawn in Fig. 3. Now the noncentrosymmetric system has an internal degree of freedom, the intracell coordinate, whose difference between two bands is represented by the shift vector R as shown in Fig. 3a. Thus the tunneling from the lower band to the upper band induces a spontaneous shift of the wave packet by R, which effectively modifies the tunneling depth as 2δ/e|E| + R. Since the direction of the shift vector is intrinsically determined by the underlying crystal, this modification of the tunneling depth by R contributes in a constructive/destructive manner in accordance with the direction of the bias, as shown in Fig. 3b,c. As can be seen in the tunneling formula (11), this modification of the tunneling depth is essentially governed by the shift vector R+− (defined with the Berry connection difference between the two bands) at the band gap minimum.

Fig. 3: Zener tunneling in inversion broken systems.
figure 3

a Real space picture in the absence of the electric field E. The energy spectrum ϵ with an energy gap 2δ is depicted as a function of position x. The positions of the wave packets in the two bands are shifted in the unit cell by the shift vector R. b, c Real space pictures of two bands in the presence of E. The shift in the intracell coordinate results in the different tunneling depth depending on the direction of E by the shift vector R. Namely, the shift of the wave packet by R upon the interband transition effectively modifies the tunneling depth as 2δ/e|E| + R.

Full size image

Application to Rice–Mele model

Let us see the emergence of the nonreciprocal (direction-dependent) tunneling probability using the Rice–Mele model (with the lattice constant a = 1)

$$H=\left(\begin{array}{ll}m&t\cos k-i\delta t\sin k\\ t\cos k+i\delta t\sin k&-m\end{array}\right),$$
(12)

a prototypical example of an inversion-broken (polar) system. Indeed, this model has a finite shift vector

$${R}_{+-}=\frac{mt\delta t\sqrt{{t}^{2}{\cos }^{2}k+\delta {t}^{2}{\sin }^{2}k+{m}^{2}}}{{m}^{2}(\delta {t}^{2}{\cos }^{2}k+{t}^{2}{\sin }^{2}k)+{t}^{2}\delta {t}^{2}}.$$
(13)

Since the integrand of Eq. (11) for the present system is invariant under k → k ± π, let us consider a slow parameter change k : 0 → −πsgn(E). The branching point of the present model is given as

$${k}_{c}^{(\!\pm\!)}=\frac{2n+1}{2}\pi \pm i{\tanh }^{-1}\sqrt{\frac{\delta {t}^{2}+{m}^{2}}{{t}^{2}+{m}^{2}}}$$
(14)

with \(n\in {\mathbb{Z}}\). By evaluating the generalized formula (11) with \({k}_{c}={k}_{c}^{(-{\text{sgn}}(E))}\), we obtain

$$P= \ \exp \left[-\frac{4\sqrt{{t}^{2}+{m}^{2}}}{e| E| }\left(K(\gamma )-E(\gamma )\right)\right.\\ -\left. \frac{2m\,{\text{sgn}}(E)}{t\delta t\sqrt{{t}^{2}+{m}^{2}}}\left({t}^{2}K(\gamma )-({t}^{2}-\delta {t}^{2})\Pi \left(\frac{\delta {t}^{2}}{{t}^{2}},\gamma \right)\right)\right]$$
(15)
$$\sim \exp \left[-\frac{\pi (\delta {t}^{2}+{m}^{2})}{e| E| t}-\,{\text{sgn}}\,(E)\frac{\pi m\delta t}{2{t}^{2}}\right]\quad (\delta t,m\ll t),$$
(16)

where K(γ), E(γ), and Π(n, γ) are the complete elliptic integrals with \(\gamma =\sqrt{(\delta {t}^{2}+{m}^{2})/({t}^{2}+{m}^{2})}\) being the elliptic modulus. In the last line, we recover the conventional exponent Eq. (2) for the ∣E∣−1 term (\(\delta \to \sqrt{\delta {t}^{2}+{m}^{2}}\) and ℏv → t). The shift vector correction is represented as −sgn(E) × (π/2)∣Imkc∣ × R+− evaluated at k = π/2. Remarkably, as we have mentioned, the non-zero shift vector not just provides an exponentially-large correction, but also a strong nonreciprocity via sgn(E).

Using the obtained generalized Landau–Zener formula, we show the tunneling probability P in Rice–Mele model in Fig. 4. Plots of P(E) as a function of the applied electric field E (Fig. 4a) show a well-known nonpertrubative behavior at E = 0 and a good agreement with tunneling probabilities obtained by numerically solving time-dependent Schrödinger equation (shown in dots) (for details, see “Methods”). We note that the oscillation with a small amplitude seen in the numerical results arises from the perturbative corrections due to the integral along the first and last vertical lines of Fig. 2. In addition, the plots show that tunneling probabilitites differ depending on the direction of E (the sign of E). This direction dependence arises from the nonzero shift vector and shows almost proportionality to the value of R+− at the band gap minimum (on the real axis of k). Figure 4b shows the nonreciprocity ratio P(−E)/P(+E) (the ratio of tunneling probabilities for negative and positive E), which quantifies the strength of nonreciprocity as a rectifying device, as a function of the strength of alternating hopping δt that introduces inversion symmetry breaking. We find that the nonreciprocity ratio P(−E)/P(+E) grows monotonically, as we increase δt and the effect of inversion symmetry breaking gets stronger. For small δt, the nonreciprocity ratio is linearly proportional to δt. In particular, large nonreciprocity of P(−E)/P(+E) ~ 2 can be achieved for a feasible value of δt ~ 0.5t, which indicates that Landau–Zener tunneling in noncentrosymmetric crystals is able to realize strong nonreciprocal functionality.

Fig. 4: Nonreciprocal tunneling in Rice–Mele model.
figure 4

a Nonreciprocal tunneling probability P(E) in Rice–Mele model, plotted as a function of electric field E. Solid lines represent P(E) from the generalized Landau–Zener formula, and dots represent P(E) obtained from numerical simulations of time-dependent Schrödinger equation. We used the parameters (δt, m) = (0.2t, 0.4t). b Nonreciprocity ratio P(−E)/P(+E) plotted as a function of alternating hopping amplitude δt that controls the magnitude of inversion breaking. We set the staggered potential m = 0.4t. t is the (uniform) hopping amplitude.

Full size image

Nonreciprocal tunneling in a continuous model

As we have mentioned, the small-amplitude oscillation of the tunneling probability seen in Fig. 4a is derived from the first and last vertical segments of the deformed contour shown in Fig. 2. This contribution survives due to the finite bandwidth, and should vanish when the energy gap is infinite at the initial and final states. To verify this, we consider a Hamiltonian in a continuous space

$$H=\left(\begin{array}{ll}m\sqrt{1+{k}^{2}}&-tk-i\delta t\\ -tk+i\delta t&-m\sqrt{1+{k}^{2}}\end{array}\right),$$
(17)

whose low-energy structure coincides with that of the Rice–Mele model (12), while the energy spectrum is unbounded. The generalized formula (11) for this Hamiltonian reads

$$P= \,\exp \left[\frac{-\pi (\delta {t}^{2}+{m}^{2})}{e| E| \sqrt{{t}^{2}+{m}^{2}}}\right. \\ \left. -\frac{2m\,{\text{sgn}}(E)}{t\delta t\sqrt{{t}^{2}+{m}^{2}}}\left({t}^{2}K(\gamma )-({t}^{2}-\delta {t}^{2})\Pi \left(\frac{\delta {t}^{2}}{{t}^{2}},\gamma \right)\right)\right].$$
(18)

It is worth noting that the nonreciprocal part (∝ sgn(E)) is identical to that in Eq. (15), while the reciprocal part is modified from Eq. (15). We compare the formula with numerically computed tunneling amplitudes and nonreciprocity ratios in Fig. 5. The numerical result no longer shows oscillation with respect to E from the perturbative correction and shows better agreement with the generalized Landau–Zener formula.

Fig. 5: Nonreciprocal tunneling in the continuous model.
figure 5

a Nonreciprocal tunneling probability P(E), plotted as a function of electric field E. Solid lines represent P(E) from the generalized Landau–Zener formula, and dots represent P(E) obtained from numerical simulations of time-dependent Schrödinger equation. We used the parameters (δt, m) = (0.2t, 0.4t). b Nonreciprocity ratio P(−E)/P(+E) plotted as a function of alternating hopping amplitude δt that controls the magnitude of inversion breaking. We set the staggered potential m = 0.4t. t is the (uniform) hopping amplitude.

Full size image

Discussion

Here we further consider the role of time-reversal symmetry T in the nonreciprocal responses. In the context of magnetochiral anisotropy17,18,19,20,21, it has been discussed that the nonreciprocal transport requires the broken T in addition to the broken inversion symmetry. This can be understood intuitively that the reversal of time corresponds to that of the current direction when there is no dissipation. The no-go theorem by Morimoto and Nagaosa34 indeed shows that the current proportional to the square of the electric field is forbidden in non-interaction systems with dc electric field, when the T-symmetry is preserved. However, it has been revealed that this is not the case when the inter-band transitions and the associated absorption of energy occur due to ac E-field, where the shift current is induced. The Zener tunneling studied in the present paper can be regarded as the “inter-band” transition under dc E-field due to the tunneling, and hence the no-go theorem, which is based on the adiabatic assumption that there occurs no inter-band transition, does not apply. This is related to the non-analytic and non-perturbative nature of the tunneling probability, which cannot be expanded in E. In a more general situation in the plane of the frequency ω and the strength E of the electric field, there are several different regions as indicated by Fig. 17 in Aoki et al.40. In this respect, the shift current and Zener tunneling are two limiting cases, i.e., ac limit of weak E-field and dc limit of strong E-field, and the crossover between these two corresponds to the Keldysh line and how the nonreciprocal responses behaves in this plane is an interesting problem to be studied in the future. Note that from the viewpoint of the symmetry, the requirement for the nonreciprocal Zener tunneling is the same as that for the shift current. Another unique feature of the nonreciprocal Zener tunneling is that the ratio of the tunneling currents for the two directions is independent of the electric field E and is of the order of unity as indicated in Eq. (16).

It is an interesting future problem to investigate further geometric corrections. There might be new effects that do not emerge until the higher-order perturbation is taken into account41,42, which can lead to even larger nonreciprocity. Since the geometric correction discussed in the present study appears in a combination of ERmn, and we have evaluated the nonperturbative contribution within the leading order of E (for the prefactor), such new effects might be found in the presence of strong E field.

In the present paper, we considered the spinless model. Incorporating the spin degrees of freedom and spin-orbit interaction, one can expect the nonreciprocal charge and spin tunneling currents. These phenomena can offer novel mechanisms for spin diode or switchable diode by a magnetic field. The large spin polarization of the electrons transmitted through the DNA molecules that has been observed experimentally43 might be understood from this point of view44.

Methods

Adiabatic perturbation theory

Equations (5) and (7) in the main text are derived using the adiabatic perturbation theory as follows. By inserting Eq. (4) to the time-dependent Schrödinger equation, we obtain

$$\sum _{m}{e}^{-i\int\nolimits_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{m}({t}_{1})+eE{A}_{mm}({t}_{1})]/\hslash }\left\{[i\hslash {\partial }_{t}{a}_{m}(t)+eE{a}_{m}(t){A}_{mm}(t)]\left|m,k(t)\right\rangle -eE{a}_{m}(t)i{\partial }_{k}\left|m,k(t)\right\rangle \right\}=0.$$
(19)

Here we have used \(\hslash {\partial }_{t}\left|m,k(t)\right\rangle =-eE{\partial }_{k}\left|m,k(t)\right\rangle .\) By taking an inner product with \(\left\langle n,k(t)\right|\), 〈n, k(t)∣m, k(t)〉 = δnm and \(\left\langle n,k(t)\right|i{\partial }_{k}\left|m,k(t)\right\rangle ={A}_{nm}(t)\) leads to

$$i{\partial }_{t}{a}_{n}(t)=\frac{eE}{\hslash }\sum _{m\ne n}{A}_{nm}(t){e}^{i\int_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{n}({t}_{1})+{{eEA}_{nn}}({t}_{1})]/\hslash }{e}^{-i\int\nolimits_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{m}({t}_{1})+eE{A}_{mm}({t}_{1})]/\hslash }{a}_{m}(t)$$
(20)
$$=\frac{eE}{\hslash }\sum _{m\ne n}| {A}_{nm}(t)| {e}^{i\int\nolimits_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{n}-{\varepsilon }_{m}+eE({A}_{nn}-{A}_{mm})]/\hslash +i\arg {A}_{nm}(t)}{a}_{m}(t)$$
(21)
$$=\frac{eE}{\hslash }\sum _{m\ne n}| {A}_{nm}(t)| {e}^{i\int\nolimits_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{n}-{\varepsilon }_{m}+eE({A}_{nn}-{A}_{mm})+\hslash {\partial }_{t}\arg {A}_{nm}]/\hslash +i\arg {A}_{nm}({t}_{0})}{a}_{m}(t)$$
(22)
$$=\frac{eE}{\hslash }\sum _{m\ne n}| {A}_{nm}(t)| {e}^{i\int\nolimits_{{t}_{0}}^{t}d{t}_{1}[{\varepsilon }_{n}-{\varepsilon }_{m}+eE({A}_{nn}-{A}_{mm}-{\partial }_{k}\arg {A}_{nm})]/\hslash +i\arg {A}_{nm}({t}_{0})}{a}_{m}(t).$$
(23)

This coincides with Eq. (5).

By integrating Eq. (5) from t0 to t, we obtain a perturbative expansion in terms of the adiabatic perturbation theory,

$${a}_{n}(t) = {a}_{n}({t}_{0})-i\frac{eE}{\hslash }\int\nolimits_{{t}_{0}}^{t}d{t}_{1}\sum \limits_{m\ne n}| {A}_{nm}({t}_{1})| {e}^{i\mathop{\int}\nolimits_{{t}_{0}}^{{t}_{1}}d{t}_{2}[{\varepsilon }_{n}-{\varepsilon }_{m}+eE{R}_{nm}]/\hslash +i\arg {A}_{nm}({t}_{0})}{a}_{m}({t}_{1})$$
(24)
$$ ={a}_{n}({t}_{0})-i\frac{eE}{\hslash }\int\nolimits_{{t}_{0}}^{t}d{t}_{1}\sum \limits_{m\ne n}| {A}_{nm}({t}_{1})| {e}^{i\mathop{\int}\nolimits_{{t}_{0}}^{{t}_{1}}d{t}_{2}[{\varepsilon }_{n}-{\varepsilon }_{m}+eE{R}_{nm}]/\hslash +i\arg {A}_{nm}({t}_{0})}{a}_{m}({t}_{0})\\ -{\left(\frac{eE}{\hslash }\right)}^{2}\mathop{\int}\nolimits_{{t}_{0}}^{t}d{t}_{1}\\ \int\nolimits_{{t}_{0}}^{{t}_{1}}d{t}_{2}\sum _{m\ne n}\sum _{l\ne m}| {A}_{nm}({t}_{1}) {A}_{ml}({t}_{2})| {e}^{i\int\nolimits_{{t}_{0}}^{{t}_{1}}d{t}_{3}[{\varepsilon }_{n}-{\varepsilon }_{m}+eE{R}_{nm}]/\hslash +i\int\nolimits_{{t}_{0}}^{{t}_{2}}d{t}_{3}[{\varepsilon }_{m}-{\varepsilon }_{l}+eE{R}_{ml}]/\hslash +i\arg ({A}_{nm}({t}_{0}){A}_{ml}({t}_{0}))}{a}_{l}({t}_{2})\left.\right].$$
(25)

We then neglect the third term O(∣A∣2) in the last line, and set n = + , m = − . We obtain Eq. (7) as

$${a}_{+}(t)=-i\frac{eE}{\hslash }\int\nolimits_{{t}_{0}}^{t}d{t}_{1}| {A}_{+-}({t}_{1})| {e}^{i\int\nolimits_{{t}_{0}}^{{t}_{1}}d{t}_{2}[{\varepsilon }_{+}-{\varepsilon }_{-}+eE{R}_{+-}]/\hslash +i\arg {A}_{+-}({t}_{0})}$$
(26)
$$=i\mathop{\int}\nolimits_{k({t}_{0})}^{k(t)}d{k}_{1}| {A}_{+-}({k}_{1})| {e}^{-i\mathop{\int}\nolimits_{k({t}_{0})}^{k({t}_{1})}d{k}_{2}[\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{eE}+{R}_{+-}]+i\arg {A}_{+-}({t}_{0})},$$
(27)

since a−(t0) = 1 and a+(t0) = 0.

Dykhne–Davis–Pechukas (DDP) method

We evaluate Eq. (7) as Eq. (11) by using the DDP method with the deformed contour shown in Fig. 2, as we sketch in the main text. The detail of the evaluation is as follows. We focus on the detoured part (the circular arc around kc), which yields a dominant contribution. As mentioned in the main text, we assume that d2 ~ iα(k − kc), \({({\partial }_{k}{\boldsymbol{d}})}^{2} \sim \beta\), and \(({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot ({\partial }_{k}^{2}{\boldsymbol{d}}) \sim \eta\) in the vicinity of kc. Then \({({\partial }_{k}{\boldsymbol{d}}\,\times {\boldsymbol{d}})}^{2}={({\partial }_{k}{\boldsymbol{d}})}^{2}{{\boldsymbol{d}}}^{2}-{({\partial }_{k}{{\boldsymbol{d}}}^{2})}^{2}/4 \sim {\alpha }^{2}/4\) leads to

$${\tilde{A}}_{+-} \sim \pm \frac{{\text{sgn}}\,({\text{Im}}\,{k}_{c})}{4i(k-{k}_{c})},$$
(28)

where, to avoid confusion, we have introduced \({\tilde{A}}_{+-}=\sqrt{{({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})}^{2}}/(2{{\boldsymbol{d}}}^{2})\) [see Eq. (9)] as an analytic function on the complex plane that coincides with  \(| {A}_{+-}|\) on the real axis. This, remarkably, does not depend on any detail (α, β, η, etc.) of the system. The sign ± should be chosen such that \({\tilde{A}}_{+-}\ge 0\) holds on the real axis (typically  +). We also expand the exponent as

$$\int\nolimits_{{k}_{0}}^{{k}_{1}}d{k}_{2}\left(\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{eE}+{R}_{+-}\right) \sim z+\int\nolimits_{{k}_{0}}^{{k}_{c}}d{k}_{2}\left(\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{eE}+{R}_{+-}\right),$$
(29)

where

$$z=\frac{4}{3}\sqrt{i\alpha }\left(\frac{1}{eE}-\frac{2\eta }{{\alpha }^{2}}\right){({k}_{1}-{k}_{c})}^{3/2}.$$
(30)

The main part of the contour (where the absolute value of the exponential factor is constant) extends from kc to directions where z is real. There are three such directions in 2π/3 intervals, as can be seen in Fig. 2. The circular arc connects two of them, so it has an angle 2π/3.

We choose either kc or \({k}_{c}^{* }\), such that \(\arg z\in [0,\pi ]\) holds on the arc [Note: In typical systems (e.g., the Rice–Mele model), α > 0 for Imkc > 0 and α < 0 for Imkc < 0 hold. In such a case we should choose \({k}_{c}^{(-{\text{sgn}}(E))}\). This makes the contour depend on the direction of E, as shown in Fig. 2a, b]. Then further we set the infinitesimal radius of the arc to be ∝ ∣E∣2/3−ϵ with ϵ > 0. By doing so we can make the contribution of detoured part dominate that of the main part (for rigorous bounds, see Ref. 3. The shift vector correction 2η/α2 does not affect the bounds as it is higher-order w.r.t. E).

We then change the integrating variable from k1 to z. The integral (7) along the detoured contour becomes

$${a}_{+}(t) \sim \pm {e}^{i\arg {A}_{+-}({t}_{0})}\exp \left[-i\int\nolimits_{{k}_{0}}^{{k}_{c}}d{k}_{2}\left(\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{eE}+{R}_{+-}\right)\right]\int\nolimits_{C}dz\frac{{e}^{-iz}}{6z},$$
(31)

where the contour C is a semicircle with a radius ∝ E−ϵ → ∞ as E → 0, covering the upper half-plane [note: the sign ± is identical for kc and \({k}_{c}^{* }\) (typically  −). This sign is derived from \({\tilde{A}}_{+-}\) and the direction of the contour]. With the help of Jordan’s lemma, the integral along C is given by the residue at z = 0, i.e., ∮ dze−iz/(6z) = iπ/3. Finally we arrive at Eq. (11) multiplied by (π/3)2.

Actually it is known that the prefactor π/3 should be corrected by taking into account all order contributions in the adiabatic perturbation3. To obtain the correct tunneling amplitude, we need to multiply the z integrand of Eq. (31) by a−(z). a−(z) with the higher-order correction can be obtained as follows: (i) differentiating Eq. (5) with respect to k leads to the differential equation for a−,

$${\partial }_{k}^{2}{a}_{-}-\left[{\partial }_{k}{\mathrm{ln}}\,| {A}_{+-}| +i\left(\frac{{\varepsilon }_{+}-{\varepsilon }_{-}}{{eE}}+{R}_{+-}\right)\right]{\partial }_{k}{a}_{-}+| {A}_{+-}{| }^{2}{a}_{-}=0.$$
(32)

(ii) In the vicinity of the branching point kc, the differential equation is simplified as

$$36{z}^{2}{\partial }_{z}^{2}{a}_{-}+36z(1-iz){\partial }_{z}{a}_{-}-{a}_{-}=0,$$
(33)

which has an asymptotic series solution,

$${a}_{-}(z)=1-\mathop{\sum }\limits_{n=1}^{N}\frac{n!{i}^{n}}{36{n}^{2}{z}^{n}}{\mathop{\prod }\limits_{k = 1}^{n-1}}\left(1-\frac{1}{36{k}^{2}}\right).$$
(34)

With this solution, the z integral is evaluated as

$$\oint dz\frac{{e}^{-iz}}{6z}{a}_{-}(z)= \frac{i\pi }{3}\left[1-\mathop{\sum }\limits_{n=1}^{N}\frac{1}{36{n}^{2}}{\mathop{\prod }\limits_{k = 1}^{n-1}}\left(1-\frac{1}{36{k}^{2}}\right)\right]$$
(35)
$$=\frac{i\pi }{3}{\mathop{\prod }\limits_{k = 1}^{N}}\left(1-\frac{1}{36{k}^{2}}\right)\mathop{\longrightarrow }\limits^{N\to \infty }2i\sin \frac{\pi }{6}=i$$
(36)

instead of iπ/3. Hence we end up with Eq. (11) with the prefactor exactly 1 instead of (π/3)2.

Comparison of Eq. (11) with the method in Joye et al.7

The formula for the tunneling amplitude in Joye et al.7 can be reproduced by evaluating Eq. (11) with eigenvectors

$$\left|\pm ,k(t)\right\rangle =\frac{(\sqrt{{{\boldsymbol{d}}}^{2}}\pm {\boldsymbol{d}}\cdot {{\boldsymbol{e}}}_{z})\left|\uparrow \right\rangle \pm {\boldsymbol{d}}\cdot ({{\boldsymbol{e}}}_{x}+i{{\boldsymbol{e}}}_{y})\left|\downarrow \right\rangle }{\sqrt{2({{\boldsymbol{d}}}^{2}\pm \sqrt{{{\boldsymbol{d}}}^{2}}{\boldsymbol{d}}\cdot {{\boldsymbol{e}}}_{z})}}.$$
(37)

For this gauge choice, we obtain

$${A}_{++}-{A}_{--}=\frac{(({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot {{\boldsymbol{e}}}_{z})({\boldsymbol{d}}\cdot {{\boldsymbol{e}}}_{z})}{\sqrt{{{\boldsymbol{d}}}^{2}}{({\boldsymbol{d}}\,\times {{\boldsymbol{e}}}_{z})}^{2}},$$
(38)
$$\arg {A}_{+-}=\frac{\pi }{2}-{\tan }^{-1}\frac{2\sqrt{{{\boldsymbol{d}}}^{2}}({\boldsymbol{d}}\,\times {\partial }_{k}{\boldsymbol{d}})\cdot {{\boldsymbol{e}}}_{z}}{({\boldsymbol{d}}\cdot {{\boldsymbol{e}}}_{z}){\partial }_{k}{{\boldsymbol{d}}}^{2}-2({\partial }_{k}{\boldsymbol{d}}\cdot {{\boldsymbol{e}}}_{z}){{\boldsymbol{d}}}^{2}}.$$
(39)

If we assume that d(kc) ⋅ ez ≠ 0 (in accordance with Ref. 7), \(\arg {A}_{+-}({k}_{c})=\pi /2\) holds in general cases, so that \(\,{\text{Im}}\,\int\nolimits_{{k}_{0}}^{{k}_{c}}dk{\partial }_{k}\arg {A}_{+-}=0\). Hence, we can replace R+− in Eq. (11) by A++ − A−−, which leads to the expression given in Refs. 6,7. While this expression is obtained in a particular gauge, our formula (11) is gauge-invariant and thus free from the assumption d(kc) ⋅ ez ≠ 0. Note that the integrand (38) has no physical meaning while the integrated value does. In particular, it diverges as \(\sim {(k-{k}_{c})}^{-1/2}\) as k → kc, and hence not suitable for numerical evaluation. Let us compare the two formulae, \(\,{\text{Im}}\,\int\nolimits_{\pi /2}^{k}d{k}_{2}{R}_{+-}({k}_{2})\) and \(\,{\text{Im}}\,\int\nolimits_{\pi /2}^{k}d{k}_{2}({A}_{++}({k}_{2})-{A}_{--}({k}_{2}))\) with Eq. (38), with a concrete example. We plot these integrals for the Rice–Mele model Eq. (12) with (δt, m) = (0.4t, 0.4t) in Fig. 6. While the integrals terminated at a generic k do not coincide, they do at the branching point kc.

Fig. 6: Comparison of the formulae.
figure 6

Comparison of the integral \(\,{\text{Im}}\,\int\nolimits_{\pi /2}^{k}d{k}_{2}{R}_{+-}({k}_{2})\) with Eq. (10) (blue) and \(\,{\text{Im}}\,\int\nolimits_{\pi /2}^{k}d{k}_{2}({A}_{++}({k}_{2})-{A}_{--}({k}_{2}))\) with Eq. (38) (magenta) for the Rice–Mele model. A±± is the Berry connection of the energy band ±, and R+− is the shift vector between these bands. The momentum k here is on the segment [π/2, kc] on the complex plane, where kc is the branching point.

Full size image

Numerical simulations

The solutions of time dependent Shrödinger equations were obtained by numerically solving differential equations with Mathematica NDSolve function that uses an adaptive procedure to determine the size of integration steps.