Quantum nonlinear ac transport theory at low frequency

Based on the nonequilibrium Green’s function (NEGF), we develop a quantum nonlinear theory to study time-dependent ac transport properties in the low frequency and nonlinear bias voltage regimes. By expanding NEGF in terms of time to the linear order in Wigner representation, we can explicitly include the time-dependent self-consistent Coulomb interaction induced by external ac bias. Hence this theory automatically satisfies two basic requirements, i.e. current conservation and gauge invariance. Within this theory, the nonlinear ac current can be evaluated at arbitrarily large bias voltages under the low frequency limit. In addition, we obtain the expression of time-dependent current under the wide band limit and derive the relation between the nonlinear electrochemical capacitance and the bias voltage, which are very useful in predicting the dynamical properties of nanoelectronic devices. This quantum theory can be directly combined with density functional theory to investigate time-dependent ac transport from first-principles calculation.


I. INTRODUCTION
Nonlinear Hall effect [1][2][3][4] is the emerging frontier of condensed matter research [5][6][7][8][9][10][11][12][13][14][15][16][17][18] .Prominent nonlinear Hall phenomena include the second order Hall effect discovered in WTe 2 [19][20][21] as well as TaIrTe 4 22 , and the third order Hall effect verified in T d -MoTe 2 23 , etc.These higher-order responses are relatively weak and their measurement requires the phase lock-in amplifier technique, where low frequency ac bias voltage is applied as the driving signal while the second [20][21][22] and third harmonic 23 responses labeling nonlinear signals are probed.However, theoretical interpretation of nonlinear Hall effects is mostly based on the semiclassical Boltzmann approach.In terms of the second 13 and third order 14 conductances expressed in nonequilibrium Green's function (NEGF), quantum transport properties of nonlinear Hall effects have been discussed, but only in the dc case.Therefore, a nonlinear ac transport theory at low frequency is on demand and definitely helps to comprehend nonlinear Hall effects.
On the other hand, time-dependent quantum transport in molecular and nanoscale devices have attracted intensive research attention both experimentally [24][25][26][27][28][29][30][31][32][33][34] and theoretically  . Time-ependent response in quantum transport is of great importance, since it can provide critical information on dynamical properties that are absent in the dc case.The dynamical property is essential in accurate design of nanoelectronic devices.Under time-dependent external fields, quantum transport investigation focuses on two different regimes, i.e., the transient regime when the external bias is immediately turn on or off and the long time ac steady-state regime.For the transient dynamics that characterizes the relaxation time of an electronic system when turned on or off, various approaches have been developed, including scattering matrix theory 37,48,50 , NEGF 38,47,67 , and timedependent Schrödinger equation 46 .In the ac steady-state regime, Büttiker et al derived a current conservation theory which explicitly includes the displacement current in terms of scattering matrix 37 .However, this theory is only applicable under near-equilibrium condition.Later, the situation has been extended to the far-from-equilibrium condition by using NEGF 68 .
According to Ref. [37], the long-range Coulomb interaction must be included in time-dependent and nonlinear dc quantum transport theories, which is necessary for satisfying two basic requirements, i.e., current conservation and gauge invariance 37,69 .The Coulomb interaction can be treated at least in the Hartree level in density functional theory (DFT).Therefore, it is convenient to adopt the first-principles approach, i.e., DFT carried out within the Keldysh NEGF (NEGF-DFT) formalism, to study the transport properties of nanoelectronic devices.In practice, quantitative predictions of their transport properties were compared with experimental results [70][71][72][73] .Although the NEGF-DFT formalism has been widely applied in predicting dc quantum transport properties, its applications to ac situations received far less attention.This is due to the lack of a quantum ac transport theory in the nonlinear bias voltage regime.Such an ac transport theory can be directly coupled with DFT to predict quantum nonlinear ac transport properties from atomic first principles, which is absent so far.Therefore, a nonlinear ac NEGF theory is urgently needed.
In this work, we develop a quantum nonlinear formalism to study ac transport properties at low frequency, where Coulomb interaction is explicitly treated within the NEGF formalism.Based on the adiabatic approximation, we expand the Floquet NEGF with respect to Frozen Green's function in Wigner representation at low frequency.Physically, Frozen Green's function describes that the potential experienced by electrons during their transport is instantaneously adjusted to the applied ac bias.Following this route, the nonlinear time-dependent current is obtained.Furthermore, a concise formula is given under the wide band limit (WBL), which can recover the dynamic conductance at small ac bias 68 .In particular, we demonstrate that the present theory can be used to study the nonlinear electrochemical capacitance in response to external bias, i.e., the C-V curve.
The paper is organized as follows.In Sec.II, we introduce the theoretical model for ac transport.Sec.III presents the derivation of time-dependent NEGF in terms of Frozen Green's function.In Sec.IV, we provide the time-dependent ac current formula and discuss two limiting cases, i.e., wide band limit and small ac bias limit.Nonlinear electrochemical capacitance is also discussed.Finally, a brief summary is given in Sec.V.

II. MODEL
In this section, we introduce the model system for ac transport.The system under investigation is a quantum device in contact with two leads extending to electron reservoirs where time-dependent ac bias is applied, as shown in Fig. 1.The Hamiltonian of this model system can be expressed in three parts Here H α is the Hamiltonian of lead α where C † kα /C kα creates/annihilates an electron in lead α.ǫ kα (t) = ǫ (0) kα + qv α (t), where ǫ (0) kα is the energy level of lead α and v α (t) = v α cos(ωt + ϕ α ) is the time-dependent external ac bias.ω is the frequency of this ac bias and ϕ α denotes the corresponding phase.The second term H C is the Hamiltonian of the isolated central region, where d † n /d n is the creation/annihilation operator of the electron in the central scattering region.Internal Coulomb potential U n under the Hartree approximation is included in the central region, which is defined as where V nm is the matrix element of the Coulomb potential.In real space representation, V (x, The exchange and correlation interactions can be treated in a similar way.The last term H T describes the coupling between the central region and leads with t kαn the coupling constant.Based on NEGF, the time-dependent charge current flowing in lead α can be written as ( = −q = 1) 43 , Here G r,a,< (Σ r,a,< ) is the time-dependent Green's function (self energy) with double time indices in principle.In order to obtain the time-dependent current I α (t), one has to know G r,a,< and Σ r,a,< .In the following section, we will discuss how to calculate them in Wigner representation by expanding them to the first order in frequency, and then evaluate I α (t).

III. GREEN'S FUNCTION AND SELF ENERGY IN WIGNER REPRESENTATION
We start from the time-dependent retarded and lesser Green's functions of the central scattering region on realtime axis 67 , and In general, the Hamiltonian H C (t) in Eqs.(7) where ρ(x, t) is the charge density and x represents the real space position in the central region.Note that Eqs. ( 7), (8), and (9) are coupled, which means that they should be solved self-consistently.The Green's functions in Eqs.(7) and ( 8) depend on two time indices (t, t ′ ), since the time translational invariance is broken due to the presence of the time-dependent ac bias voltage applied in the lead.In the following, we first express the two-time function A(t 1 , t 2 ) in Wigner coordinates, where a slow classical timescale t = (t 1 + t 2 )/2 and a fast quantum timescale τ = t 1 −t 2 are introduced 67 .Thus, A(t 1 , t 2 ) is transformed into A(t, τ ) and an integral expression such as gives 67 where A(t, E) = dτ e iEτ A(t, τ ) and is the gradient operator and arrow indicates the direction of differentiation.When the external ac bias is slow enough compared to characteristic timescales of the electron in the system, the gradient operator in Eq. ( 12) can be expanded order by order in time or energy.Here, we keep only the first order of the Taylor expansion For simplicity, abbreviations for derivatives ∂ ∂E = ∂ E and ∂ ∂t = ∂ t are adopted in the following derivation.Thus Wigner transformation of a two-time function gives 67 where Poisson bracket is defined as Therefore, by taking Wigner transformation and using Eq. ( 14), the corresponding Floquet Green's functions of Eqs. ( 7) and ( 8) up to the first order become (detailed derivation is presented in Appendix A) and Here we have used ) and its Wigner transformation is expressed as In order to calculate the Green's function G r/< (t, E), the self energy Σ r/< (t, E) should be expanded up to the first order in time.We show that the self energy up to the first order is the same as that of the equilibrium one.
We start from the retarded Green's function g r α (t, E) of lead α first, since the self energy Σ r α (t, E) is defined as Σ r α (t, E) = H Cα g r α (t, E)H αC with H Cα describing the coupling between the central region and lead α.Similar to Eq. ( 16), the retarded Green's function of lead α can be expressed as where η is an infinity small positive number.Let's introduce the Frozen Green's function of the lead, Then Eq. ( 19) becomes Since the Poisson bracket is actually calculating the first order term, g r α in it should be replaced by the Frozen Green's function g r α,f directly.And the term where we have used the relation 21) and ( 22), we have Then the self energy Σ r α (t, E) is equal to Furthermore, we calculate the lesser self energy up to the first order in time.As shown in Ref. [68], the lesser Green's function of lead α can be obtained from retarded and advanced Green's functions where g r,a kα is defined as By taking the Wigner transformation, Eq. ( 25) can be written as where 20) and (23).According to the definition, the lesser self energy can be expressed as where ) is the Fermi distribution function of lead α with applied bias.The linewidth function Γ α is defined as It is also equal to Therefore, the self energy in Eqs. ( 16) and ( 17) can be calculated by using Frozen Green's function of the lead.Finally, we introduce the Frozen retarded Green's function of the central scattering region Then Eqs. ( 16) and ( 17) become and Furthermore, we express Floquet Green's function in terms of Frozen Green's function up to the first order in time, and where we have introduced lesser Frozen Green's function After some derivations, we have the following relations and where the relation Notice that the time-dependent Coulomb interaction U (t) is obtained through solving the Poisson equation Eq. ( 9).The current conservation law is expressed as, with total charge of the system Since we wish to expand the time-dependent current I α (t) up to the first order in frequency, the first order of lesser Green's function g < can be safely neglected during calculation of Q(t), i.e., Thus the corresponding Poisson equation can be solved.Therefore, the Coulomb potential U (t) and G < f (t, t) should be solved self-consistently.In the following section, we will show that once the Coulomb potential is self-consistently obtained, current conservation α I α (t) = 0 is automatically guaranteed.

IV. TIME-DEPENDENT CURRENT FORMULA
Similarly, the time-dependent current I α (t) in Eq. ( 6) can also be expressed in Wigner representation where {...} is the Poisson bracket defined in Eq. ( 15).Green's functions and self energies inside the bracket is short for operators in Wigner representation, e.g., G r,a,< ≡ G r,a,< (t, E), etc.Note that the time-dependent Coulomb potential is explicitly included in the frozen Green's function through Eq. ( 31) and hence the timedependent current.Based on Eqs. ( 34) and ( 35), the current can be divided into two parts: α .Here, I (0) α corresponds to the adiabatic current and α is the first order correction due to frequency.With the help of the relation α can be written as where Correspondingly, the first order correction current is Plugging Eqs. ( 36) and (37) into Eq.( 44), the explicit expression of I α can be obtained after some straightforward algebra.Then we have From Eq. ( 45) can be further expressed as where g a ≡ (g r ) † and we have used integration by parts for energy and introduced two terms The physical meaning of [dN α /dE] xx is similar to injectivity [74][75][76] : the local charge at position x that an electron injects from lead α and emits to all leads due to the external bias.
Since current conservation is one of the basic requirements in quantum transport theory, we demonstrate that the current shown in Eq. ( 46) is conserved when the Coulomb interaction is included.Firstly, it is straightforward to show the adiabatic current satisfies α Then from Eq. ( 46), we have α where total charge Q(t) is defined in Eq. ( 40).In the central scattering region, the Poisson equation in Eq. ( 41) should be solved with proper boundary conditions.The natural and reasonable conditions are zero electric fields on the boundaries leading to constant electric potential U there.Thus we have Tr ∇ 2 U (x, t) = 0 due to these boundary conditions of the Poisson equation.Physically, the total charge in the central region Q(t) is zero and so is its derivative, which is known as the charge neutrality condition.Then we arrive at α I α (t) = 0.
A. Wide band limit Wide band limit (WBL) is usually adopted to simplify the analysis of transport properties.In this subsection, we will derive formulas for the time-dependent current under WBL, where self energies are approximately independent of both energy and time, i.e., Therefore, we have and Note that the lesser self energy Σ < is still time-dependent due to the presence of Fermi distribution function f α (E, t) in it.Under WBL, the zeroth order current I (0) α has the same form where Γ α is energy-independent.By using Eqs.( 51), (52), and ( 53), the first order current I α (t) in Eq. ( 46) can be greatly simplified as We can also check the validity of current conservation for the zeroth and first order currents in Eqs. ( 54) and (55).The current conservation of zeroth order can be easily obtained, and the summation of the first order cur-rents in Eq. ( 55) To summarize, the time-dependent current expressed in Eq. ( 55) under WBL is also conserved.

B. Small ac bias limit
At the small ac bias limit, we can derive analytic expressions of the dynamic conductance and timedependent currents, which gives intuitive understanding on quantum ac transport properties.Meanwhile, these analytic results can be compared with the previous work Ref. [68] which focused on quantum ac transport theory at finite frequency and small bias limit.
In general, time-dependent currents in Eqs. ( 54) and ( 55) are in nonlinear response to bias voltage.In the following, we expand the time-dependent current I α (t) order by order in the bias voltage.First, we need to separate the Hamiltonian into two parts: H C (t) = H C0 + U (t), where H C0 is time-independent and U (t) is the time-dependent Coulomb potential.U (t) can be expanded in terms of the bias magnitude v α (t = 0) = v α 68 , with U eq the equilibrium potential in the absence of external bias.U 1 (t) and U 2 (t) are the first and second order corrections due to the presence of external bias, respectively.Correspondingly, u α (t) and u αβ (t) are the first and second order characteristic potentials 37,68,77 .Under Thomas-Fermi approximation, the Poisson-like equations for these characteristic potentials are defined as where Here dn α /dE is the time-dependent injectivity, and The procedure for calculating the first and second order characteristic potentials is shown in Appendix B.
For simplicity, WBL is adopted in the following discussion.The Frozen Green's function in Eq. ( 31) can be expressed as where equilibrium Green's function is defined as ] −1 and self energy Σ r is independent of energy; H ′ C0 = H C0 − U eq and U n is the n-th order correction to the time-dependent Coulomb potential.Now we expand the nonequilibrium Green's function in terms of the first order bias.Thus the Frozen Green's function in Eq. ( 59) can be further written as Correspondingly, the zeroth order current I (0) α (t) in Eq. ( 54) is given by Here we have used the Taylor expansion of Fermi function up to the first order in bias and f 0 is the equilibrium Fermi distribution.In general, by taking Fourier transformation, we can obtain the frequency-dependent current I α (Ω), where Ω is the response frequency.From the definition of dynamic conductance This result agrees with Eq. ( 35) of Ref. [68] when Ω = 0. Furthermore, we can calculate the first order dynamic conductance G (1) αβ (Ω).Since only the first order bias volt-age is considered, Eq. ( 55) is simplified as where Fourier transformation of Eq. ( 63) gives Thus the first order dynamic conductance G (1) which is exactly the same as that in Ref. [68].Notice that the first order dynamic conductance is actually the emittance describing low frequency response of the system 37 .
Here the consistency between Ref. [68] and the present work at the small ac bias limit is demonstrated.It is worth mentioning that, in Ref. [68], the Coulomb potential is treated in a perturbative way, which is suitable for describing ac transport at finite frequency and small ac bias limit.If one tries to simulate the finite ac bias case with the theory developed in Ref. [68], the Coulomb potential has to be expanded in terms of the characteristic potential order by order, which is very complicated especially for nonlinear terms.In contrast, the present work can deal with the Coulomb potential directly by solving the Poisson equation Eq. (41).Compared with Ref. [68], the present work has advantages in investigating nonlinear ac transport with finite ac bias voltages at the low frequency limit.Ref. [68] and the present work focus on different regimes of quantum ac transport and complimentary to each other.The application of our quantum nonlinear ac theory to nonlinear electrochemical capacitance is shown below.

C. C-V curve and nonlinear electrochemical capacitance
In this subsection, we discuss the nonlinear emittance in a system where dc transport is not allowed, for example, the magnetic tunneling junction (MTJ).Then the nonlinear emittance is equivalent to the electrochemical capacitance [78][79][80][81] .
One the other hand, we could assume that the scattering region of a two-probe system can be roughly divided into two regions, i.e., Ω I and Ω II , and expand total charge of the system in terms of bias voltage where α = I, II.Since the bias-dependent Q α can be calculated though expanding Frozen lesser's Green's function G < f in terms of v β (t), we have where Eq. ( 59) is used.Finally, the linear capacitance C αβ (t) and nonlinear capacitance C αβγ (t) are equal to and Here WBL is not assumed in deriving Eq. ( 69).After obtaining the characteristic potentials u β (t) and u βγ (t), one can calculate the voltage-dependent linear and nonlinear electrochemical capacitances.First principles calculation of the time-dependent current or electrochemical capacitance proceeds as follows.First, at a given time t, the lead Hamiltonians and selfenergies are obtained from conventional DFT.Second, one needs to solve the Poisson equation Eq. ( 41) selfconsistently within NEGF-DFT.Once the self-consistent accuracy is reached, we obtain the nonequilibrium Hamiltonian H C (t) and the retarded and lesser Green's functions defined in Eqs.(34) and (35).Finally, we can calculate the time-dependent current from Eqs. ( 43) and (44)  as well as the linear and nonlinear electrochemical capacitances from Eqs. ( 68) and ( 69) to investigate quantum ac transport properties.

V. SUMMARY
In summary, we have developed a quantum nonlinear ac transport theory based on NEGF.At low frequency, through expanding the time-dependent NEGF in terms of Frozen Green's function, we can explicitly include the time-dependent self-consistent Coulomb interaction in the time-dependent current, which is nonlinear in bias voltage.Therefore, this nonlinear response theory automatically satisfies both current conservation and gauge invariance conditions.More importantly, it can be directly combined with DFT to calculate ac transport properties from atomic first principles.As a demonstration, we discuss how to calculate the ac current and nonlinear electrochemical capacitance versus the bias voltage.This quantum nonlinear ac transport theory can be easily extended to multiterminal systems for describing nonlinear Hall responses at low frequency.Furthermore, inelastic process, such as the electron-phonon/photon interaction, can also be included in the present theoretical formalism as effective self-energies to discuss their effects in quantum ac transport.
Then by taking Fourier transform ( dτ e iEτ ) of each term in Eq. ( 7) and using Eqs.( 14) and ( 18 Thus, we can easily arrive at Eq. ( 16) and obtain the lesser Green's function defined in Eq. (17).

APPENDIX B
In this appendix, we demonstrate how to calculate the characteristic potentials defined in Eq. (58).We first solve the first order potential.From Eq. ( 58), we have where the injectivity is defined as 37,68 dn α (t, x) dE = − cos(ωt) with dn(t, x)/dE the local density of states at time t.
In general, the Poisson equation Eq. ( 72) needs to be solved within the NEGF-DFT framework and u α is selfconsistently obtained.For the simplest case, we can adopt the quasineutrality approximation 37 , which means that the charge density at each point is zero so that −∇ 2 u α (t, x) = 0. Then analytic expression of the first order characteristic potential is shown as When the first order potential u α is obtained, we can calculate the second order injectivity dñ αβ /dE: where the second order derivatives with respect to energy can be numerically facilitated with the finite difference method.Then substituting it into Eq.( 58), we have At the small ac bias limit, the second order potential u αβ contributes to ac transport properties less than the first order one, and hence the quasineutrality approximation is appropriate here.†) These authors contributed equally to this work.* jianwang@hku.hk

FIG. 1 .
FIG. 1. Schematic plot of a two-terminal nanoelectronic device consisting of a central scattering region (red line box) and two leads.The left and right leads extend to electron reservoirs at infinity, where time-dependent bias voltages vαcos(ωt + ϕα) are applied and electric currents are probed.
∇ 2 u αβ (t, x) = 4π dñ αβ (t, x) dE − 4π dn dE u αβ (t, x),(77)which also needs to be solved.Notice that here dn/dE = α dn α /dE is the total density of states.Again we can use the quasineutrality approximation, which leads to