Manipulation of magnetization and spin transport in hydrogenated graphene with THz pulses

Terahertz (THz) field pulses can now be applied in scanning tunneling microscopy (THz-STM) junction experiments to study time-resolved dynamics. The relatively slow pulse compared to the typical electronic time-scale calls for approximations based on a time-scale separation. Here, we contrast three methods based on non-equilibrium Green’s functions: i) the steady-state, adiabatic results, ii) the lowest-order dynamic expansion in the time variation, and iii) the auxiliary mode propagation method without approximations in the time variation. We consider a concrete THz-STM junction setup involving a hydrogen adsorbate on graphene where the localized spin polarization can be manipulated on/off by a local field from the tip electrode and/or a back-gate affecting the in-plane transport. We use steady-state non-equilibrium Green’s function theory combined with density functional theory to obtain a Hubbard model for the study of the junction dynamics. Solving the Hubbard model in a mean-field approximation, we find that the near-adiabatic first-order dynamic expansion in the time variation provides a good description for STM voltage pulses up to the 1 V range.

1 Introduction

Recently, it has become possible to extend scanning tunneling microscopy (STM) studies to the time domain by applying strong sub-cycle near-field electromagnetic pulses in the THz regime inside the STM, so-called THz-STM.. This enables studies of time-resolved dynamics of the transport in junctions with sub-Å and sub-ps resolution [1–3]. Due to the strong field confinement at the STM tip electrode, the effective voltages between tip and sample can be on the order of 1 V [4].

The STM-THz experiments pose interesting questions regarding the theoretical description of the dynamics of the electronic subsystem on the time scale of typical atomic vibrations. The case of atomic-scale junctions in STM under strong bias and coupling to the tip requires the consideration of coupling to semi-infinite electrodes. To this end, the non-equilibrium Green’s function (NEGF) methods [5–7] have been a popular choice, where the electrodes are included via self-energies. However, it is relevant to consider simplifications to the full time dynamics since the electron dynamics most often nearly adiabatically follow fields in the THz range. This can be accomplished by the Wigner representation, which involves a time-scale separation into a slow THz time scale (T) and a fast electronic time scale (τ), typically around a femtosecond. This has been considered recently by Honeychurch and Kosov [8] where the NEGF Kadanoff–Baym equations (KBEs) are expressed in the so-called Wigner form with an expansion in terms of the derivative of the slow time, T, from hereon called “dynamic expansion” (DE).

From a numerical point of view, one important advantage of the DE approach is that each discrete time step will be independent from each other: quantities like electron density and current are expressed by quantities independent of the other time steps allowing for computational parallelization over time as opposed to the time propagation of the KBEs. This is important for efficient calculations based, e.g., on first-principles methods such as time-dependent density functional theory based on the NEGF (TD-DFT-NEGF).

In this paper, we apply the DE method to a realistic STM junction which has received considerable attention, namely, STM on a hydrogen adsorbate on graphene. This system is magnetic, but the magnetism can be tuned by the applied tip or gate field, resulting in a change in the in-plane spin transport. Here, we will consider the junction subject to a THz pulse and use it to benchmark the DE against the full time-dependent NEGF calculation with the auxiliary mode (AM) method and steady-state field approximation (zero-order DE). We first use steady-state DFT-NEGF to examine the bias- and gate-dependent electronic structure and in-plane transport in the system. We extract parameters from here to obtain an effective mean field, a Hubbard model for the junction, which we use for benchmarking and which serves as a toy model of a full TD-DFT-NEGF calculation, highlighting the numerical considerations in a self-consistent treatment of time-dependence within the DE approach. Our results show how the computationally efficient DE method is able to capture the main features in the THz dynamics of the system, in particular that the THz-bias pulse gives rise to fast switching dynamics of the magnetism, generating higher harmonic dynamics in the occupations and currents which, on the other hand, are not captured well by the steady-state adiabatic approximation.

2 System and DFT calculations

Magnetism of graphene can be created by point defects like adsorption of isolated hydrogen (H/Gr) or a vacancy, owing to the creation of unpaired π electrons around the defects [9]. The electronic density of graphene is highly tunable by electrostatic gating [10] which makes it interesting to consider how this magnetism may be tunable either by a global back-gate or a local gating such as the voltage on the STM tip, see Figures 1A,B.

FIGURE 1

FIGURE 1. (A) Schematic of the three-terminal hydrogen on the graphene (H/Gr) device setup with a charge back-gate plane. μ_tip and μ_Gr are chemical potentials of the gold tip and graphene electrodes. Here, W is the width of the device. (B) Local atomic structure of the hydrogen atom on graphene showing an out-of-plane buckling of the C atom below H (sp³ hybridization). (C) Density of states projected on the H atom, the midgap states, with g = 0 and g = −1 gating in units of 10¹³cm⁻², respectively. Solid red and dashed blue lines represent spin up and down populations, respectively. E_F refers to the Fermi energy of the system in equilibrium, i.e., without bias.

First, we consider steady-state DFT-NEGF calculations where we use a three-terminal device setup with a gate plane, encompassing left and right graphene electrodes (with width W), and an additional gold tip electrode above the H atom, as schematically shown in Figure 1A. The calculations were performed using the SIESTA/TranSIESTA package [11–14] with the GGA-PBE [15] functional for exchange correlation, a DZP atomic orbital basis set and an electronic temperature of 50 K (for further details, see the work of Gao et al. [16]). To model the gate-induced doping of graphene, a gate plane was placed 15 Å underneath the graphene. The gate carries a charge density of n = g × 10¹³ e/cm², where g defines the gating level, with g < 0 (g > 0) corresponding to n(p) doping [17].

The local atomic structure of H/Gr in Figure 1B displays an out-of-plane buckling of the C atom, leading to a transition from sp² to sp³ hybridization. In Figure 1C, we show the calculated density of states around the Fermi energy (E_F) projected on the H atom (PDOS) without gating, g = 0, and with gate-induced n-doping, g = −1. The H impurity resonances at the center of the graphene’s pseudogap, the midgap peaks, imply a strong interaction between the s orbital of the H atom and the p_z orbital of graphene. Here, the H adatom forms a σ bond with the carbon atom next to it, and π bonds are broken. The bonding states are located at −8 eV, far below the Fermi level.

We use a tip-H distance of 4.5 Å, where there is only a weak connection between the tip and H to keep the spin moment at 1 μ_B as the initial state at zero bias (0 V) and gate. Without gating, the magnetic moment can be switched on/off with the tip voltage-induced doping below the tip, as observed by the disappearance of the spin splitting in the PDOS. On the other hand, we can also turn the magnetism “off” for an n-doped graphene by applying gating (g = −1), where we see a non-spin split, fully occupied peak for zero tip bias (0 V). In this situation the magnetism/spin splitting reappears when we apply a 0.8 V tip bias counteracting the gate doping. Thus, it is possible to manipulate the spin by either the STM tip bias or by a global field from a gate plane [16]. The main driver for this behavior is the electronic occupation of the carbon below the hydrogen.

Next, we address the effects of the field and spin manipulation on the in-plane transport in H/Gr. To this end, we consider the energy-dependent scattering cross section, σ(E) [18], corresponding to the “shadowing” length (area) caused by the point defect in 2D (3D) conductors. We can estimate σ(E) from the electron transmission functions for the pristine and defected systems, T₀ and T, respectively [19, 20].

$\sigma \left(E\right)=W\frac{T_0\left(E\right)-T\left(E\right)}{T\left(E\right)},(1)$

where W is the width of the device, as shown in Figure 1A. In Figure 2A, we see two peaks around E_F caused by the resonant, spin-split H resonances, c.f. the PDOS in Figure 1C, yielding a significant in-plane transport spin polarization, P_σ = (σ_↑ − σ_↓)/(σ_↑ + σ_↓). The cross section we obtain here is close to the one obtained for H on graphene nanoribbons [18]. In accordance with the PDOS in Figure 1C, we observe in Figure 2A how the spin splitting and transport spin polarization vanish when applying a tip voltage. On the other hand, as shown in Figure 2B, we show how the in-plane transport spin polarization can be turned “on” using the counteracting tip bias for the gate-induced n-doped system. In this case, we also observe a dip in the cross section around −0.3 eV below E_F, corresponding to the Dirac point of the doped graphene substrate. The H resonance is pinned around the Fermi level at equilibrium (0 V), but it can be shifted by the tip bias. This shows how external biases may be used to turn spin polarization “on/off” in the H/Gr system and how this has implications not only for the tunnel current between the tip and sample (PDOS) but also for the in-plane transport in graphene.

FIGURE 2

FIGURE 2. In-plane spin-dependent elastic scattering cross section (σ) as a function of energy. In the bottom panel, the corresponding spin polarization (P_σ) without (A) g = 0 and (B) with g =−1 gate-induced n doping is shown. Solid red and dashed blue lines represent spin-up and spin-down states, respectively.

3 Model and time-dependent methods

We consider next the response of the H/Gr system to a time-dependent variation of the tip and gate potentials with the waveform shown in Figure 3A chosen to mimic THz pulses in typical experiments [1, 4], $V(t)=A\cos (2\pi {\nu }_{0}t-\varphi )\exp \left(-\frac{t^2}{2{\zeta }^2}\right)$. The time-resolved shape and power spectrum of the pulse are shown in Figure 3. We note that a realistic experimental pulse integrates to zero over time, which is not fulfilled by the simple analytic expression specified above. However, we do not expect this deviation to have a significant effect on the dynamics. The central (carrier) frequency ν₀ can be adjusted to investigate how the dynamics depend on the speed of external driving, thereby testing the adiabaticity in the system. The width is adjusted along with the frequency such that the dimensionless product ζν₀ = 0.3 remains constant.

FIGURE 3

FIGURE 3. THz pulse used to drive the system in the time-domain (A) and frequency domain (B). The pulse has the form $V(t)=A\cos (2\pi {\nu }_{0}t-\varphi )\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{t^2}{2{\zeta }^2}\right)$. For the numerical investigations carried out in this work, we fixed the parameters at $\zeta =\frac{3}{10}{\nu }_0^{-1}$, $\varphi =-\frac{\pi }{4}$, and A = 1 eV. These parameters were chosen to mimic pulses used in the experiment. The carrier frequency ν₀ is varied between 1 and 10 THz; in the figure, a pulse with ν₀ = 3 THz is shown.

We will use a minimal model for the H/Gr system with parameters chosen according to the aforementioned DFT calculations and model the system in the typical electrode–device–electrode setup where the time-dependence of the potential is considered in the device region, while the electrodes are assumed to be perfectly screened so that the time dependence in these regions consists only of a rigid shift of the energies and chemical potentials following the pulse [5]. The device region consists of the hydrogen atom and the sp³ carbon atom in the graphene on which it is adsorbed. This system is modeled using an electronic two-level model for each spin, corresponding to one orbital per atom per spin. The coupling to the graphene sheet and the STM tip electrodes is included using self-energies. We note that the spin-polarization in the H/Gr in reality involves carbon atoms further away from the adsorbate [9], but we will neglect these effects here.

The system is described by the total Hamiltonian

$\widehat{\mathcal{H}}\left(t\right)=\widehat{H}\left(t\right)+{\widehat{H}}_{\text{gr}}\left(t\right)+{\widehat{H}}_{\text{tip}}\left(t\right)+{\widehat{V}}_{d,\mathrm{g}\mathrm{r}}+{\widehat{V}}_{d,\mathrm{t}\mathrm{i}\mathrm{p}},(2)$

where ${\widehat{H}}_{\text{gr}}$ describes the graphene sheet, ${\widehat{H}}_{\text{tip}}$ the STM tip, $\widehat{H}$ the device region, and ${\widehat{V}}_{d,\mathrm{g}\mathrm{r}}$ and ${\widehat{V}}_{d,\mathrm{t}\mathrm{i}\mathrm{p}}$ the graphene-device and tip-device coupling, respectively. The rigid shift of the energies in the screened electrodes caused by the electric field in the THz pulse is Δ_αV(t), where α ∈ {gr, tip}, and we have

${\widehat{H}}_{\alpha }\left(t\right)={\displaystyle \sum _{k,\sigma }}\left[{\epsilon }_{\alpha k}+{\mathrm{\Delta }}_{\alpha }V\left(t\right)\right]{\widehat{c}}_{\alpha k\sigma }^{†}{\widehat{c}}_{\alpha k\sigma }.(3)$

Here, αk is a generalized state label referring to the eigenstates of the isolated graphene sheet and tip, respectively. The mean-field Hubbard Hamiltonian $\widehat{H}$ describing the device is

$\widehat{H}\left(t\right)={\displaystyle \sum _{\sigma }}\left[{\epsilon }_H+V_{\text{DC}}+U{n}_{H\bar{\sigma }}\left(t\right)\right]{\widehat{c}}_{H\sigma }^{†}{\widehat{c}}_{H\sigma }+{\epsilon }_C\left(t\right){\widehat{c}}_{C\sigma }^{†}{\widehat{c}}_{C\sigma }+\left(v_{CH}{\widehat{c}}_{C\sigma }^{†}{\widehat{c}}_{H\sigma }+\mathrm{h}.\mathrm{c}.\right),(4)$

where ${\widehat{c}}_{C,H}^{(†)}$ are annihilation (creation) operators for the hydrogen and adsorption-site orbitals, ɛ_H and ɛ_C are the on-site energies for these orbitals, v_CH is the hopping matrix element between these, and σ ∈ {↑, ↓} refers to the spin. The Hubbard interaction strength is given by U, and the opposite-spin H occupation is given by $n_{H\bar{\sigma }}=⟨{\widehat{c}}_{H\bar{\sigma }}^{†}{\widehat{c}}_{H\bar{\sigma }}⟩$, where $\bar{\sigma }$ refers to the spin opposite to σ. The parameter V_DC designates a DC potential due to the tip bias, as shown in Figures 1, 2. The device couples to the tip by tunneling through the hydrogen orbital

${\widehat{V}}_{\text{d,tip}}={\displaystyle \sum _{k\sigma }}{v}_{\text{tip},k}{\widehat{c}}_{H\sigma }^{†}{\widehat{c}}_{\text{tip},k\sigma }+\mathrm{h}.\mathrm{c}.,(5)$

and to the graphene sheet through the C atom on the adsorption site

${\widehat{V}}_{\text{d,gr}}={\displaystyle \sum _{k\sigma }}{v}_{\text{gr},k}{\widehat{c}}_{C\sigma }^{†}{\widehat{c}}_{\text{gr},k\sigma }+\mathrm{h}.\mathrm{c}.,(6)$

such that v_tip,k is the hopping matrix element between the tip state k and the hydrogen orbital and v_gr,k is the hopping matrix element between the graphene state k and the device carbon orbital. These parameters enter the numerical model only indirectly through the electrode self-energies to be specified below.

The values of the model parameters were estimated to mimic the results of the DFT calculation, and we find that the values ɛ_H = −1.7 eV, ɛ_C(t) = 0 eV + Δ_grV(t), v_CH = 3.25 eV, U = 6.5 eV, Δ_gr = 1, and Δ_tip = 0 resulted in a qualitative agreement between the models when comparing PDOS on the respective atoms and their spin splitting. The effect of the THz pulse is to shift the graphene electronic energies relative to the tip. There is no spin-flipping mechanism such that the spin-up and spin-down electrons effectively form two separate systems, interacting only through the time-dependent mean-field occupation of the opposite spin on the H atom, which must be solved self-consistently.

The dynamical evolution of the system for a given time-dependence may be described using the NEGF formalism. We now express all operators as matrices in the device state space, consisting of the hydrogen and carbon orbital for each spin. The retarded, advanced, and lesser Green’s functions in the device region are governed by the KBEs [5],

$\begin{array}{c}\left[i{\partial }_t-H\left(t\right)\right]G^{r/a}\left(t,t^{\prime }\right)=\mathbf{1}\delta \left(t-t^{\prime }\right)+\int \mathrm{d}{t}_1{\mathrm{\Sigma }}^{r/a}\left(t,t_1\right)G^{r/a}\left(t_1,t^{\prime }\right),\end{array}(7)$

for the retarded/advanced functions and

$\begin{array}{c}\left[i{\partial }_t-H\left(t\right)\right]G^{<}\left(t,t^{\prime }\right)=\int \mathrm{d}{t}_1\left[{\mathrm{\Sigma }}^r\left(t,t_1\right)G^{<}\left(t_1,t^{\prime }\right)+{\mathrm{\Sigma }}^{<}\left(t,t_1\right)G^a\left(t_1,t^{\prime }\right)\right]\end{array}(8)$

for lesser Green’s function. In these equations, Σ = Σ_gr + Σ_tip is the self-energy related to the coupling of the device to the graphene sheet and the tip. These self-energies may be found from Green’s function operator g_α of the electrodes in the absence of coupling to the device:

${\mathrm{\Sigma }}_{\alpha }^{r/a/<}\left(t,t^{\prime }\right)={\displaystyle \sum _k}{v}_{\alpha ,k}{g}_{\alpha ,k}^{r/a/<}\left(t,t^{\prime }\right)v_{\alpha ,k}^{*}.(9)$

In the present work, the tip is modeled in the wide-band limit using an energy constant, ${\mathrm{\Sigma }}_{\text{tip}}^{r/a}=\mp i{\mathrm{\Gamma }}_{\text{tip}}/2$, with the weak tip-hydrogen coupling Γ_tip = 0.1 eV, while the self-energy of the graphene sheet is found from a tight-binding model of the graphene with the nearest-neighbor hopping parameter t = −2.7 eV and zero on-site energy. From this model, the self-energy accounting for the presence of the rest of the graphene sheet on the carbon atom under the hydrogen was found using a k-integral to obtain the real-space self-energy in the primitive unit cell of graphene [21]. Subsequently, this self-energy is further downfolded to the C atom entering the C–H system. This yields a self-energy which reproduces the Dirac-cone DOS for graphene.

We note that all quantities in Eqs (7) and (8) are matrices in the device state-space. Thus, in our model, we obtain two sets of 2 × 2 matrices, one set for each spin. To solve the KBEs, we use two different numerical approaches. The first is a method based on the separation of time scales [8], where the resulting equations of motion are expanded in an asymptotic series, which we denote as the DE. This represents an approximation to the dynamics, which becomes exact when the external time-dependence varies on a much slower time scale than the electronic degrees of freedom in the system. The second method (AM) recasts the KBEs into a set of coupled ordinary differential equations (ODEs) which depend only on a single time variable. The system dynamics can then be found by propagating these ODE’s forward in time from a specified initial condition [22, 23]. This method provides an exact propagation of the dynamics, but requires that the self-energies be approximated by a series of Lorentzian functions. In the following, we will provide some more details on these two schemes.

3.1 Dynamical expansion

We introduce a central time variable T and a difference time variable τ,

$T=\frac{t+t^{\prime }}{2}\tau =t-t^{\prime }.(10)$

In the absence of external driving, the system is invariant with respect to translations in time, and thus, in this case, the quantities in the KBEs will only be a function of the difference time τ. The introduction of an external time-dependent driving breaks time-translation invariance and, thus, introduces a dependence on the central time variable T. Therefore, the variables (T, τ) are well suited to separate the time scales associated with internal dynamics of the system (τ) and the external driving (T). After the introduction of the new time variables, we cast the KBEs in the so-called Wigner space representation by performing a Fourier transform over τ. Thus, we obtain a new set of Green’s functions and self-energies:

$\tilde{G}\left(T,\omega \right)={\int }_{-\infty }^{\infty }\mathrm{d}\tau G\left(T,\tau \right)e^{i\omega \tau },(11)$

where the superscript ∼ indicates that a quantity is in the Wigner representation. After transforming all self-energies and Green’s functions to the Wigner representation, the KBEs of motion take the form

$\begin{array}{c}\left[\frac{i}{2}{\partial }_T+\omega -e^{-\frac{i}{2}{\partial }_{\omega }^G{\partial }_T^H}H\left(T\right)\right]{\tilde{G}}^{r/a}=1+e^{-\frac{i}{2}\left({\partial }_T^{\mathrm{\Sigma }}{\partial }_{\omega }^G-{\partial }_{\omega }^{\mathrm{\Sigma }}{\partial }_T^G\right)}{\tilde{\mathrm{\Sigma }}}^{r/a}{\tilde{G}}^{r/a},\end{array}(12)$

and

$\begin{array}{c}\left[\frac{i}{2}{\partial }_T+\omega -e^{-\frac{i}{2}{\partial }_{\omega }^G{\partial }_T^H}H\left(T\right)\right]{\tilde{G}}^{<}=e^{-\frac{i}{2}\left({\partial }_T^{\mathrm{\Sigma }}{\partial }_{\omega }^G-{\partial }_{\omega }^{\mathrm{\Sigma }}{\partial }_T^G\right)}\left({\tilde{\mathrm{\Sigma }}}^{<}{\tilde{G}}^a+{\tilde{\mathrm{\Sigma }}}^r{\tilde{G}}^{<}\right),\end{array}(13)$

where the differential operators in the exponentials act only on the function whose superscript they bear. For notational convenience, we have suppressed the arguments (T, ω) on the self-energies and Green’s functions. We note that Eqs (12) and (13) are still formally exact. If the time-dependence is taken to be adiabatic, all central-time derivatives in Eqs (12) and (13) vanish. In this limit, the exponentials act as the identity operator, and we obtain the well-known steady-state KBEs in the frequency domain:

$\left[\omega -H\left(T\right)\right]{\tilde{G}}_0^{r/a}=1+{\tilde{\mathrm{\Sigma }}}^{r/a}{\tilde{G}}_0^{r/a}.(14)$

$\left[\omega -H\left(T\right)\right]{\tilde{G}}_0^{<}={\tilde{\mathrm{\Sigma }}}^{<}{\tilde{G}}_0^a+{\tilde{\mathrm{\Sigma }}}^r{\tilde{G}}_0^{<},(15)$

where the 0 subscripts indicate that these are adiabatic Green’s functions. If we assume that the external perturbation acts on a much slower time scale than any internal dynamics in the system, G and Σ will be slowly varying functions of the central time T. In that case, we can expand the exponential operators in Eqs (12) and (13) and keep only the lowest order terms on the grounds that higher-order terms contain higher powers of the central time derivative ∂_T, which must be small by virtue of the slow variation. This allows us to expand the Green’s functions in a formal series.

$G=G_0+G_1+G_2+\cdots ,(16)$

where G_i satisfies an equation of motion containing terms only of ith order in the central time derivative ∂_T. Thus, G₀ describes the adiabatic response of the system (Eqs (14)–(15)) to the time-dependent perturbation, while G₁ describes the first-order dynamical corrections and so on. The explicit expressions for G^r/a and G^< resulting from this expansion at the zeroth and first orders of approximation are shown in Supplementary Material.

Once the Green’s functions governing the system have been found, the time-resolved current into the device from the graphene sheet and tip, respectively, is given by [5]

$J_{\alpha }\left(t\right)=\mathrm{Tr}{C}_{\alpha }\left(t\right),(17)$

where the current matrix C_α is

$C_{\alpha }\left(t\right)={\int }_{-\infty }^t\mathrm{d}{t}^{\prime }\left[G_{\alpha }^{>}\left(t,t^{\prime }\right){\mathrm{\Sigma }}_{\alpha }^{<}\left(t^{\prime },t\right)-G_{\alpha }^{<}\left(t,t^{\prime }\right){\mathrm{\Sigma }}_{\alpha }^{>}\left(t^{\prime },t\right)\right]+\mathrm{h}.\mathrm{c}.(18)$

In the Wigner space representation, the integral in Eq. (18) becomes an exponential operator, and we find the equivalent expression in terms of Wigner space quantities as

$\begin{array}{c}{C}_{\alpha }\left(T\right)={\int }_{-\infty }^{\infty }\frac{\mathrm{d}\omega }{2\pi }{e}^{-\frac{i}{2}\left({\partial }_T^G{\partial }_{\omega }^{\mathrm{\Sigma }}-{\partial }_{\omega }^G{\partial }_T^{\mathrm{\Sigma }}\right)}\left[{\tilde{G}}^{>}\left(T,\omega \right){\tilde{\mathrm{\Sigma }}}_{\alpha }^{<}\left(T,\omega \right)-{\tilde{G}}^{<}\left(T,\omega \right){\tilde{\mathrm{\Sigma }}}_{\alpha }^{>}\left(T,\omega \right)\right]+\mathrm{h}.\mathrm{c}.\end{array}(19)$

We expand the current and current matrix in a series in the same manner as we did for the Green’s functions.

$J_{\alpha }=J_{\alpha ,0}+J_{\alpha ,1}+J_{\alpha ,2}+\cdots C_{\alpha }=C_{\alpha ,0}+C_{\alpha ,1}+C_{\alpha ,2}+\cdots ,(20)$

where J_α,i = Tr C_α,i, and C_α,i satisfies an equation containing terms of order i in ∂_T. These equations, at each order, are found from the series expansion of the exponential operator in Eq. (19) and the expansion Eq. (16) of the Green’s functions by collecting all terms of the same order in ∂_T. They are given explicitly in Supplementary Material.

It is worth emphasizing that the method does not require any propagation of the equations of motion because the dynamics are included as time-dependent corrections to the adiabatic solution. It is, therefore, well suited to problems where exact propagation may not be computationally feasible, such as problems that require the consideration of dynamics over very long time scales, and, in particular, the computations may be parallelized trivially over the discrete time steps of the central time variable, T.

3.1.1 Self-consistent solution

The device Hamiltonian is an input parameter to the dynamical expansion and, thus, must be known a priori. This represents a problem for working with systems with interactions or electron–electron repulsion, as the effective device Hamiltonian will then depend on the state of the system. One approach to this problem is the introduction of correlation self-energies, which must then also be approximated by an additional expansion in ∂_T [24]. In this work, we use instead an iterative procedure, seeking a self-consistent dynamical electron density. This approach is in the spirit of DFT, and we envision that future implementations of the method may work with DFT programs to obtain self-consistent dynamical corrections to time-dependent problems. In our approach, we start from a self-consistent adiabatic solution on a pre-defined time-grid. The self-consistent adiabatic Hamiltonian obtained from this solution is then used as the input Hamiltonian to the DE model, from which a new dynamical Hamiltonian is obtained on the same time grid. This process is repeated until the self-consistency of the dynamical Hamiltonian is achieved. Obtaining convergence for all time steps is not trivial except at very low frequencies ν₀ of the external perturbation and was not obtained at the highest frequencies investigated. This is related to the fact that the higher-order harmonics generated by the switching of the magnetization state are highly non-adiabatic and, therefore, unsuitable for description by an expansion such as the DE. Nonetheless, for frequencies below ν₀ ≈ 5 THz, the dynamical corrections were in good agreement with the exact time propagation (AM) scheme even when the self-consistent cycle did not formally converge for all time steps. Further details of the numerical implementation may be found in Supplementary Material.

3.2 Exact density matrix propagation (AM)

The KBE, as shown in Eq. 8, can be dealt with under the assumption that the spin-independent electrode broadening matrices ${\mathrm{\Gamma }}_{\alpha }=i({\mathrm{\Sigma }}_{\alpha }^r-{\mathrm{\Sigma }}_{\alpha }^a)$ are a weighted sum of Lorentzians in the form [22, 23, 25, 26]

${\mathrm{\Gamma }}_{\alpha }\left(ϵ\right)={\displaystyle \sum _{l=1}^{N_l}}\frac{{\gamma }_{\alpha l}^2}{{\left(ϵ-{ϵ}_{\alpha l}\right)}^2+{\gamma }_{\alpha l}^2}{W}_{\alpha l},(21)$

with a corresponding retarded self-energy

${\mathrm{\Sigma }}_{\alpha }^r\left(ϵ\right)=\frac{1}{2}\mathbb{H}\left[{\mathrm{\Gamma }}_{\alpha }\right]\left(ϵ\right)-\frac{i}{2}{\mathrm{\Gamma }}_{\alpha }\left(ϵ\right).(22)$

In Eq. (21), W_αl is a fitting coefficient in a matrix sense, and in Eq. (22), $\mathbb{H}$ denotes the Hilbert transform. This leads to a particularly simple form of the lead self-energies and allows the KBE to be solved exactly, in terms of the density matrix ρ(t) = −iG^<(t, t) and a collection of auxiliary modes that contain the necessary information of the two-time Green’s function. The coupled set of ODEs for propagating the density matrix then becomes [22, 23, 25, 27–30]

$i{\partial }_t\rho \left(t\right)=\left[H\left(t\right),\rho \left(t\right)\right]+i{\displaystyle \sum _{\alpha }}{C}_{\alpha }\left(t\right),(23)$

where the current matrix C_α is defined in Eq. (18)

The memory integral is handled through contour integration and identifying so-called auxiliary modes to obtain the auxiliary mode approach to the problem [23]. The current matrix C_α(t) has several associated equations of motion with coefficients related to the weights and poles in Eq. (21) (23). The system of ODEs, as shown in Eq. (23), can then be solved using a standard ODE solver method. In this particular case, an adaptive Runge–Kutta–Fehlberg fourth-order method (RKF45) was used [31]. This serves as a basis for comparison, as this method is only dependent on having a suitable fit of Γ_α instead of having an approximation of how fast the system can change. Because of the setup, getting the fit of Γ_α for each electrode is just fitting two scalar functions to the form in Eq. (21), which is easily carried out using an equally spaced grid and minimizing the least square error with respect to the fitting coefficient on each Lorentzian.

These steps, from obtaining fits to the lead Γ_α’s to propagating the density matrix and its auxiliary modes in time in a numerically exact way, were implemented in a code that will be published in the future. The fits used can be seen in Figure 4A together with a comparison of the transmission function (tip to graphene) with numerically “exact” self-energies and fitted self-energies in Figure 4B. For the fits, as shown in Figure 4, 100 Lorentzians were used.

FIGURE 4

FIGURE 4. (A) Lorentzian fits to Γ_gr (blue) and Γ_tip (red) used for the AM approach. (B) Tip-graphene transmission function calculated with the numerically “exact” self-energies and the fitted Γ’s without the Hubbard term (U = 0). The transmission was calculated using TBtrans [13].

4 Results

To test our methods, we calculate self-consistently the time-dependent occupation of spin-up and spin-down electrons on the H atom and central C atom, along with the time-resolved spin-polarized current injected into the graphene sheet. A net current may be injected since electrons can tunnel from the STM tip through the hydrogen atom and into the graphene sheet. These properties are calculated using the two methods outlined previously. For our purposes, the AM method may be regarded as numerically exact. In the DE method, the equations of motion are expanded to first order; i.e., we include only first-order dynamical corrections. In the fully adiabatic limit, the self-consistent solution to the electronic dynamics may be found semi-analytically, and for comparison, we include such a calculation in all figures shown.

In the absence of the THz pulse, the system relaxes to an equilibrium configuration which is similar to the ground-state found in the DFT calculations in Section 3. In particular, the ground state is spin polarized when V_DC = 0, but the polarization vanishes upon the inclusion of a sufficiently large bias. We analyze in the following the dynamics for V_DC = 0.5 eV and V_DC = 0.0 eV. These two situations correspond to the situations in the top (0.4 V) and bottom (0.0 V) panels in the left (g = 0) part of Figure 1C, respectively. This way of applying a DC bias to the system mimics that used in the experiment by Cocker et al. [32].

In Figure 5, we show the time-dependent occupation of the hydrogen (A and B) and carbon (C and D) orbitals as determined self-consistently by the DE and the AM method when the system is perturbed by the THz pulse. We take the central frequency of the pulse to be ν₀ = 3 THz and its amplitude to 1 eV, as shown in the experiment by Peller et al. [4]. We will show in the following that the DE method breaks down for this system between 5 and 10 THz, and so 3 THz is one of the highest frequencies for which we would expect the method to work. The gray dotted line marks the adiabatic self-consistent solution, which neglects all dynamics associated with the finite response time of the electronic system. The main new feature when dynamics are included is a delay in the switching time from the polarized to the non-polarized state, and vice versa. The DE and AM methods are in reasonably good agreement with respect to this time delay, in particular for the polarized → non-polarized (P → NP) transition for the DC biased system in Figures 5A,C. Note that for the NP → P transition to occur, the spin-up/down symmetry of the system must be spontaneously broken. Physically, this is achieved by any random (e.g., thermal) fluctuation of the system, but in the numerical implementation, a small spin-up/down asymmetry is introduced explicitly to break this symmetry. We implemented this by introducing a small energy splitting between spin-up and -down electrons, ɛ_H → ɛ_H + ɛ_σ, where ɛ_↑/↓ = ∓1 μeV, i.e., a much smaller energy scale than any other scales in the system. This parameter, therefore, controls the spontaneous symmetry breaking during the NP → P transition without, otherwise, affecting the dynamics.

FIGURE 5

FIGURE 5. (A) Time-dependent occupation of the hydrogen orbital (⟨n_H⟩(t)) as determined self-consistently by the DE and the AM method when the system is perturbed by a THz pulse at a central frequency of 3 THz (middle panels) for a DC bias where the system is non-polarized (V_DC = 0.5 eV). (B) H-occupation for a DC bias yielding a polarized system (V_DC = 0 eV). (C) Carbon p_z-orbital occupation (⟨n_C⟩(t)) for non-polarized and (D) polarized cases. In each case, the lines labeled n_ad mark the adiabatic solution to the problem.

The sign of ɛ_σ was chosen to make the system eventually return to the initial polarization state for the unbiased configuration. One could equally well have chosen the opposite sign, which would make the pulse switch the polarization state of the H adatom for sufficiently low central frequencies. In the absence of an external magnetic field splitting the degeneracy between the spin states and if the central frequency of the pulse is sufficiently low to completely quench the polarization, the final polarization state would be independent of the initial state.

In Figure 5C, we can observe that as the hydrogen atom acquires a spin polarization, the carbon atom acquires the reverse polarization. This happens even though the model does not include electron–electron repulsion in the carbon orbital. It is known that the spin polarization of a hydrogen adsorbate on graphene induces a sublattice-dependent spin polarization in the graphene, extending for several unit cells around the adsorption site [9, 16]. Our observation is qualitatively in line with these results, though the effect could be modeled more accurately by using a spin-dependent coupling self-energy for the graphene-device contact. One could also include more carbon atoms in the device region and, thus, explicitly model their time-dependent spin polarization. The AM results for the carbon orbital were shifted by a constant vertical off-set to make all methods agree in the steady-state regime. If such an offset is not included, a small discrepancy is observed due to the Lorentzian fit to the self-energy used in the AM method in place of the exact self-energy. The carbon orbital, which couples directly to the graphene sheet through the self-energy, is particularly sensitive to any variation in this parameter. We remark that the quality of the fit may always be improved by including more Lorentzian functions, until a satisfactory fit is achieved.

Figure 6 shows the time-resolved spin-polarized current J_s = J_↑ − J_↓ injected into the graphene due to the THz pulse for the DC biased and unbiased configurations. Large peaks are observed when the magnetization state switches, corresponding to the abrupt change in occupation of the hydrogen atom. These peaks are completely absent in the adiabatic calculation. The DE result displays rapid small-amplitude oscillations shortly before the switch occurs. This is a numerical artifact which reflects the fact that the self-consistent cycle did not converge for these particular time steps. This is not surprising as the DE approach is based on the assumption that the time-evolution is near adiabatic, which is far from the case near the transition, see the discussion in Supplementary Material. Despite this challenge, the overall prediction by the DE method is in good agreement with the exact (AM) result, especially for the biased (NP → P → NP) configuration.

FIGURE 6

FIGURE 6. Spin-polarized current J_s = J_↑− J_↓ into the graphene layer for (DC) biased (A) and unbiased (B) configurations. A positive sign of J_s means that there is a net spin-up current into the graphene sheet. Both time-dependent methods predict large spikes in the current when the polarization state of the hydrogen atom changes, in particular during the transition from non-polarized to polarized currents. This feature is absent in the adiabatic calculation.

The spectrum of the spin-polarized current is shown in Figure 7. The DE and AM spin-polarized currents contain significant Fourier components at much higher frequencies than the adiabatic current, as one may have anticipated from the peaks observed in Figure 6. They both also predict oscillations in the spectra at high frequencies (most evident in Figure 7A), although the amplitudes and exact locations of these oscillations are different. This finding indicates that the DE method captures well the qualitative effects of the dynamics, but it may not be quantitatively correct on all accounts. All three methods predict currents whose spectra contain higher frequencies than the incident THz pulse, attesting to the highly non-linear behavior of the system.

FIGURE 7

FIGURE 7. Spectrum of the spin-polarized current J_s = J_↑ − J_↓ into the graphene layer for (DC) biased (A) and unbiased (B) configurations. Due to the spike in the current associated with the spin-polarization phase transition, the spectrum of the current contains much larger frequencies than that of the external pulse, an indication that the problem is highly non-linear. This is especially evident in the biased case, when the non-polarized → polarized transition is particularly rapid. Note the different scales of both axes in the two figures; the black curve showing the spectrum of the external pulse in arbitrary units is the same in both figures.

As the DE method is exact in the adiabatic limit, a natural question is how fast the externally imposed time-dependence may be before the method breaks down. Figure 8 shows the total spin injected into the graphene sheet $Q_s({\nu }_0)=\frac{1}{2}{\int }_{-\infty }^{\infty }{J}_s(t,{\nu }_0)\mathrm{d}t$ due to a single THz pulse as a function of the central frequency ν₀ of the pulse.

FIGURE 8

FIGURE 8. Total spin injected into the graphene layer multiplied by pulse frequency plotted for various frequencies for (DC) biased (A) and unbiased (B) configurations. The total spin injected is calculated as $Q_s=\frac{1}{2}\int (J_s)\mathrm{d}t$ with J_s given by three different methods. Note the different scales of the y-axes on the two figures. In the adiabatic case, the transferred charge is simply inversely proportional to the frequency, and so ν₀Q_s,ad = constant. Multiplying by the frequency, thus, emphasizes the dynamical variation of the transferred charge as the frequency of the external potential is varied.

Neglecting dynamical effects (i.e., the adiabatic case), the transferred spin is inversely proportional to the central frequency of the pulse because the current is frequency-independent. To highlight the difference due to the dynamical corrections, we plot ν₀Q_s which is constant for the adiabatic calculation and is shown in Figure 8 as a gray dashed line. At low frequencies, the DE and AM methods are in good agreement, predicting that the transferred charge is roughly independent of ν₀, which is observed in Figure 8 as a linear relationship between ν₀Q_s and ν₀. At low frequency, the system has time to completely polarize/depolarize, resulting in an approximately frequency-independent current associated with the filling and emptying of the device orbitals. As the frequency increases the DE method begins to deviate. From the figures, we assess that ν₀ ≈ 5 THz is the largest frequency for which the DE works well. We remark that at high frequencies, the self-consistent cycle used in the DE method does not converge well. For these frequencies, the dynamical density with the smallest error obtained in the self-consistent cycle is included in the figure. The dots show the frequencies at which calculations were made, and the lines interpolate between these points. In the biased case (Figure 8A), the dynamical behavior predicted by the AM method changes significantly around ν₀ = 10 THz. This change is related to the fact that at high frequencies, the electronic system cannot adjust to follow the pulse. Thus, when the system starts out in the NP state, it does not have time to polarize. For a fully non-polarized system, the spin-polarized charge will always vanish due to the spin-up/down symmetry. This explains the sudden decrease in the transferred charge at this point. This highly non-adiabatic effect is not at all captured by the DE method.

5 Conclusion and outlook

We have tested the computationally efficient DE/time-scale separation method, introduced in [8] on a model of an experimentally relevant system, namely, mean-field, open-system calculations of charge, and spin dynamics in a THz-STM junction involving a hydrogen adsorbate on graphene. We demonstrated that for low THz frequencies, the DE method gives a good description of the non-adiabatic dynamics compared to full propagation using the AM method. This calculation may, furthermore, be seen as a toy model for a self-consistent TD-DFT-NEGF setup, where dynamical corrections to the density due to an external field are taken into account via the dynamical expansion method. Some of the numerical challenges and methods are discussed in Supplementary Material.

We have shown how the THz pulse generates spin-dynamics in the hydrogen on the graphene THz-STM system with higher harmonic frequencies appearing due to the highly non-linear system. This dynamics may affect both the tunnel current from the tip electrode as well as the spin-dependent scattering of carriers in graphene, as observed by the spin-dependent scattering cross section of the hydrogen adsorbate.

The magnetic structure in the graphene induced by the polarized hydrogen atom is long ranged [9]. As an outlook, it would be interesting to extend the computational region and possibly to include dynamics among multiple localized spins. We also note that the carrier-envelope phase (CEP) may be changed so the THz waveform changes in a continuous way, e.g., from maximum being positive to negative [33], as a further way of tuning the dynamics.

The expansion in powers of the central time derivative ∂_T used in the DE may in principle be continued up to arbitrary order to include dynamical effects to higher degrees of accuracy. It would be interesting to extend the numerical implementation beyond the first-order dynamical corrections. In a preliminary study, we found that it was difficult to converge the self-consistent cycle within the second-order DE, and we did not pursue this further. An open question for the DE method is the issue of convergence. It is not trivial to quantify how slow the external time variation in a given system should be in order for the lowest-order terms to approximate the full dynamics. A further investigation into this issue would be highly interesting, in particular, in combination with the inclusion of higher-order terms in the formal series expansion.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

JKS implemented the DE method and performed the DE and some of the AM calculations. ABL implemented the AM method and performed some of the AM calculations. FG performed the DFT-NEGF calculations. All authors contributed to the article and approved the submitted version.

Funding

JKS, FG were supported by Villum fonden (Project No. VIL 13340), ABL by the Independent Research Fund Denmark (Project No. 0135-00372A), and APJ by the Center for Nanostructured Graphene (CNG), which is sponsored by the Danish National Research Foundation (Project No. DNRF103).

Acknowledgments

The authors thank Profs. Peter Uhd Jepsen and Peter Bøggild for stimulating discussions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2023.1237383/full#supplementary-material

References

1. Ammerman SE, Jelic V, Wei Y, Breslin VN, Hassan M, Everett N, et al. Lightwave-driven scanning tunnelling spectroscopy of atomically precise graphene nanoribbons. Nat Commun (2021) 12:6794. doi:10.1038/s41467-021-26656-3

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Cocker TL, Jelic V, Hillenbrand R, Hegmann FA. Nanoscale terahertz scanning probe microscopy. Nat Photon (2021) 15:558–69. doi:10.1038/s41566-021-00835-6

CrossRef Full Text | Google Scholar

3. Wang L, Xia Y, Ho W. Atomic-scale quantum sensing based on the ultrafast coherence of an h2 molecule in an stm cavity. Science (2022) 376:401–5. doi:10.1126/science.abn9220

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Peller D, Roelcke C, Kastner LZ, Buchner T, Neef A, Hayes J, et al. Quantitative sampling of atomic-scale electromagnetic waveforms. Nat Photon (2021) 15:143–7. doi:10.1038/s41566-020-00720-8

CrossRef Full Text | Google Scholar

5. Haug H, Jauho AP. Quantum kinetics in transport and optics of semiconductors. Berlin, Germany: Springer (2008).

Google Scholar

6. Stefanucci G, Van Leeuwen R. Nonequilibrium many-body theory of quantum systems: a modern introduction. Cambridge, United Kingdom: Cambridge University Press (2013).

Google Scholar

7. Bonitz M. Quantum kinetic theory. Berlin, Germany: Springer (2016).

Google Scholar

8. Honeychurch TD, Kosov DS. Timescale separation solution of the Kadanoff-Baym equations for quantum transport in time-dependent fields. Phys Rev B (2019) 10:245423. doi:10.1103/PhysRevB.100.245423

CrossRef Full Text | Google Scholar

9. González-Herrero H, Gómez-Rodríguez JM, Mallet P, Moaied M, Palacios JJ, Salgado C, et al. Atomic-scale control of graphene magnetism by using hydrogen atoms. Science (2016) 352:437–41. doi:10.1126/science.aad8038

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Novoselov KS, Geim AK, Morozov SV, Jiang D, Zhang Y, Dubonos SV, et al. Electric field effect in atomically thin carbon films. Science (New York, N.Y.) (2004) 306:666–9. doi:10.1126/science.1102896

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Brandbyge M, Mozos JL, Ordejón P, Taylor J, Stokbro K. Density-functional method for nonequilibrium electron transport. Phys Rev B (2002) 65:165401. doi:10.1103/physrevb.65.165401

CrossRef Full Text | Google Scholar

12. Soler JM, Artacho E, Gale JD, García A, Junquera J, Ordejón P, et al. The siesta method for ab initio order-n materials simulation. J Phys Condensed Matter (2002) 14:2745–79. doi:10.1088/0953-8984/14/11/302

CrossRef Full Text | Google Scholar

13. Papior N, Lorente N, Frederiksen T, García A, Brandbyge M. Improvements on non-equilibrium and transport green function techniques: the next-generation transiesta. Comp Phys Commun (2017) 212:8–24. doi:10.1016/j.cpc.2016.09.022

CrossRef Full Text | Google Scholar

14. Papior N. Zerothi/sisl: v0.13.0. Zenodo (2022) 12:2. doi:10.5281/zenodo.597181

CrossRef Full Text | Google Scholar

15. Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Phys Rev Lett (1996) 77:3865–8. doi:10.1103/physrevlett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Gao F, Zhang Y, He L, Gao S, Brandbyge M. Control of the local magnetic states in graphene with voltage and gating. Phys Rev B (2021) 103:L241402. doi:10.1103/PhysRevB.103.L241402

CrossRef Full Text | Google Scholar

17. Papior N, Gunst T, Stradi D, Brandbyge M. Manipulating the voltage drop in graphene nanojunctions using a gate potential. Phys Chem Chem Phys (2016) 18:1025–31. doi:10.1039/c5cp04613k

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Saloriutta K, Uppstu A, Harju A, Puska MJ. Ab initio transport fingerprints for resonant scattering in graphene. Phys Rev B (2012) 86:235417. doi:10.1103/physrevb.86.235417

CrossRef Full Text | Google Scholar

19. Markussen T, Rurali R, Cartoixa X, Jauho AP, Brandbyge M. Scattering cross section of metal catalyst atoms in silicon nanowires. Phys Rev B (2010) 81:125307. doi:10.1103/physrevb.81.125307

CrossRef Full Text | Google Scholar

20. Uppstu A, Saloriutta K, Harju A, Puska M, Jauho AP. Electronic transport in graphene-based structures: an effective cross-section approach. Phys Rev B (2012) 85:041401. doi:10.1103/physrevb.85.041401

CrossRef Full Text | Google Scholar

21. Papior N, Calogero G, Leitherer S, Brandbyge M. Removing all periodic boundary conditions: efficient nonequilibrium green’s function calculations. Phys Rev B (2019) 100:195417. doi:10.1103/physrevb.100.195417

CrossRef Full Text | Google Scholar

22. Croy A, Saalmann U. Propagation scheme for nonequilibrium dynamics of electron transport in nanoscale devices. Phys Rev B - Condensed Matter Mater Phys (2009) 80:245311. doi:10.1103/PhysRevB.80.245311

CrossRef Full Text | Google Scholar

23. Popescu BS, Croy A. Efficient auxiliary-mode approach for time-dependent nanoelectronics. New J Phys (2016) 18:093044. doi:10.1088/1367-2630/18/9/093044

CrossRef Full Text | Google Scholar

24. Kershaw VF, Kosov DS. Non-equilibrium Green’s function theory for non-adiabatic effects in quantum transport: inclusion of electron-electron interactions. J Chem Phys (2019) 150:074101. doi:10.1063/1.5058735

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xie H, Jiang F, Tian H, Zheng X, Kwok Y, Chen S, et al. Time-dependent quantum transport: an efficient method based on liouville-von-neumann equation for single-electron density matrix. J Chem Phys (2012) 137:044113. doi:10.1063/1.4737864

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zheng X, Chen G, Mo Y, Koo S, Tian H, Yam C, et al. Time-dependent density functional theory for quantum transport. J Chem Phys (2010) 133:114101. doi:10.1063/1.3475566

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zheng X, Wang F, Yam CY, Mo Y, Chen G. Time-dependent density-functional theory for open systems. Phys Rev B (2007) 75:195127. doi:10.1103/physrevb.75.195127

CrossRef Full Text | Google Scholar

28. Jauho AP, Wingreen NS, Meir Y. Time-dependent transport in interacting and noninteracting resonant-tunneling systems. Phys Rev B (1994) 50:5528–44. doi:10.1103/physrevb.50.5528

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Hu J, Luo M, Jiang F, Xu RX, Yan Y. Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems. J Chem Phys (2011) 134:244106. doi:10.1063/1.3602466

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Jin J, Zheng X, Yan Y. Exact dynamics of dissipative electronic systems and quantum transport: hierarchical equations of motion approach. J Chem Phys (2008) 128:234703. doi:10.1063/1.2938087

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Hairer E, Nørsett SP, Wanner G. Solving ordinary differential equations. 1, Nonstiff problems. Berlin, Germany: Springer-Vlg (1993).

Google Scholar

32. Cocker TL, Peller D, Yu P, Repp J, Huber R. Tracking the ultrafast motion of a single molecule by femtosecond orbital imaging. Nature (2016) 539:263–7. doi:10.1038/nature19816

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Ludwig M, Kazansky AK, Aguirregabiria G, Marinica DC, Falk M, Leitenstorfer A, et al. Active control of ultrafast electron dynamics in plasmonic gaps using an applied bias. Phys Rev B (2020) 101:241412. doi:10.1103/PhysRevB.101.241412

CrossRef Full Text | Google Scholar