Green's function methods for single molecule junctions

We present a brief pedagogical review of theoretical Green's function methods applicable to open quantum systems out of equilibrium in general, and single molecule junctions in particular. We briefly describe experimental advances in molecular electronics, then discuss different theoretical approaches. We then focus on Green's function methods. Two characteristic energy scales governing the physics are many-body interactions within the junctions, and molecule-contact coupling. We therefore discuss weak interactions and weak coupling, as two limits that can be conveniently treated within, respectively, the standard nonequilibrium Green's function (NEGF) method and its many-body flavors (pseudoparticle and Hubbard NEGF). We argue that the intermediate regime, where the two energy scales are comparable, can in many cases be efficiently treated within the recently introduced superperturbation dual fermion approach. Finally, we review approaches for going beyond these analytically accessible limits, as embodied by recent developments in numerically exact methods based on Green's functions.


I. INTRODUCTION
Since the first theoretical proposal to use single molecules as electronic devices 1 and the first experimental realization of a single molecule junction, 2,3 the field of molecular electronics has made tremendous progress. This is evidenced by improved nanoscale fabrication techniques and a significant increase in the variety of signals measurable in single molecule junctions. 4 As experimental techniques developed, the focus of research has shifted over the years from measurements of elastic coherent transport in junctions 2,5,6 to studies of quantum coherence effects 7-10 and to improving the stability and reproducibility of measurements. [11][12][13][14][15] Increasingly sensitive measurements allowed for detection of inelastic effects in the off-resonant tunneling regime, leading to the appearance of inelastic electron tunneling spectroscopy (IETS). [16][17][18][19][20][21][22][23][24][25][26] IETS is particularly useful as a diagnostic tool, since it provides vibrational fingerprints of molecules within the junction. 17 It is also instrumental for imaging of molecular structure and chemical bonding, 21,23 as well as for probing intermolecular interactions and anharmonic overtones. 25,26 Inelastic measurements in the resonant regime (RI-ETS) can also be performed, where the inelastic signal reveals itself as a peak in conductance that indicates positions of vibronic molecular levels. [27][28][29][30] The electronic population fluctuations within the molecule in this regime lead to stronger electron-vibration interactions, which then leads to destruction of coherence 31 and to reorganization of polarization (polaron formation) in the junction. 32 The former also controls the transition from tunneling-dominated transport to the hopping a) Electronic mail: gcohen@tau.ac.il b) Electronic mail: migalperin@ucsd.edu regime, [33][34][35][36][37][38][39][40][41][42][43] while the latter may result in, e.g., conformational changes in molecular structure. This has been considered as a possible mechanism for highly nonlinear current-voltage characteristics such as negative differential resistance (NDR) and hysteresis. 32,[44][45][46][47][48][49][50][51][52][53][54] A related body of research involves molecular heating/cooling driven by the electronic flux [55][56][57] and thermoelectricity. [58][59][60][61][62][63][64][65] Heating is directly relevant to questions regarding junction stability and the reproducibility of measurements. Nanoscale thermoelectrics raises interesting potential technological possibilities, such as utilization of quantum effects 66,67 for constructing highly efficient nanoscale devices for energy transduction.
Historically, IETS was the first spectroscopic tool to be applied to single molecule junctions. However, developments in laser techniques also allowed for performing optical spectroscopy measurements in currentcarrying molecular junctions. This combination of optical spectroscopy and molecular electronics lead to the appearance of a new field of research coined molecular optoelectronics. 68 In particular, surface enhanced Raman spectroscopy (SERS) [69][70][71][72] in junctions, besides providing (complementary to IETS) information on molecular vibrations, allows for estimating bias-induced heating of electronic and vibrational degrees of freedom. 73,74 Tipenhanced Raman spectroscopy (TERS) 75,76 yields information on molecular structure. [77][78][79][80][81] Furthermore, optical emission in biased junctions was observed as biasinduced luminescence. [82][83][84][85][86][87] In this context, electroluminescence was employed to study energy transfer [88][89][90] and as a tool for molecular imaging with submolecular resolution. 91 It was also employed for designing energetically efficient light emitting diodes. 92 Recently, strong light-matter interaction was measured in single molecule nanocavities. 89,[93][94][95][96][97][98] While measurements in the strong coupling regime were performed in the absence of electron flux, extensions to current-carrying molecular junc-tions are expected soon. Such developments will be an important step forward in the quest for optical control and characterization of molecular junctions.
Molecular spintronics is another rapidly growing field of research closely related to construction of nanodevices for quantum information. 99 Here, the focus is on the flow of spin through the junction, rather than that of charge. Measurements of spin polarized currents, 100,101 spin-flip IETS, 102,103 spin interactions, 104-106 molecular spin selectivity, [107][108][109][110][111] spin valves, 112,113 and photospintronic effects 114,115 have all been reported in the literature.
In strongly interacting junctions at low temperatures, spin exchange interactions can dress unpaired spins in the molecule with a cloud of opposite spins in the contacts. The result is a correlated many-body electronic state comprising electrons in a physical region that can be much larger than that of the molecule. This correlation reveals itself as a zero-bias conductance peak, known as the Kondo or Abrikosov-Suhl resonance, and has been observed in transport measurements in single atom 116 and single molecule [117][118][119][120][121] transistors. The behavior of the Kondo effect far from equilibrium (beyond linear response) is still not fully understood, but its signature has been observed in carbon nanotubes 122 and πconjugated (OPV-5) molecules. 123 Vibrational sidebands to the Kondo peak, 124,125 as well as an IETS-RIETS transition in the presence of a Kondo resonance, 124,126,127 have also been reported in the literature.
Finally, modern measurement techniques allow access not only to the electronic flux, but also to its noise or variance, and sometimes higher moments. [128][129][130][131][132] These quantities are related to the cumulants of the full counting statistics (FCS) generating function of electronic transport. [133][134][135][136] Noise and higher moments can provide information not accessible via flux measurements, e.g. the number of transport channels in the junction and the effective charge of carriers. 136 This paradigm has important implications: for example, noise measurements could potentially yield information on entanglement between degrees of freedom in the system [137][138][139] or serve as an important ingredient in experimental validation for theories of nanoscale quantum thermodynamics. 140 Today, molecular electronics is a thriving and interdisciplinary area of research involving researchers from fields as diverse as condensed matter physics and statistical mechanics; nonlinear optics, plasmonics and nanoscience; quantum chemistry and chemical dynamics; and a wide variety of engineering disciplines. Single molecule junctions are used as testbeds for the study of fundamental physical properties of matter at nanoscale. They are also widespread in more applied contexts, as models for quantum devices enabling energy transduction and storage, and as potential elements in proposed architectures for classical and quantum computers. [141][142][143][144][145] Throughout the evolution of molecular electronics, the tremendous progress in experimental techniques has presented theorists with increasingly challenging and com-plex questions. Molecular junctions are open quantum systems that can be driven far from their equilibrium state, and are characterized by a plethora of interactions at widely varying relative magnitudes. This makes them attractive candidates for exploring a rich variety of fundamental quantum many-body physical problems in limits that remain poorly understood. This physics is very different from better understood limits that can be treated by standard approximations built around mean field arguments, which implies that rigorous theoretical treatments must rely on the development and implementation of specialized and advanced new methods. Here, we present a short overview of recent developments in the Green's function techniques for single molecule junctions. Our consideration is mostly focused on method development. We note that a complementary review focused on applications was recently published in Ref. 146.
The rest of our paper is structured as follows: in Section II we give a short overview of (non-Green's function based) theoretical methods utilized in transport problems. Section III specifically focuses on Green's function methods for molecular electronics. We discuss the standard nonequilibrium Green's function (NEGF) approach formulated around the noninteracting limit, then its many-body flavors, pseudoparticle (PP-) and Hubbard NEGF. We then highlight recent developments in the superperturbation approach to transport, briefly review the state-of-the-art on time dependent problems, and go on to describe numerically exact quantum Monte Carlo approaches to accessing nonequilibrium Green's functions. Section IV summarizes the review.

II. OVERVIEW OF THEORETICAL METHODS
Before discussing theoretical methods, we begin by defining a concrete model of a molecular junction. We consider a molecule M coupled to a set of baths B. The molecule may also be subjected to external (classical) driving. Baths usually encompass fermionic electronic contacts or leads, modeled as reservoirs of free charge carriers, each of which is taken to be in its own equilibrium state. They can also include bosonic environments representing phonons and/or the optical modes of a quantized radiation field (see Fig. 1). The Hamiltonian is taken to beĤ whereĤ M (t) andĥ B are the molecular and bath Hamiltonians, respectively.V M B (t) describes coupling between the molecule and baths. We note that the molecular Hamiltonian may describe several degrees of freedom (e.g., electronic and vibrational) and coupling between them, and can be interacting. The baths are taken to be noninteracting, in the sense that their Hamiltonian is quadratic in each bath's set of fermionic/bosonic creation and annihilation operators. The molecule-bath couplinĝ V M B (t) can also be interacting, for example when it describes electron-phonon coupling. For most purposes in the field, it is sufficient to consider a bath Hamiltonian h B that is both noninteracting (i.e. describes a normal metal or semiconductor) and time independent. We note that generalization to include time-dependent driving in the baths (e.g. ac bias) is straightforward and for Green's function techniques follows the celebrated work by Jauho, Wingreen and Meir. 147 A few simple cases are of special interest. Below, if the molecular HamiltonianĤ M (t) and the couplingV M B (t) are both quadratic, we refer to the entire system as noninteracting. One canonical case that will be referred to below is the Anderson impurity model. Here, only interacting term is of the form Un M ↑nM ↓ , wheren M ↑(↓) is a a spin up (down) population operator on the molecule; and there is no inter-spin coupling. Another is the Holstein impurity model, where a single spinless electronic band is coupled to a single vibration or to a bath of phonons by a term of the form k σ∈↑,↓ λ σnσ b k +b † k . Here, b † k (b k ) creates (annihilates) a phonon in mode k. If both interactions are present, the model may be called an Anderson-Holstein model.
Theoretical approaches to open quantum systems far from equilibrium differ in their treatment of both the system and bath degrees of freedom, but most rely either formally or in practice on time evolution from a convenient initial state to a more interesting, possibly stationary nonequilibrium state. In particular, wavefunction methods consider the unitary time evolution generated by the Schrödinger equation-complemented with appropriate boundary and initial conditions-in the entire system, comprising both the molecular and bath regions. Density matrix approaches integrate out bath degrees of freedom, leading to a generally non-Markovian effective dynamics of the reduced density matrix as a function of time. This encapsulates the effect of the baths into a modified nonunitary equation of motion known as a generalized quantum master equation. Finally, Green's function methods consider the effective equations of motion of correlation functions rather than those of singletime properties. At the very least, two-time correlations are used, and higher order, many-time correlations may also be evaluated and employed. All three formulations are equivalent in the exact case, but most studies rely on approximations; here, these different languages lend themselves to very different schemes. In what follows we will discuss the relative merits of Green's function approaches in the study of molecular junctions. First, we briefly consider several important points regarding wavefunction and density matrix methods, without providing an exhaustive review. Green's function approaches, our main focus, will be treated in more detail in Section III.

A. One-body wavefunction methods
In stationary transport problems wavefunction formulations are often combined with scattering theory. Here, the set of scattering states are obtained by solving the Lippmann-Schwinger equation. The solution is most often performed for an effective single-particle wavefunction. This means that the many-body wavefunction is assumed to be a Slater determinant of single-particle orbitals, often referred to as a single reference state. For noninteracting systems, where the many-body problem can be reduced to single particle description, such considerations are exact; in interacting systems, they become approximate. For example, it is common to use the solutions of the Hartree-Fock or Kohn-Sham equations as effective single-particle orbitals. The celebrated Landauer-Büttiker formalism and its many extensions describes the application of noninteracting wavefunction methods to transport. In one set of representative examples from the literature, a combination of density functional calculations for molecular electronic structure and wavefunction-based scattering theory was used in the simulation of elastic transport in junctions. [148][149][150] Single particle wavefunction methods, especially in the scattering state representation, are remarkably simple and computationally efficient. Where the single particle picture remains accurate, they can be generalized in a straightforward manner to treating a large variety of interactions and observables. These attractive properties have resulted in a massive popularity, and such methods are now ubiquitous.
For interacting systems in the Fermi liquid regime, at linear response from zero temperature equilibrium, an exact noninteracting picture is available that allows e.g.
for the evaluation of currents and shot noise. 151 While a few other examples of this kind may exist, in general many-body interactions can be handled only approximately within the noninteracting framework. For example, inelastic transport can be addressed by treating the electron-vibration coupling perturbatively within the Born approximation; 152 and inelastic current noise has also been treated within a similar set of considerations. 153,154 A method treating electronvibration interactions exactly in single electron scattering problems was introduced in Refs. 155 and 156 and later applied to the description of cooperative effects in inelastic tunneling. 157 Single particle scattering theory was also employed in the derivation of current-induced forces on molecular nuclei [158][159][160] and in the quantum thermodynamics context. [161][162][163] The relationship between single particle scattering-matrix and nonequilibrium Green's function formalisms for noninteracting systems is explained in detail in Refs. 164 and 165. Another set of approaches based to some degree on the ability to understand single-particle scattering states relies on an exact mapping of steady-state nonequilibrium to an effective equilibrium state, 166 which was shown to be equivalent to the Zubarev's nonequilibrium statisti-cal operator method. 167 This mapping was utilized in studies of transport in the resonant level, 168 Anderson impurity 169 and Anderson-Holstein 170 models. An interesting development in this regard is that nonequilibrium steady state transport properties can be obtained from imaginary time quantities, albeit with the need for analytical continuation. 171,172 B. Many-body wavefunction methods Where interactions are strong enough that the singleparticle picture no longer holds, it is necessary to explicitly consider wavefunctions in the full many-body Hilbert space. The size of this (multireference) space increases exponentially with the number of electrons in the system, which is always infinite in scenarios describing electronic transport across junctions. Therefore, use of many-body wavefunctions necessarily requires efficient truncated representations of the Hilbert space in practice. A variety of celebrated numerical wavefunction schemes have employed some variation of this idea to address transport problems, including the numerical renormalization group (NRG) method, 173,174 matrix product state (MPS) techniques, 175 multiconfiguration Hartree-Fock (MCTDH) and its later, [176][177][178] and more recent modifications. 179,180 We will very briefly go over some of these numerically exact approaches.
NRG relies on extracting a small, effective many-body Hamiltonian from the infinite interacting problem by way of a normalization procedure. Nonequilibrium problems can be addressed directly in steady state by considering scattering states, 174 or by following the time evolution from some initial state. 173 Variations of this idea have been applied to steady state 181,182 and transient properties 173,183,184 of strongly correlated junctions.
Applications of MPS techniques to transport began to appear in the literature before that term came to be widely used, initially in the form of time-dependent density matrix renormalization group (DMRG) calculations. [185][186][187][188][189] Shot noise and FCS within the self-dual interacting resonant level model (an exactly solvable model for interacting transport) have also been evaluated. [190][191][192] In the language that has emerged over the years, these simulations rely on the MPS, an efficient and accurate ansatz for the ground state of onedimensional systems, and understanding this has lead to concrete algorithmic improvements. [193][194][195][196] It is possible to propagate the complete system (M + B) in time after an initial quench, such as the sudden application of a bias voltage; or to directly find an MPS representing a nonequilibrium steady state 197,198 . The main limitation of MPS methods is their inability to describe states with long-ranged entanglement along the chain. This often makes it difficult to reach long simulation times, though this restriction can be greatly ameliorated by a good choice of basis. [199][200][201][202][203][204] Higher dimensional systems also pose a challenge to MPS methods and other tensor network techniques, even in equilibrium.
Multilayer, multiconfiguration time-dependent Hartree (ML-MCTDH) theory employs a time-dependent multiconfigurational expansion of wavefunctions that are propagated in time by using the Dirac-Frenkel variational princle. The method, originally formulated in Ref. 205 for distinguishable particles, was extended to account for bosonic and fermionic statistics in Ref. 176. The method has been applied to heat transport 206 and electronic transport in junctions with electron-electron 207 and electron-phonon interactions. 177,178,208,209 In the latter case, the possibility of bistability in nonequilibrium quantum systems coupled to vibrational degrees of freedom 210 was explored. 211,212 C. Density matrix formulations The density matrix ρ of the full system, including both the molecular region M and the baths B, evolves in time according to the Liouville-von-Neumann equation. However, when the wavefunction is determined by N coefficients, the density matrix is determined by N 2 ; and in a many-body system, N is exponential in the physical system size. This means that dealing with density matrices of the full system is generally more expensive than dealing with wavefunctions. However, it is possible to consider the reduced density matrix σ ≡ Tr B {ρ}, which has a dimensionality determined entirely by the size of the molecular region M . We will broadly refer to methods formulated in terms of the reduced density matrix as "density matrix methods".
Diagonal elements of σ yield information about probabilities to observe molecular states, while off-diagonal elements relate to bath-induced coherences between such states. Its exact dynamics is given by the Nakajima-Zwanzig-Mori generalized quantum master equation (GQME), 213 but are often treated perturbatively in V M B . 214, 215 The lowest (second) order expansion in the molecule-bath coupling, along with a Markovian assumption, is known in the quantum transport context as the Born-Markov or Redfield approximation. 216 The latter is not adequate in representing bath-induced coherences, 217,218 and is problematic in the presence of degeneracies. 219 The Redfield QME is often written within an additional rotating wave/secular approximation (RWA), where dynamics of the populations (or the Pauli equations) is decoupled from that of the coherences. 216 This has the advantage of guaranteeing positive definite reduced dynamics. The Pauli equations are formulated in the molecular eigenbasis with rates given by the Fermi golden rule.
Markovian Pauli equations have been successfully employed in far too many contexts to provide an exhaustive list here. One interesting example is the modeling of transport in donor-bridge-acceptor molecular complexes; [220][221][222][223][224] another is the description of inelastic transport. 225,226 Many Redfield QME studies also go be-yond the level of Pauli equations, providing some account of quantum coherence. [227][228][229] Nevertheless, while attractive due to their simplicity, Markovian Redfield QMEs fail to account for molecule-bath correlations and hybridization, and therefore cannot be reliably employed in all regimes.
We note in passing that similar considerations employing bare perturbation theory in the molecule-bath coupling are at the heart of nonlinear optical spectroscopy studies. [230][231][232][233][234] While completely justifiable in isolated systems or for strictly Markovian open system dynamics, utilization of bare perturbation theory in open systems with essentially non-Markovian dynamics is known to violate conservation laws, 235,236 and bare perturbative corrections can lead to qualitatively incorrect predictions. 237 Recently, a conserving flavor of nonlinear optical spectroscopy methods based on GF analysis was presented in Ref. 238. We also note that treatments of dynamics involving a combination of non-Markovian and kinetic schemes have been employed in the literature. [239][240][241][242] Such mixing requires care, as it may otherwise lead to qualitative failures in the analysis 243 .
Using the rigorous GQME as a theoretical framework, it is possible to perform expansions in the molecule-bath couplingV M B to higher orders. As with any perturbation theory, the semianalytical implementation becomes expensive for high-order self energies, but is extremely efficient and powerful when low orders are appropriate. Rules for the construction of Feynman diagrams and their resummation or dressing, as well as an application of the resulting theory to transport processes, were presented in Ref. 244. This machinery was then successfully applied to the description of various transport problems. [245][246][247] A simplification resulted in formulation of a kinetic equation approach to transport, 248 where rates were calculated within fourth order perturbation theory. Similar fourth order considerations were employed in Ref. 249, and later developments allowed for a reduction in the number of diagrams. 250 We note that the Liouville space formulation is complicated by the need to define diagrams with respect to ordering on both the Keldysh contour and the real time axis. Green's function analysis, where ordering along the real axis can be dropped, allows some simplifications with respect to this. 251

D. Numerically exact approaches to the density matrix
Most numerically exact methods and many approximate ones employ time propagation and are limited to short times, in the sense that the computational cost of simulating dynamics scales more than linearly with time. Another application of the exact GQME leverages the rapid decay of the memory kernel in time that is characteristic to certain physical regimes. 252 If the memory kernel for the reduced dynamics of the molecular region M can be simulated up to the time where it becomes negligible, the GQME can then be used to evaluate dy-namics to any later time. An extreme example of this is the Markovian limit, where evaluation of the kernel over an infinitesimal timescale suffices to fully describe the local population dynamics at any later time. Evaluating the memory kernel is nontrivial, and one successful scheme, pioneered in Refs. 253 and 254 for the spin-boson model, relies on the evaluation of a set of auxiliary observables and the solution of an integral equation involving them. This scheme was applied to transport through a junction with electron-electron interactions using real time Monte Carlo; 252,255 to inelastic vibrational transport using MCTDH; 208,211 and as a way of enhancing semiclassical simulations that are more accurate at short timescales. 215,256 Alternative schemes for obtaining the memory kernel without the need for evaluating auxiliary observables later emerged, 257,258 though with an added requirement for higher-order derivatives of the dynamics that can sometimes be numerically problematic. A discrete version of these ideas can be expressed in terms of transfer tensors, or dynamical maps that have a matrix product form. 259 Next, we mention two very successful numerically exact techniques that explicitly employ the density matrix picture in their construction: the hierarchical equation of motion (HEOM, sometimes also called the hierarchical quantum master equation or HQME in the literature), and the iterative summation of path integrals (ISPI).
The HEOM was originally introduced to address dissipative chemical dynamics, 260,261 but was later generalized to models of quantum transport. 262,263 Beginning with an exact solution at the atomic limit, the idea is to write a systematic expansion in terms of a series of increasingly high-order auxiliary density matrices that contain information about mixed molecule-bath quantities. This hierarchy can be efficiently evaluated to very high orders for certain forms of the lead density of states such as, for example, a linear combination of several Lorentzian functions. The expansion is eventually truncated at some finite order when a desired level of accuracy is reached. In most cases, many of the auxiliary density matrices have little or no effect and can be dropped.
Early versions of HEOM were limited to high temperatures and specific band structures, 264 but the approach has since undergone rapid progress. [264][265][266][267][268] Several works have considered time-dependent transport characteristics, [268][269][270] and HEOM has the important advantage over most numerically exact techniques that computational scaling in time is linear. New developments enable trading some of this advantage for greatly improved access to more general band structures and lower temperatures. 271 It should soon become clear to what degree these novel techniques enable access to correlated physics at low temperatures.
The ISPI method, 272,273 like the HEOM, traces its roots to a predecessor in chemical dynamics. [274][275][276] This method is based on expressing the action of the time evolution operator e −iĤt on an initial density matrix and measured operator as a path integral over a set of in-termediate molecular states at a set of discrete times. The contribution of each path, which is analogous to the Feynman-Vernon influence functional, is then truncated beyond a finite memory time. Using this assumption, the ISPI method expresses time dependent observables in such a way that the dominant computational cost becomes exponential in the ratio between the memory time and the time discretization interval. This becomes numerically exact for small time discretization when the truncation is applicable.
ISPI methods have been used to explore fermionic transport through quantum junctions in the presence of many-body interactions, in the Anderson impurity model and its spinless counterpart, 277,278 explored interference effects in the presence of a magnetic field 279 and employed to benchmark simulation methods. 280 Systems simultaneously coupled to both fermionic and bosonic reservoirs were also considered. 281,282 The method was recently generalized to allow for calculating the FCS of thermal transport in the spin-boson model. 283 We note in passing that NRG and DMRG are also able to provide access to density matrices, though for most applications this is more expensive than access to wavefunctions. In particular, the coupling of effective Markovian baths has been shown to be a promising route for obtaining accurate results at long dynamical timescales or steady state. 179,195,284,285

E. Semiclassical approaches
In an effort to describe molecules coupled to anharmonic environments, a semiclassical approach to quantum transport was proposed in Ref. 286. Classical dynamics can be simulated at polynomial complexity in both time and system size, whereas (exact) quantum dynamics always carries an exponential cost. The main idea is therefore to enable simulation of quantum dynamics by mapping a second-quantized Hamiltonian onto a corresponding classical system. The classical model is then solved using standard molecular dynamics techniques.
While the mapping suggested by Swenson et al. is approximate, the approach was shown to correctly reproduce the time dependent characteristics of the resonant level model for elastic 286 and inelastic 287,288 transport. The method was also applied to transport through chaotic cavities. 289,290 Interesting developments include a classical mapping for Hubbard operators, 291 an exact mapping for fermions, 292 and an accurate quasiclassical map that captures both noninteracting fermion dynamics and Coulomb blockade effects. 293 These methods are relatively young in the field and appear very promising, but how they end up fitting into the puzzle remains to be seen.

F. Renormalization group methods
Renormalization group methods can provide highly accurate results under certain considerations, especially when the physics is dominated by modes in a restricted energy range, such as those associated with the Kondo effect. The NRG, mentioned in our earlier discussion of wavefunction methods, is a numerically exact method based on such ideas. Other RG methods employed in transport simulations are the perturbative RG, functional RG and similarity RG techniques. The perturbative renormalization group (RG) method starts from a perturbative expansion up to the leading order in a logarithmic correction. [294][295][296][297] The Keldysh formulation of functional RG (fRG) has also been widely applied to transport problems. [298][299][300][301][302][303] Limits of validity of these powerful approaches is of great interest to the theoretical transport community. One interesting example of such investigations is Ref. 302, where it was found that different results are predicted by the equilibrium and nonequilibrium formulations of fRG.
Somewhat special among RG approaches is the flow equation method or similarity RG, [304][305][306] where the renormalization flow is implemented by a sequence of infinitesimal unitary transformations that leads to a diagonalized Hamiltonian taking all energy scales into account. This is in contrast to usual RG approaches, where the flow leads to the truncation of higher energies; however, the method remains perturbative by nature, and cannot be applied at all parameter regimes. The flow equation method has been successfully applied to the transient and steady state behavior of the nonequilibrium Kondo, 307,308 spin-boson, 309,310 resonant level 311 and Anderson impurity 312 models.
The GQME-based perturbation theory has also been reframed in terms of a real time renormalization group method. 313,314 Real time renormalization group in the framework of kinetic equation was presented in Refs. 315 and 316.

G. Comparative work
Some of the approaches we have discussed are clearly suitable for certain tasks, and analytical considerations can often guide us in our choice. Nevertheless, it is often difficult to understand which method should be used or what reliable benchmarks are available for a particular model and parameter regime, even in surprisingly simple cases. Several studies have taken on the important work of comparing and sometimes contrasting the different approaches described in this section, as well as some of the Green's function approaches described below. 317

III. GREEN'S FUNCTION METHODS IN MOLECULAR ELECTRONICS
We now discuss the application of Green's function methods to transport through molecular junctions. First, we consider the standard diagrammatic technique and its extensions. This is followed by review of numerically exact approaches based on Green's functions.
In quantum many-body theory, the term "Green's function" (GF) actually refers to correlation functions between sets of two or more spatial locations or orbitals, taken at different times. The simplest one-body, or twotime GF takes on the following form: Here,d † m (d m ) creates (annihilates) an electron on molecular orbital m, T c is the contour ordering operator and the τ 1,2 are times on the Keldysh contour. 325 These objects are of interest for three main reasons. First, empirically, few-body (few-time) GFs tend to correspond the kinds of observables that are easier to measure in experiments. For example, particle density, energy density and their fluxes can all be expressed via the simplest GFs, which have only two times. These properties are easily measured in some very complex systems, whereas the many-body wavefunction is generally experimentally inaccessible.
Second, in all but the smallest systems, even if the wavefunction could be measured in complete detail, there is not enough data storage in the world to contain even a rough discretized description of it. The amount of information in a wavefunction increases exponentially with the number of degrees of freedom. On the other hand, GFs contain an amount of information that scales only polynomially with the system size, with two-time GFs scaling quadratically. In the noninteracting case efficiently expressed wavefunctions, density matrices and GFs all provide essentially equivalent information at polynomial scaling with the system size and simulation time. Choosing between them is therefore often a purely aesthetic matter. In the presence of many-body correlations, however, the latter two languages allow us to concentrate on a small subset of relevant degrees of freedom while avoiding any direct reference or access to the exponential Hilbert space of the leads. While density matrix methods remain exponential in the dimension of the region of interest, GFs of any given order are described by correlation functions, the number of which is always polynomial.
Third and finally, GF theory embodies a powerful methodological framework for constructing approximation schemes. Concepts like resummation and conserving approximations make GFs extremely useful. GFs also naturally allow for consistently and seamlessly treating some regions (say, the leads) as noninteracting or analytically solvable, while numerically accounting for detailed structure and interactions in other regions.
The nonequilibrium Green's function (NEGF) technique, which is built around (actually or effectively) noninteracting single-particle orbitals, has long been a tool of choice for ab initio simulations in molecular electronics. The noninteracting NEGF is a convenient starting point for more accurate or general GF methods. This requires identifying a small expansion parameter. Two common choices are related to two generic energy scales that exist in essentially all single molecule junctions (see Fig. 1 for an illustration). The first is the strength of many-body interactions U : for example, intra-molecular electronelectron attraction, or coupling between electronic and vibrational degrees of freedom. The second, Γ, is the coupling strength, or hybridization or escape rate, between the molecule and electronic leads.
When many-body interactions are much smaller than the hybridization, standard NEGF theory exactly solves for the GFs of the noninteracting problem modeled bŷ becomes a small parameter, and perturbation theory in V can be efficiently employed. The approaches we previously mentioned all rely on this idea.
The opposite limit can be addressed by exactly solving the leads and molecule separately: Here, the molecule at the so-called "atomic limit" is assumed to be tractable because it is a physically small region, comprising only a few interacting orbitals. The hybridization,V =V M B , is then treated perturbatively. The pseudoparticle and Hubbard NEGF techniques to be discussed below are based on a hybridization expansion.
In the absence of a small parameter, it may be possible to take advantage of hybrid schemes like the superperturbation theory or dual fermion technique, to be discussed below. The conditions for the applicability of such methods are more involved and not yet fully understood. If approximate schemes cannot be established as accurate, one must rely on substantially more expensive numerically exact schemes, to which we refer at the end of the Section. Below we discuss the various approximations with regard to applications in molecular electronics.

A. Nonequilibrium Green's functions (NEGF)
The NEGF method is an established and widely used technique that is exact and efficiently solvable at the noninteracting limit, but also includes well understood machinery for taking into account weak many-body interactions within diagrammatic perturbation theory. [325][326][327][328][329][330] It relies on the fact that the GF of Eq. (2) obeys a Dyson equation: First, the noninteracting GF, G m1m2 (τ 1 , τ 2 ), is obtained in the absence of many-body interactions. One immediate advantage of the GF methodology is that it is possible to obtain all GF elements in the molecular region M without solving for the dynamics of the full system. The leads can be handled separately, and their GFs are usually solved for in the continuum limit. The effect of the leads then enters as a self-energy term Σ hyb m m (τ , τ ) in the Dyson equation for the molecular GFs, and the effective dimension of the GF becomes the number of orbitals in the molecular region.
Second, interactions are taken into account diagrammatically, also taking the form of a set of many-body or interaction self-energy terms, Σ int m m (τ , τ ). Selfenergies are usually functionals of the GF, such that the solution of the Dyson equation often involves finding the fixed point of a self-consistency relation. These self-energies are by definition zero in all noninteracting regions.
Given the GF of the junction, responses to external perturbations (fluxes, noises and higher cumulants) are obtained from, e.g., the celebrated Jauho-Wingreen-Meir formula 147 and related expressions. 331 In ab initio simulations, an effective noninteracting model of the molecular region is usually constructed from the Kohn-Sham orbitals of a density functional theory (DFT) simulation. DFT, widely utilized for electronic structure simulations, naturally combines with NEGF, because both methods are formulated in the language of quasiparticles or single-particle orbitals. Indeed, the combination of these methods, known as NEGF-DFT, is frequently utilized in the theoretical description of molecular electronics. 190,332 For example, NEGF-DFT was successfully applied to studies of elastic transport in noninteracting systems, 8,15,[333][334][335][336] as well as inelastic transport mostly in the off-resonant tunneling regime. [337][338][339][340][341] We note in passing that approximate NEGF based theoretical schemes for treating resonant IETS were also suggested in the literature. 342,343 The method was further employed in studies of energy transport, 56,57,67,344 for evaluating shot noise, 345 and for studies of disorder in nonequilibirum systems. 346 Besides global responses, it provides access to local characteristics: bond [347][348][349] and current density [350][351][352][353] fluxes, as well as local noise spectroscopy, 354 have all been described within NEGF-DFT.
Notably, the use of the Kohn-Sham Hamiltonian as an effective noninteracting description does not have firm theoretical justifications; nevertheless, it works well in many cases. It has been suggested in the literature that in certain nonequilibrium situations the essentially single-particle nature of this approach may lead to an incomplete description of transport. 355,356 Finally, we remark that while static DFT is most often used, timedependent DFT has also been employed within NEGF calculations. [357][358][359][360] When a molecule is coupled to metallic electrodes to form a junction, its frontier molecular orbitals are renormalized-i.e., shifted and broadened-by the proximity of the substrate. This is due to a combination of delocalization of electrons across the junction and screened polar interactions between electrons in the different regions. These effects are not accounted for in any simulations of the isolated molecule, or by DFT simulations of the full system. To account for such effects, first principle simulations within the GW approximation 361,362 were proposed. [363][364][365] The GW approximation is a selfconsistent GF approximation in principle, but most commonly used implementations employ one iteration of GW to make corrections to the quasiparticle orbitals obtained from DFT. The combination of DFT, GW and NEGF was successfully used in a number of studies. [366][367][368][369][370][371][372][373] Nevertheless, the mean-field nature of GW and the lack of complete screening in junctions may lead to failure of the method when applied to transport simulations. 374 GW formulations remain very expensive. We briefly note here two promising advances that have not yet been used in the context of transport, but would be of interest in reducing the numerical cost: first, range-separated hybrid density functionals can be tuned to provide band gaps on par with the accuracy of GW. 375,376 Second, stochastic formulations of GW allow simulations including thousands of electrons. 377 Significant effort has been devoted to developing NEGF-based theoretical methods for understanding optical measurements in current-carrying molecular junctions. 378 For example, NEGF methods for linear optical response, 379 current-induced light, 380 Raman spectroscopy [381][382][383][384] and multi-dimensional spectroscopy 218 have been formulated. We note that traditional approaches to nonlinear optical spectroscopy, 230 which employ bare perturbation theory, do not satisfy conditions for building conserving approximations in the NEGF diagrammatic approach. 237 Recently, a fluxconserving version of nonlinear optical spectroscopy for application in open molecular systems was suggested. 238 Note also that within diagrammatic expansions, one can only treat relatively weak light-matter interactions, when both light and matter degrees of freedom are treated quantum mechanically. On the other hand, there is no such restriction on the strength of interaction with a classical radiation field. In the latter case, the Maxwell-NEGF method was utilized for simulation of junctions driven by external classical fields. 385,386 NEGF has also often been applied to molecular spintronics. A few examples are studies of magnetic field effects on electronic conduction, 387 spin pumps, [388][389][390] spin-flip IETS, 391-394 transport in helical molecules (e.g., DNA), [395][396][397][398][399][400][401] quantum interference in spintronic devices, 402 thermospintronics. [403][404][405] Green's function methods recently also were employed in studies of formation and dynamics of skyrmions. [406][407][408][409] NEGF has further been applied to study shot noise. For example, shot noise noise spectroscopy of spin current 410,411 and within inelastic transport [412][413][414][415][416][417][418] have been discussed in the literature. The theory of light emission as a probe of electron shot-noise was put forward, 419 and current fluctuations in the transient regime 420,421 have been discussed. The notion of local noise spectroscopy (LNS) (contrary to total junction noise) has also been introduced. 354 This concept is useful in theoretical consideration of, e.g, local molecular light emission patterns, heating maps, or mapping out local interaction within the junction (see Fig. 2). While in single molecule junctions the technique has purely theoretical applications, in extended systems (such as graphene nanoribbons) LNS predictions are measurable. In addition to noise, higher moments and FCS have been considered. As in flux modeling, study of FCS and waiting time distribution (WTD) for systems with strong manybody interactions is challenging within NEGF. [422][423][424][425][426] At the same time, rigorous treatment of FCS within NEGF for noninteracting systems is well established in both the steady-state 217,427 and transient [428][429][430] regimes.
In summary, NEGF has many important advantages when applied to simulation of responses in single molecule junctions. It has a well-established route to account for interactions within controlled diagrammatic perturbation theory. The conserving character of the resulting approximations 235 functions for the self energies. 325,431 In the case of bilinear molecule-contacts coupling, the embedding or hybridization can be exactly taken into account (i.e. the corresponding term in the self-energy is exact). NEGF-DFT allows for convenient pairing between ab initio electronic structure calculations and the evaluation of transport properties. FCS can be evaluated exactly for noninteracting and approximately for interacting systems, allowing for modeling of noise measurements in junctions. Nevertheless, within the NEGF diagrammatic formulation it is challenging to treat strong many-body interactions. [432][433][434] Also, information regarding molecular many-body states (as opposed to single-body orbitals) is more difficult to access, as it requires higher-order GFs. Strong many-body interactions and state-specific information are of value for, e.g., the characterization of optoelectronic devices, where in order to measure optical response the molecule should be relatively weakly coupled to contacts. A more convenient description of this regime is given by diagrammatic expansions formulated directly in terms of the interacting molecular many-body states. Below we discuss GF techniques conveniently applicable in this regime.

B. Pseudoparticle NEGF
Pseudoparticle (or auxiliary-operator) NEGF was first utilized in simulations of transport within the Anderson impurity model in Refs. 435-438. More recently, such techniques have become useful in the context of driven materials, as solvers for the effective impurity models appearing in nonequilibrium dynamical mean field theory. 439 We note in passing that the dynamical mean field theory, originally formulated for bulk materials, can also be employed in the treatment of molecular electronics including large interacting regions. 440,441 The relationship between PP-NEGF and standard NEGF is established via a spectral decomposition of electron creation (annihilation) operatorsd † m (d m ) in the many-body state basis {|S } of the isolated molecular region M . Here, it is assumed to be a complete local basis:d † m = S1,S2

PP-NEGF introduces second quantization in an extended
Hilbert space, where pseudoparticle operatorsp † S (p S ) create (destroy) many-body state |S starting from an unphysical vacuum state |vac Similarly to the NEGF, Eq. (2), the central object of interest is the single pseudoparticle GF The formulation of diagrammatic technique in the PP-NEGF follows closely that of standard NEGF, with the difference that the NEGF expansion in weak many-body interaction is substituted by a PP-NEGF expansion in weak system-bath coupling. The PP-NEGF therefore expands around the nonequilibrium atomic limit. A more important difference is that the PP-NEGF is formulated in an extended Hilbert space, the physical subspace of which is defined by the normalization condition Sp † Sp S = 1.
The necessity to impose the normalization condition on unrestricted solution is the weakest point in the PP-NEGF methodology. Contrary to standard NEGF, where GFs are correlation functions between excitation operators, PP-NEGF considers correlations between states themselves. In this sense it is closer in spirit to QME formulations, but allows for a more natural treatment of non-Markovian effects. PP-NEGF was shown to be useful in studying transport through strongly correlated systems, 442,443 in modeling inelastic transport in the resonant tunneling regime with strong electron-vibration interactions or in presence of anharmonicities, 444 in treating electronic and energy transport on the same footing, 445,446 and in first principles simulations of Raman spectroscopy in single molecule junctions 447 (see Fig. 3).
The ability to access state-specific information also makes PP-NEGF useful in understanding nonadiabatic molecular dynamics in the regime of weak but nonvanishing molecule-contact electron exchange coupling. In particular, in Ref. 448 PP-NEGF was utilized to derive equations for the nuclear dynamics that reproduce the surface-hopping formulation in the limit of small coupling, and Ehrenfest dynamics when information on different charging states of the molecule is traced out. We note that closely related to PP-NEGF methods are those where a reduced molecular propagator between contour times, Tr B ρ B e −i z 2 z 1Ĥ (z) , replaces the GF as the fundamental object. Expansions of this type appear as initial approximations in numerically exact schemes. 255,[449][450][451] Propagator formulations trade some of the useful machinery of GF theory for the ability to treat intra-molecular interactions exactly while expanding in the molecule-bath hybridization, without the need for artificially extending the Hilbert space as in PP-NEGF. Propagator methods do enable a type of modified diagrammatic resummation, and can be used to construct conserving approximations. For instance, they have been used to explore transport in junctions with both electron-electron and electron-phonon interactions within the Anderson-Holstein model, 452 as well as in a study (with comparisons to exact results) of the splitting of the Kondo conductance peak resulting from a nonequilibrium bias voltage. 453 To summarize, PP-NEGF enjoys all the standard tools of quantum field theory. This includes well-developed rules for the diagrammatic technique together with established ways of building conserving approximations. Like reduced density matrix techniques, PP-NEGF al-lows access to state-specific information while presenting a simple and controlled approach to account for timenonlocal (non-Markovian) effects. Unfortunately, formulation in the extended Hilbert space and the necessity to impose the normalization condition Eq. (7) makes the method more difficult to reliably apply to full counting statistics (FCS). Only noise (second cumulant) estimates within perturbative expressions have so far been reported the literature. 454,455 Practical applications using PP-NEGF are often limited to the lowest (second) order in the diagrammatic expansion-the non-crossing approximation-which may fail qualitatively and break symmetries even in rather simple limits (see Fig. 4). We now turn to discuss another many-body flavor of the NEGF with several advantages.

C. Hubbard NEGF
Historically, Hubbard Green's functions were introduced for the description of electronic correlations in narrow bands. 456 Their main usage was in the formulation of perturbation theory around the atomic limit. Hubbard NEGFs are a nonequilibrium analog of the formulation, where the isolated molecule is the starting point and (similarly to PP-NEGF) a diagrammatic perturbative expansion is considered in the molecule-contacts coupling. The main idea is to consider correlation functions between Hubbard operators: These take the form Similarly to PP-NEGF, Hubbard GFs are constructed around the nonequilibrium atomic limit and yields molecular state-specific information. The spectral decomposition (4) allows for a standard GF to be calculated from the Hubbard GF, but the inverse process is not possible. Contrary to PP-NEGFs, Hubbard GFs live within a physical Hilbert space, allowing for example simulation of FCS. While standard NEGF contains only information about correlations between transitions and PP-NEGF contains that of states, Hubbard NEGFs include both. Some of the first applications of Hubbard NEGFs to transport were carried out within the equation-ofmotion (EOM) approach, where auxiliary field techniques allow for closed-form equations in terms of the GFs and their functional derivatives to be written. 457,458 The method was useful in modeling elastic [459][460][461] and inelastic 462 transport, and was also utilized in ab initio simulations. 463 Despite these successes, the somewhat arbitrary nature of the auxiliary fields makes it difficult to generalize these methods into systematic expansions.
Diagrammatic technique for Hubbard GFs was developed in studies of strongly correlated equilibrium lattice models. 464,465 The technique is based on an observation that (anti-)commutator of two Hubbard operators yields not a number, but another Hubbard operator; and on the equilibrium form of the density operator. Given these, it can be shown that any Hubbard operator within an n-time correlation function can be permuted with other Hubbard operators and with the density operator under the trace, leading to an expression containing only n − 1time correlators. Iterating this procedure generates an expression analogous to Wick's theorem and allows for diagrammatic expansions to be formulated: by resumming diagrams, one then arrives at a Dyson-like equation for so-called locators. Finally, locators be multiplied by time-non-local factors accounting for the spectral weight and renormalization of quasiparticles to obtain the Hubbard GFs.
In Ref. 466 this idea was extended to Keldysh GFs. The generalization assumes (as in NEGF) that in some far-off past time the molecule and contacts were separated and that after transients die out, the system does not remember its initial state. Thermal equilibrium for the separated molecule and contacts are therefore assumed for the far-off initial state. With this, each term in the Taylor expansion of the S-matrix can be evaluated by using the standard Wick's theorem for bath operators, and the generalized Wick's theorem for Hubbard operators within the molecule. Later, in Ref. 251, it was pointed out that the generalized Wick's theorem can equivalently be formulated for any local density operator that is a function of the molecular Hamiltonian only. For example, such a density operator results from the Markov Redfield/Lindblad QME, and therefore the Hubbard NEGF technique can be formulated as an expansion starting from an unbroadened QME solution.
Hubbard GF diagrammatics can be applied to FCS. 467 It was shown that it is an inexpensive method capable of reproducing satisfactory noise characteristics in molecular junctions over a wide range of parameters. In this sense, Hubbard NEGFs perform comparably to standard NEGF at the weakly interacting limit and to QME for strong many-body interactions. In the Coulomb blockade regime the second order Hubbard NEGF theory is able to reproduce experimentally observed shot noise suppression, whereas NEGF and QME treatments at the same perturbation order fail (see Fig. 5).
Another recent application of Hubbard NEGFs is in the derivation of current-induced nuclear forces: in particular, the nontrivial electronic friction force was evaluated. 468,469 Ref. 468 derives a general expression for the friction, applicable to any form and/or strength of electron-nuclei coupling. The Hubbard NEGF result was benchmarked with respect to exact simulation in a simple model, and performed well. The expression obtained appears to be defined by retarded projection of the Hubbard GFs, and contains previous known results as its limiting cases. For example, the celebrated Head-Gordon and Tully (HGT) expression for electronic friction is shown to be given by a (diagonal) subset of the contributions given by Hubbard GF theory, taken at the noninteracting quasi-particle limit. The Hubbard GF accurately reproduces the HGT result in the limit of weak molecule-contact coupling. Based on this developmemt, Ref. 469 discussed the possibility of engineering molecular friction in single molecule nanocavity junctions.
The Hubbard NEGF theory has proven particularly useful in simulations of molecular optoelectronic devices. 378 Here, an insulating layer is placed between the molecule and contacts in order to prevent nonradiative energy transfer between them. The corresponding model parameters are therefore close to the atomic limit. Formulation in the many-body state description is helpful  for avoiding spurious orbital-shifting and renormalization effects upon charging and/or excitation. This picture also lends itself to a simple explanation of the difference between transport and optical gaps in junctions, where the two gaps correspond to transitions between different pairs of molecular many-body states.
Recently, the Hubbard NEGF was employed in ab initio simulations of transport and optical response in molecular junctions. Ref. 41 modeled an experiment dealing with photocurrents induced in nitroazobenzene molecular junctions, and Refs. 92 and 470 respectively, simulated bias-induced phosphorescence and luminescence (see Fig. 6).
In summary, much like PP-NEGF, the Hubbard NEGF is formulated in the many-body state basis of the isolated molecule. This makes it useful for treating systems with strong many-body interactions within the molecule and weakly coupled to the baths, such as molecular optoelectronic devices. The formulation can exactly account for all intra-molecular many-body interactions, and provides a way of incorporating high-level quantum chemistry methods (CI, CCSD, etc.) into the realm of open nonequilibrium molecular junctions. Recently developed diagrammatic technique for Hubbard NEGFs allows one to take into account system-bath interactions in a controlled manner. Simultaneously, access to correlation functions of excitations makes Hubbard NEGF similar to the standard NEGF, which may account for the surprisingly high level of accuracy observed at strong systembath coupling. Also, contrary to the PP-NEGF, the Hubbard NGEF is constructed within the physical Hilbert space, making it usable for evaluating FCS. The main difficulty with the Hubbard NEGF is the absence of a clear way to construct an analog of the Luttinger-Ward functional. Because of this, so far no clear rules for constructing conserving approximations within the Hubbard NEGF have been found.

D. Dual fermion approach
We have discussed NEGF theory, which is formulated around the noninteracting or quasiparticle limit, and therefore applicable to weak many-body interactions, U Γ. We have also discussed PP and Hubbard NEGF; these are formulated around the atomic limit, and suitable for weak system-baths coupling, U Γ. The parameter regime where the two characteristic energy scales are similar, U ∼ Γ, has to be treated within an approach which does not rely on existence of small parameter (see Fig. 1). Rigorous treatment of such regimes therefore requires numerically exact approaches. However, such approaches are typically computationally expensive and limited to very simple, minimal models like the Anderson impurity model. Thus, especially for ab initio purposes, inexpensive approximate methods that remain reasonably accurate in the absence of a small parameter are in high demand.
A new method that is potentially applicable to realistic simulations in the absence of a small parameter is the dual fermion (DF) approach. Historically, the DF method was formulated for equilibrium lattice models as a way to account for nonlocal corrections to dynamical mean field theory. [471][472][473][474] However, it was quickly realized that it can serve as a generic "superperturbative" method for solving quantum impurity problems, both in 475 and out 476,477 of equilibrium.
The formulation of the DF method for solving an impurity problem or molecular junction comprises two main parts. In the first step, a (solvable) reference system containing the many-body interactions is chosen. The physical system's action, S[d, d], is written in terms of the known action of the reference systemS[d, d] and the difference between the exact hybridization self energy Σ B , and that of the reference systemΣ B : Here, the subscript indices 1, 2 indicate both quantum numbers and contour times, and a sum over repeated subscripts is implied; G 0 is the noninteracting GF of the molecular region M ; and S int [d, d] is the action resulting from the many-body interactions. The reference action is interacting, and therefore Wick's theorem does not hold and the standard tools of many-body perturbation theory cannot be directly applied to it. To resolve this, in the second step, a Hubbard-Stratonovich transformation is employed. This introduces a new auxiliary quasiparticle, the "dual fermion", and makes it possible to formulate a perturbation theory in the coupling between the dual and real fermions. When the real fermions are traced out, the free dual fermion theory is an exact effective theory reproducing the interacting reference action at its zeroth order. Corrections to this take the form of effective many-body interactions between the dual fermions, which can now be treated perturbatively. In Ref. 476, which considered the nonequilibrium Anderson impurity model (see Fig. 7a), the reference system was chosen to be the interacting impurity site along with one to three bath sites. This necessarily entails an approximate treatment of the impurity-bath hybridization. The finite, closed reference system used in Ref. 476 (see Fig. 7b) does not account for the dissipative dynamics of the open quantum system, or for the nonequilibrium population dynamics, even at the classical level. It therefore makes it difficult to capture the steady state. Ref. 478 suggested that an auxiliary open system with Markovian leads, described by a Lindblad/Redfield QME, be used as the reference system instead (see Fig. 7c). This aux-DF method allows for the problem to be formulated directly in the nonequilibrium steady state and offers substan-tially improved accuracy in comparison to the closed system DF method. It also significantly reduces the numerical cost of simulations by avoiding the long time propagation that is required in the original formulation.
The idea of constructing accurate QMEs by explicitly taking into account the dynamics of some bath sites is known in the literature as the auxiliary master equation approach (AMEA), and has been employed as an impurity solver in nonequilibrium dynamical mean field theory. 195,[479][480][481] This method becomes exact at the limit of an infinite number of auxiliary sites, 482,483 but remains approximate in practice, where only finite number of auxiliary sites can be treated. The accuracy of AMEA depends on a fitting procedure for the effective hybridization function Σ of the original problem in the reference system. A large number of auxiliary modes is often needed to obtain a high level of precision. There is therefore a trade-off between accuracy and efficiency: more auxiliary modes enable a better fit, but make the solution of the auxiliary QME more expensive. In comparison, the alternative and complementary strategy embodied by the aux-DF method is to perform a perturbative expansion in the difference between the hybridization functions of the physical and reference systems. This can be done for a smaller number of modes, while relying on relatively poor fitting, because the superperturbation theory provides further improvement over the bare AMEA result. While more expensive than AMEA for a given number of auxiliary sites, aux-DF therefore also yields substantially more accurate results at that same limit (compare solid and dotted lines in Fig. 7d).
Test simulations of the Anderson impurity model performed in Ref. 478 were benchmarked against numerically exact tDMRG and Inchworm continuous time quantum Monte Carlo (CT-QMC) results and demonstrated high accuracy (see Fig. 7d). The Inchworm method will be discussed below. Note that aux-DF requires only a fraction of the computational cost of the numerically exact techniques to achieve results that (for the parameters tested) are of comparable accuracy. Like the DF method in equilibrium, aux-DF correctly describes the limits of weak and strong molecule-bath coupling, and allows for interpolation between them. The superperturbation theory converges rapidly at both limits because it becomes exact when either the hybridization or the many-body interactions are small. Indeed, for weak interactions or strong molecule-bath coupling, the vertices can be neglected and DF reproduces standard NEGF perturbation theory. For strong interactions (or the atomic limit, where molecule-bath coupling is weak) the DF method is also accurate, due to exactly accounting for the local many-body interactions in the reference system. Currently, the main obstacle on the road to an ab initio implementation of the aux-DF methodology is in the need to find a numerically inexpensive way to evaluate vertices in the solution of the reference system.
We also note that very recently, a generalization of the method able to address not just electron, but also energy transport was presented and applied to the Anderson-Holstein model: the auxiliary QME-dual boson technique (aux-DB). 484 This allows to account for electronelectron, electron-phonon or electron-photon interactions while treating particle and energy fluxes on an equal footing.

E. Time-dependent processes with Green's functions
Molecular electronics is usually focused on steady state properties, which tend to be more experimentally accessible. Nevertheless, measuring and modeling time-dependent and transient phenomena is an important challenge when considering time-dependent quantum transport, time-resolved photoabsorption, pumpprobe type experiments, Auger decay processes and other ultrafast phenomena. 421,[485][486][487][488][489] As such problems typically require GFs explicitly depending on two times, they tend to be more computationally challenging in GF methods: in time-local methods like wavefunction and density matrix approaches, the state is always described by a single time index. For noninteracting systems, some of the most efficient methods for obtaining GFs actually rely on wavefunction dynamics. 490,491 The upside is that GFs also provide rich time-nonlocal information.
One exception to the above rule is problems with periodic driving, which can be relatively easily treated in GFs using Floquet theory 399,492-499 and similar considerations. 383,500,501 If two-time GFs do need to be considered, a straightforward approach is to employ discretization on a two-dimensional time grid. Accounting for the fast decay of GFs with separation of their time variables then allows some simplifications. Such considerations have been applied to the equations of motion of GFs in both their integral (Dyson) or differential (Kadanoff-Baym) form. [502][503][504][505] Other works perform the time discretization directly on the Keldysh contour. [506][507][508][509] To perform long-time simulations at better computational scaling, approximations can be used. For example, the generalized Kadanoff-Baym ansatz (GKBA) and other semianalytical ideas have been used. [510][511][512][513][514][515][516][517] These effectively reduce the computation to single-time evolution. For noninteracting NEGF in particular, this can be done in several highly efficient ways without the need for approximations. 490,518-521 A particularly interesting idea that should be noted in this context was proposed in Ref. 522, where the problem of solving the Dyson equation is mapped onto a noninteracting auxiliary problem with additional degrees of freedom.

F. Quantum Monte Carlo methods
The past decade has seen a great deal of progress in numerically exact GF methods for simple models of single molecule junctions, in particular with regard to continuous time quantum Monte Carlo methods (CT-QMC). All such approaches presently remain far too expensive to be usable in an ab initio context, where a large interacting basis must be taken into account and interactions exist throughout both the molecular and lead regions. Nevertheless, in simple models, they allow for reliable access to genuinely strongly correlated physics, where there is no small parameter, and where interpolation between analytically solvable regimes is insufficient.
The CT-QMC approaches we discuss below have equilibrium predecessors formulated in imaginary (rather than real or Keldysh) time. 523 Real time Quantum Monte Carlo approaches to even the most minimal models used to describe nonequilibrium electronic transport through junctions were commonly thought to be limited to short propagation times by the dynamical sign problem: an exponential decrease in the signal-to-noise ratio as a function of time. For the spin-boson model, this view has been challenged by path integral Monte Carlo methods formulated in terms of density matrices and influence functionals. 524,525 Two sets of rather different CT-QMC methods have addressed fermionic transport, and will be briefly reviewed below. One relies on an expansion in many-body interactions and the other on one in molecule-lead hybridization, but both include ideas that are independent of the choice of expansion.

Interaction expansion
The real time interaction expansion for nonequilibrium impurity models was introduced in Ref. 526, and is based on an earlier equilibrium method. 527 Confusingly, it is sometimes referred to as the "weak-coupling" approach. The term "weak-coupling" here refers to a physical regime of the impurity problem at weak interaction strength, not to the strength of the molecule-lead coupling or hybridization. However, the interaction expansion is (a) numerically exact, and therefore not limited to a particular regime; and (b) is in fact more efficient in the presence of weak interactions and strong molecule-lead coupling. Despite featuring a dynamical sign problem, the method was successfully used to explore the transient and steady state current-voltage characteristics of a junction with electron-electron interactions 526,528 and further employed within dynamical mean field theory. 497 More recent advances in real time interaction expansion Monte Carlo stem from an interesting realization about starting from the CT-QMC interaction expansion on the Keldysh contour of Ref. 526, but then summing explicitly over Keldysh branch indices to obtain a moment expansion. 529,530 It turns out that this leads to two crucial simplifications: first, the approximate dynamics simulated by the Monte Carlo process become exactly unitary; and second, the dynamical sign problem plaguing any naive interaction expansion essentially vanishes. The cost of the summation over branch indices is exponential in the order, such that the expansion must be truncated at a finite but high order in the many-body interaction (so far, orders around 10-15 have been used). An altered scheme allows for removing this explicit exponential cost at the cost of reintroducing a sign problem, and may lead to a set of improved intermediate techniques. 531 The bare interaction expansion (regardless of how efficiently it is evaluated) has a very limited convergence radius. However, it was shown that this radius can be extended with the aid of analytical continuation techniques. 529 In later work, it was demonstrated that understanding the analytical structure of the impurity problem can be leveraged to construct conformal transformations that provide highly accurate analytical continuation procedures that continue to work deep within the strongly correlated regime. 532 Additionally, statistical noise from the Monte Carlo procedure is larger at large interaction strengths, but can be reduced by postselecting samples that obey known, exact physical constraints at the limit of infinite interactions. 532

Hybridization expansion
The use of the molecule-lead hybridization expansion in quantum transport began a year earlier than that of the interaction expansion: it was introduced in Ref. 533, which considered electron-phonon interactions. This was based not on CT-QMC, but on an earlier algorithm by Hirsch and Fye. 534 Here, times along the contour were discretized, and Monte Carlo steps comprised changing the occupation at one discrete interval. This was rapidly followed by a continuous time approach, which eliminates the need for discretization; 535,536 and by applications to electron-electron interactions. 526,537,538 An interesting recent development shows that the CT-QMC problem for a model with two spin channels can be mapped onto one with a single channel, resulting in a significant reduction in the sign problem. 539 At any finite time, the bare hybridization expansion converges at all interaction strengths and no analytical continuation is required. Nevertheless, the method is strongly limited by the dynamical sign problem to times that are on the order of the inverse molecule-bath coupling.
The hybridization expansion Monte Carlo methods are formulated in terms of reduced propagators, much like the propagator techniques mentioned in Sec. III B. As illustrated in Fig. 8a, the expansion for the propagator (thick horizontal line, with contour time going from left to right) can be expressed in diagrammatic notation as a sum over products of isolated molecular propagators (thin horizontal line) and hybridization lines given by molecule-lead GFs (curved wiggly lines). GFs of any order, as well as density matrix elements in either the molecule or leads, can be expressed in terms of modified propagators with additional operators inserted at the appropriate Keldysh times.
It is possible to perform partial resummations over a contour "self-energy"; for example, at lowest order this provides the noncrossing approximation for the propagator, as illustrated in the top half of Fig. 8b and denoted by a "bold" or partially thickened gray line. Bold-line propagators of low orders can be obtained semi-analytically, and then used to construct a dressed expansion containing only corrections to the underlying theory. These corrections can be summed in a Monte Carlo procedure expressed entirely in terms of the renormalized bold-line propagators (lower half of Fig. 8b), as first shown in equilibrium. 449 When applied to real time population dynamics, this idea was shown to greatly reduce the dynamical sign problem, allowing access to timescales several times longer than those accessible by the corresponding bare expansion. 255,450 The bold-line method was then extended to nonequilibrium GFs and used to explore the spectral and current-voltage characteristics of junctions with electron-electron interactions. 451,541,542

Inchworm expansion
The next advance in the field, dubbed "Inchworm Monte Carlo", is not specific to the hybridization expansion; however, this was its first application. 540 The main idea is first to take advantage of the fact that propagators over short contour time intervals are easily evaluated; and second to note that propagators obey a kind of contour causality, in the sense that propagators over long time intervals can be efficiently expressed in terms of propagators over shorter time intervals (see Fig. 8c). A formally exact series for the exact propagator over some interval (t, t ) can be written in terms of all propagators over intervals (s, s ) with t ≤ s, s ≤ t ↑ ≤ t . At an intermediate stage of the algorithm, the Inchworm time t ↑ , denoted by the arrow in Fig. 8c, is the time up to which all propagators are known. By performing the summation over all diagrams contributing to the expansion, propagators up to the longer time t are obtained. A series of such "inching" steps eventually generates the propagator over the entire contour and physical observables, butcrucially-propagators over long intervals are only calculated when propagators over slightly shorter intervals are already known and taken advantage of. Therefore, each inching step has a relatively small sign problem, which in practice does not grow exponentially with time. 540 As an example (see Fig. 8d), a set of trivial diagrams that are generated by low order self energies can easily be shown to generate an exponential dynamical sign problem in both the bare (left) and bold-line (right) expansions. 540 Diagrams of this form are summed implicitly within the Inchworm expansion, such that this particular source of sign problems is completely removed.
The Inchworm hybridization expansion can be truncated in diagram order for improved efficiency, though this is not necessary in all regimes: an Inchworm expansion of order n converges to a dressed expansion for the propagator self energy up to that order, and often converges very rapidly. No analytical continuation is necessary. When higher order contributions (number of vertex pairs n 10, corresponding to order 20 self-consistent perturbation theory) are needed, fast summation techniques become useful. 543 Nevertheless, the method is based on a resummed perturbation theory. If the order needed for convergence at a given set of parameters is too high, it eventually fails. Furthermore, the computational cost of existing implementations is effectively quadratic in the propagation time, whereas the interaction expansion of Ref. 529 can provide results at very long times with no additional expense.
Whereas Ref. 540 only considered population dynamics in interacting quantum junctions, the method was quickly extended to currents and GFs, 544 and used to explore the transient and steady state properties of voltagedriven interacting junctions in the strongly correlated regime. 453 The Inchworm method also proved to be suitable to the numerically exact evaluation of FCS. 545 This was used to explore the effect of junction geometry on electronic current and its fluctuations, 546 and provides access to thermal transport properties and entropy generation. 547 We briefly note that the Inchworm idea is not limited to the real time hybridization expansion in quantum transport: it has also been applied within dynamical mean field theory; 548 shown to be effective for combating the dynamics sign problem in the spin-boson model using two different diagrammatic expansions; 549-551 and very recently, to be applicable to the imaginary time sign problem in multiorbital impurity models. 552

IV. CONCLUSIONS
We presented a short pedagogical review of methods utilized in the study of quantum transport through molecular junctions. Three broad categories were considered: wavefunction based methods, density matrix and Green's function (GF) formulations. While the approaches are equivalent when treated exactly, in practice exact treatment is only possible for noninteracting systems where the single particle picture is accurate. For interacting problems, each category of methods leads to different approximations and numerically exact schemes, with widely varying advantages and disadvantages. In wavefunction methods, observables and dynamics can be expressed by quantities taken at a single time, but the amount of information scales exponentially with the size (e.g. number of orbitals) of the complete system. This includes both the molecule and contacts. In density matrix methods, at the cost of introducing either non-Markovian dynamics of single-time observables or a Markovian approximation, the spatial region contributing to the exponential scaling is reduced to that of the interacting molecular region. In GF methods, at the cost of having to address time-nonlocal dynamics of two-or multitime correlation functions, the exponential spatial scaling can be completely removed. Methods from all categories have seen rapid progress in recent years, and modern approaches can efficiently mix and match between categories to best suit specialized problems. Theoretical considerations were illustrated with numerical examples in model and ab initio simulations, and promising avenues for future research were pointed out.
The main focus of this review is on GF approaches. GFs often provide an efficient and accurate way to simulate experimentally measurable observables. This includes electronic and energy currents carried by either fermionic charge carriers or bosonic modes. It can also include current fluctuations and full counting statistics. On the theoretical side, GF methods provide a methodological framework for constructing controlled approximation schemes around solvable limits. We discussed and compared standard NEGF, which expands around the noninteracting limit; its many-body flavors (PP-and Hubbard NEGF), formulated around the atomic limit; and superperturbation theory (dual-fermion and dual-boson) that can provide accurate predictions in both limits, and bridges between them.
We also discussed numerically exact GF schemes based on Quantum Monte Carlo algorithms, that can provide reliable results at all parameter regimes, though only for very simple models and at a large computational cost. Until recently it was thought that such methods are limited to very short propagation times in nonequilibrium systems, due to dynamical sign problems. However, new methodological frameworks introduced in the last few years, such as Inchworm Monte Carlo, have shown that this sign problem can often be circumvented. Such diagrammatic Monte Carlo methods are based on perturbative expansions similar or identical to the ones used in approximate GF techniques. However, the Monte Carlo algorithms allow for systematically summing all corrections in an unbiased manner up to very high, or some-times effectively infinite, orders.
In the context of molecular electronics, we noted that diagrammatic perturbation theory around the noninteracting NEGF is convenient for treating molecular junctions with relatively weak many-body interactions, such as intra-molecular electron-electron interactions. This is often the case when considering, e.g., transport through π-conjugated organic polymers. On the other hand, the PP-and Hubbard NEGF diagrammatic expansions are perturbative in the molecule-bath coupling strength rather than the interaction. This makes them useful, e.g, in the presence of strong electron-electron interactions within the molecule, as is typical in molecular optoelectronic devices. Dual techniques, a newer alternative, are relatively inexpensive universal impurity solvers that can remain accurate even in regimes where neither simple perturbative expansion works well. Each of these methods has, or could potentially have, a numerically exact technique associated with it that can tackle nonperturbative regimes.
The theory of nonequilibrium transport through single molecule junctions is at a turning point. The community has developed an impressive array of numerical and analytical techniques that can address a wide variety of physical problems, and progress is still being made. We hope this review will be helpful to readers in navigating through the forest of available methods and using them to push our scientific understanding and technological mastery of these systems to new heights. Data sharing is not applicable to this article as no new data were created or analyzed in this study.