Towards Noise Simulation in Interacting Nonequilibrium Systems Strongly Coupled to Baths

Progress in experimental techniques at nanoscale makes measurements of noise in molecular junctions possible. These data are important source of information not accessible through average flux measurements. The emergence of optoelectronics, the recently shown possibility of strong light-matter couplings, and developments in the field of quantum thermodynamics are making measurements of transport statistics even more important. Theoretical methods for noise evaluation in first principles simulations can be roughly divided into approaches for weak intra-system interactions, and those treating strong interactions for systems weakly coupled to baths. We argue that due to structure of its diagrammatic expansion, and the use of many-body states as a basis of its formulation, the recently introduced nonequilibrium diagrammatic technique for Hubbard Green functions is a relatively inexpensive method suitable for evaluation of noise characteristics in first principles simulations over a wide range of parameters. We illustrate viability of the approach by simulations of noise and noise spectrum within generic models for non-, weakly and strongly interacting systems. Results of the simulations are compared to exact data (where available) and to simulations performed within approaches best suited for each of the three parameter regimes.

difficulty treating strong system-bath couplings. In this respect, an important extension of the latter approach is formulation of the real time perturbation theory (and similar approaches) [44][45][46][47][48][49][50] , which allowed evaluation of shot noise accounting for system-bath coupling within a systematic perturbation expansion 51 .
In this study, we use the recently proposed nonequilibrium diagrammatic technique for Hubbard Green functions (Hubbard NEGF) 52 to simulate FCS and noise spectrum in molecular junctions. We note in passing that earlier studies of Hubbard NEGF 53,54 were restricted to the evaluation of two-time correlation functions only, while simulation of noise spectrum requires evaluation of multi-time correlation functions. The latter is only possible within diagrammatic consideration (for a detailed comparison between earlier studies and diagrammatic approach to Hubbard NEGF, see ref. 52). Contrary to the standard NEGF, which considers multi-time correlation functions of elementary (de-)excitation operators ˆ † d i (d i ) -where i is a single electron orbital, the Hubbard NEGF deals with correlation functions of projection operators = X S S S S 1 2 1 2 -where |S 1,2 〉 are many-body states of the system. Similarities between the two approaches (illustrated by spectral decomposition = ∑ˆd S d S X i S S i SS , 1 2 1 2 1 2 ) indicate that the structure of the diagrammatic technique for the Hubbard NEGF should (to some extent) resemble that of the standard NEGF (although expansions themselves are different: Hubbard NEGF considers perturbation series in system-bath coupling, while NEGF expands in small intra-system interactions). A similar way of treating self-energies in the two techniques makes expectation that Hubbard NEGF will be successful in accounting for relatively strong system-bath couplings plausible. At the same time, formulation in the language of many-body states makes the Hubbard NEGF similar to the QME approaches. Thus one may expect that accounting for strong intra-system interactions is feasible within the Hubbard NEGF. Note that similar to the real time perturbation theory of refs 44-46, 51 the Hubbard NEGF expansion is systematic. Note also that lowest order Hubbard NEGF diagrams are still simple enough to be applicable in first principles simulations. This paves a way to ab initio calculations of noise beyond the weak interaction limit of ref. 42.
Below, we present noise simulation in model systems within the Hubbard NEGF, comparing results to those of other approaches. The presentation is followed by a discussion of the approach indicating its strengths and weaknesses, and pointing to directions for further research. Finally, we briefly introduce the Hubbard NEGF diagrammatic technique; details of the technique and its application to noise simulation can be found in ref. 52 and in the Supplementary Information.

Results
We consider Anderson impurity as a junction model. Its Hamiltonian is The first line utilizes second quantization, while the second employs many-body states. Here, σ † d and σ † c k are creation operators for electron with spin σ in the system or contact state k, respectively; = σ σ σˆ † n d d and = σ σ σˆ † n c c k k k , |S〉 = |0〉, |a〉, |b〉, |2〉 are empty, spin-up, spin-down, and double-occupied many-body states with energies E 0 = 0, E a = ε ↑ , E b = ε ↓ , E 2 = ε ↑ + ε ↓ + U; m = 0a, b2, 0b, a2 are single electron transitions between pairs of many-body states (0a and b2 are spin up transitions, 0b and a2 are spin down), and L and R are left and right contacts. Within the model, U represents intra-system interaction, system-baths couplings are characterized by escape rate matrices, (K = L, R), which are energy-independent in the wide band approximation.
We simulate the current and noise spectrum at the left interface for a steady-state transport situation Green function techniques use exact (Meir-Wingreen) expression for the current, and approximate (perturbative in system-baths couplings) derivation similar to that of ref. 30 for the noise spectrum. In the Hubbard NEGF, this derivation leads to diagrams presented in Fig. 1d. FCS is introduced as usual by dressing transfer matrix elements σ σ V k 2 ) with a contour-branch-dependent counting field λ for k ∈ L. Zero frequency noise (second cumulant of the FCS) is obtained from a derivative of the dressed current, (see, e.g., ref. 55 for details). Simulations are performed within the Hubbard NEGF, standard NEGF, and Lindblad/ Redfield QME for non-interacting (U = 0), weakly interacting (U ≤ Γ 0 ), and strongly interacting ( Γ  U 0 ) cases (system-bath coupling strength Γ 0 is used as a unit of energy). Note that Lindblad/Redfield QME does not yield noise spectrum; current and zero-frequency noise were simulated within the FCS (see ref. 56 for details). Results of the simulations are presented in units of I 0 = eΓ 0 /ħ for current and S 0 = e 2 Γ 0 /2ħ for noise. Non-interacting system. We first present results of simulations for non-interacting case (U = 0). Note that all the standard NEGF results (including expression for noise spectrum) are exact in this case. We consider two sets of parameters representing degenerate and non-degenerate two-level systems. Parameters for the degenerate In both models, temperature was taken k B T = 1/3, Fermi energy chosen as the origin (E F = 0), and bias applied symmetrically (μ L,R = E F ± |e|V sd /2). Simulations were performed on an energy grid spanning from −120 to 120 with step 0.01. Tolerance for convergence of the Hubbard NEGF scheme was 0.001 for the difference in density matrix values (state populations and coherences) in consequent steps of iteration. The counting field for numerical evaluation of FCS was chosen as λ = 0.01. Figure 2 presents results of FCS simulations for the degenerate model, performed within the Hubbard NEGF (dashed line) and standard NEGF (solid line). The Hubbard NEGF employs self-consistent FCS as described in ref. 38. FCS results were checked against analytical expressions known for noninteracting systems. Standard NEGF (exact) results are reproduced very accurately by the approximate Hubbard NEGF simulations for both current (Fig. 2a) and zero-frequency noise (Fig. 2b). At resonant bias, |e|V sd = 10Γ 0 , conductance is maximum (inset in Fig. 2a), while noise derivative has a dip (inset in Fig. 2b), as is expected from the Landauer expression for shot noise. Results for the non-degenerate model are similar. The noise spectrum simulated for the two models within the Hubbard NEGF also reproduces exact (NEGF) results accurately (see Fig. 3).
Note that pseudoparticle NEGF (PP-NEGF), another popular Green function methodology using system many-body states, is widely employed as a standard impurity solver in nonequilibrium dynamical mean field theory simulations 57 . The inability of the methodology to yield information on full counting statistics presumably comes from its formulation within extended (unphysical) Hilbert space. At the same time, PP-NEGF can be used for simulations of current and noise spectrum. In Supplementary Information, we show that at the same (second order) level of diagrammatic perturbation expansion in system-baths couplings, PP-NEGF reproduces current, but fails to reproduce noise spectrum. The method is unable to yield correct spectral function even for non-interacting systems, which leads to failure of the noise spectrum simulation.
Weakly interacting system. We now turn to weakly interacting cases, where Γ  U/ 1 0 and perturbative expansion in this small parameter within standard NEGF should yield reasonable results 58 , while Lindblad/ Redfield QME is not expected to be accurate. Figure 4 shows the results of FCS simulations for both models performed within the Hubbard NEGF (dashed line), standard NEGF (solid line), and Lindblad/Redfield QME (dotted line) for a set of intra-system interaction strengths U. Here both Green function techniques rely on self-consistency in numerical evaluation of FCS. As expected, results are quite similar for both Green function techniques and differ significantly from the QME calculations. Naturally, Hubbard and standard NEGF results coinciding at weak U (see Fig. 4a and e) start to depart from each other at stronger interaction values. The difference is more pronounced in zero-frequency noise (main panels) than in the current (insets). There is an interesting distinction between noise results for the two models: Comparing Fig. 4c and g shows that for U = Γ 0 , the Hubbard NEGF noise result for degenerate model develops a dip at Γ e V 10 sd 0 , while no dip is observed in the non-degenerate case. The reason is difference in energetics of the two models: Single electron transitions 0a, b2, 0b, a2 for the non-degenerate model occur at −5, −4, 5, 6, while for the degenerate model the transitions are at 5, 6, 5, 6 (all energies are in units of Γ 0 ). Thus by the time bias reaches a resonant value of |e|V sd = 10Γ 0 (when physical electronic population of the system becomes significant) channel b2 in the non-degenerate model is open so that state |2〉 can be populated. Therefore no Coulomb blockade is observed. The degenerate model is Coulomb blockaded until |e|V sd = 12Γ 0 . Such suppression of shot noise due to charging effects was observed experimentally 3 . Note that the effect is seen already at U = Γ 0 /2 (compare Fig. 4b and f), where perturbation expansion in U/ Γ 0 should still be relatively accurate. However standard NEGF does not reproduce the effect because of a mean field character of treatment of U within second-order diagrammatic expansion. At the same time, the Hubbard NEGF (also considered at the second order) is accurate enough to account for the suppression. We note that diagrammatic Hubbard NEGF is a perturbative (in system-bath coupling) method, and as such is not capable to treat strong system-bath correlations. In particular, the Kondo regime is beyond the capabilities of the method. At the same time, there are many cases (including most experimental measurements of noise in molecular junctions) where strong system-bath coupling is accompanied by non-negligible intra-system interactions (e.g. electron-vibration coupling in the molecule), and yet is outside of the Kondo regime. Theoretical simulations of noise in such systems are complicated within both standard NEGF and QME approaches. The Hubbard NEGF methodology presented yields a practical tool for first principle simulations in such systems.

Discussion
Results of simulations show that the Hubbard NEGF is an inexpensive method capable to reproduce satisfactory noise characteristics of junctions over a wide range of parameters. Indeed, one has to go only to lowest (second) order in perturbative expansion in system-bath coupling to get reliable results. We think that this universality is due to basic components of the Hubbard NEGF formulation. First, as a methodology using many-body states of the system, the Hubbard NEGF accounts for strong intra-system interactions. This also makes it similar to QME considerations, which even at the lowest (second) order in system-bath coupling are successfully utilized for FCS simulations in interacting systems weakly coupled to their surroundings. In particular, such QME methods are the cornerstone of quantum thermodynamics and nonlinear optical spectroscopy considerations. At the same time, similar to the standard NEGF and contrary to both QME and PP-NEGF formulations, the Hubbard NEGF considers time correlations between single electron transitions. Note that the standard NEGF yields exact FCS results for non-interacting systems, and was successfully applied in first principles simulations in the regime of strong coupling to baths and weak intra-system interactions. These similarities with different techniques, each most suitable at opposite limits, make the Hubbard NEGF a good candidate for a relatively inexpensive approach capable of predicting noise properties in transport junctions over a wide range of parameters. It is interesting to note that in terms of structure of diagrammatic expansion, PP-NEGF, which considers correlation functions of many-body states, is closer to QME considerations. Thus, failure of the approach to predict zero-frequency noise in non-interacting systems at the same (second) order of perturbation theory where Hubbard NEGF was successful is expected. Besides, the very formulation of the PP-NEGF in extended (unphysical) Hilbert space presumably does not allow FCS formulation within the method. The latter was shown to work well in the case of Hubbard NEGF.
We show that the Hubbard NEGF yields accurate results over a broad range of parameters already at the lowest (second) order in the system-bath coupling. For non-interacting cases, Hubbard NEGF yields results similar to the standard NEGF, which is exact in this limit. For weak intra-system interactions, Hubbard NEGF is better than both the (same order of perturbation theory) standard NEGF and QME approaches. At the same time, for strong intra-system interactions, real-time perturbation theory works better than the present Hubbard NEGF consideration. Thus, careful (diagram to diagram) comparison of the two approaches is one of goals for future studies.
With respect to relevance of the methodology to experiments, we note that almost all noise measurements at molecular junctions are restricted to either non-interacting or relatively weak electron-vibration interactions 2,7,8,59 and are treated theoretically within second-order perturbation theory 31,33,34,42 . Hubbard NEGF is a scheme applicable in such experimentally relevant situations, and is capable of going beyond lowest-order perturbative schemes in intra-system interactions. It also may serve as a useful theoretical tool for investigation of (soon expected to be reported) measurements of strong light-matter interactions [15][16][17][18]     where τ 1,2 are contour variables, T c is contour ordering operator, m 1,2 are single electron transitions between many body states, and = X S S m 1 2 with |S 1 〉 containing one electron less than |S 2 〉. The Green function is obtained by solving a modified version of the Dyson type equation (see Fig. 1a  functions (dressed second-order diagrams for vertex Δ and self-energy Σ are shown in Fig. 1b and c, respectively). Self-consistency of the approach comes from the fact that both self-energy and strength operators