Effects of electronic correlations and magnetic field on a molecular ring out of equilibrium

We study effects of electron-electron interactions on the steady-state characteristics of a hexagonal molecular ring in a magnetic field, as a model for a benzene molecular junction. The system is driven out of equilibrium by applying a bias voltage across two metallic leads. We employ a model Hamiltonian approach to evaluate the effects of on-site as well as nearest-neighbour density-density type interactions in a physically relevant parameter regime. Results for the steady-state current, charge density and magnetization in three different junction setups (para, meta and ortho) are presented. Our findings indicate that interactions beyond the mean-field level renormalize voltage thresholds as well as current plateaus. Electron-electron interactions lead to substantial charge redistribution as compared to the mean-field results. We identify a strong response of the circular current on the electronic structure of the metallic leads. Our results are obtained by steady-state Cluster Perturbation Theory, a systematically improvable approximation to study interacting molecular junctions out of equilibrium, even in magnetic fields. Within this framework general expressions for the current, charge density and magnetization in the steady-state are derived. The method is flexible and fast and can straight-forwardly be applied to effective models as obtained from ab-initio calculations.


I. INTRODUCTION
Miniaturization as a performance-enhancing concept in microelectronics may be advanced by the introduction of nanoscale molecular devices in what is today known as the concept of molecular electronics [1]. Recent years fostered fascinating advances in experimental control of fabrication [2], assembling [3], as well as contacting [4,5] on a molecular level. These achievements in combination with ever-improving measurement techniques [6] lead to a plethora of important insights into the basic mechanisms of electrical transport across molecular junctions [7]. Understanding these transport characteristics is a major focus of today's experimental as well as theoretical ventures and establishes the very basis for possible future device engineering.
A relevant and still relatively simple molecular junction is comprised of a ring-shaped molecule in between two metallic leads. For the particular molecule we have in mind the aromatic complex benzene (C 6 H 6 ). Such setups have been realized using the mechanically controllable break junction technique (MCBJ) [5] and are stable on time scales required for transport experiments. Measurement of the transport characteristics has been achieved for benzene bound by thiol anchor groups to gold electrodes [8][9][10], as well as for benzene directly connected to platinum leads [11]. These experiments typically grant access to the current-voltage characteristics, conductance, higher derivatives of the current, or shot noise, albeit often in a statistical way [9]. It is likely that molecular ring junctions will find technological applications in the foreseeable future in the form of singleelectron transistors [12,13], in quantum interference (QI) based * martin. nuss@tugraz.at Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. electronics [14], or as data storage devices [15]. Naturally such junctions are dominated by quantum mechanical effects in an out-of-equilibrium context, which promotes predicting the outcome of experiments to a highly nontrivial theoretical task.
For benzene-based molecular junctions, general understanding of the noninteracting device is available in literature [16,17]. Electronic transport in π -conjugated systems is special due to QI effects [18,19]. Magnetic fields add to the rich Aharonov-Bohm physics [20,21] of quantum ring structures by inducing, for example, persistent currents [22]. Electronic correlations are important [23] due to the confined geometry and have been recently studied in equilibrium and linear response using the Hirsch-Fye quantum monte carlo (QMC) method [24], as well as dynamical mean-field theory (DMFT) [24,25]. While the basic features of QI effects can be understood from noninteracting calculations, the interplay of QI and electronic correlations can become nontrivial [26][27][28][29][30]. The remarkable property of negative differential conductance has been reported and explained in devices considering electron-electron interactions [27,31] using generalized master equation approaches. In addition, Green's function techniques have been applied within various approximations, especially within a combination with ab initio techniques [32]. In gated devices even nonperturbative many-body features like the Kondo effect [33] have been reported [3,4,30,34,35].
Investigating transport characteristics for molecular junctions out of equilibrium [5,36] is a very active field in modern theoretical physics. Although much progress has been made lately [37][38][39][40], reproducing experiments [10,41] in a qualitative and quantitative way remains an elusive goal. In going beyond semiclassical treatment [42] and combining nonequilibrium Green's function (NEGF) methods [43] with density functional theory (DFT) [44,45], much effort is devoted to the study of molecular junctions in a combined DFT+NEGF scheme [46][47][48][49][50][51][52][53][54][55][56]. Such ab initio approaches can be considered today's gold standard for weakly correlated molecular junction calculations. However, they tend to overestimate currents in comparison to highly accurate data [57] of experiments, especially in interesting cases of low device-lead coupling [58]. This overestimation can often be attributed to too simplistic treatment of electron-electron interactions [59,60]. In the worst case, when many-body interaction effects beyond mean field become important, these methods might even predict nonphysical results such as low conductivity when the device would really show transparent conduction due to many-body effects [61].
Therefore it is desirable to extend these techniques with methods to tackle the effects of strong correlations in an out-of-equilibrium context [57,[62][63][64][65][66][67]. Many-body effects of electron-electron interaction can be studied by making use of simplified model Hamiltonians. By using a modelbased approach, however, results may strongly depend on the model parameters, which are usually notoriously difficult to estimate. For benzene, physically relevant interaction parameters have been obtained in literature by fitting experimental spectra [68,69] which provide the basis for many theoretical works [26,27]. Such parametrizations are typically not unique, since considering different interaction matrix elements may yield again a good agreement with a specific set of experimental data but might promote very different physical mechanisms. Another approach, yielding model parameters, is to determine them in an ab initio way from first-principles calculations [70,71]. Again, these parameters are subject to changes due to screening effects, when the molecular device is finally coupled to leads or embedded in an environment [72,73].
In this work we study the many-body effects of electronelectron interaction beyond mean field in the steady state (sts) of a ring-shaped molecular junction under voltage bias in a magnetic field. A simple Pariser-Parr-Pople model [74,75] description for the atomic carbon p z orbitals is employed for a charge-neutral device coupled to two metallic leads. We address the question of how electron-electron interactions influence device properties in a magnetic field and identify their fingerprints. We do so by assessing naturally important quantities for molecular junctions, which are sts observables under applied bias voltage. We focus on studying the dependence of the sts currents, sts charge density, as well as sts magnetization on electron-electron interactions and a magnetic field, as well as their dependence on devicelead coupling. The purpose of this work is to concentrate on the qualitative effects of electronic correlations, rather than seeking quantitative agreement with experiments. A detailed, quantitative comparison to experiments would of course require the inclusion of additional electronic bands, mechanical vibrations, as well as charging and temperature effects. On the other hand, correlation effects would be very difficult to treat at a comparable level of accuracy in such a rich model. In this work we focus on the interaction effect on the most simple model for a benzene junction to identify basic mechanisms. However, the presented approach can easily be generalized to include anchor groups or specific lead geometry.
To solve the interacting system out of equilibrium, we use cluster perturbation theory [76][77][78], which amounts to treating the intercluster hopping within first-order strong coupling perturbation theory. This approach is, in principle, refinable by considering larger cluster sizes. Here, we target directly the sts, which is formulated in the framework of sts cluster perturbation theory (stsCPT). We work out general expressions for the sts density matrix and all bond currents in presence of a magnetic field in terms of stsCPT single-particle Green's functions. This framework provides a flexible, versatile, and easy to use method for evaluating sts observables for interacting molecular junctions out of equilibrium in magnetic fields. It can easily be generalized to other junctions and can straightforwardly be combined with ab initio calculations.
Our results indicate that while properties of the electronic structure of the bare device, including symmetry considerations and degeneracies, have a dominant influence [27] on the sts behavior, electron-electron interactions effectively renormalize conduction characteristics like threshold voltages and plateau currents in magnetic fields. We find that the sts charge density becomes strongly renormalized beyond mean field and identify signals in the sts magnetization. This paves the way for an even more sound and trustworthy interpretation of results from sophisticated ab initio-based methods, as well as model-based calculations at fixed model parameters. Furthermore, we reexamine effects introduced by magnetic fields in context with electron-electron interactions, lead-induced broadening, as well as lead electronic structure.
This paper is organized as follows. We present the theoretical model in Sec. II and introduce the formalism for our calculations in Sec. III. We outline how to systematically improve results obtained with the presented method and compare obtained data with existing calculations in Sec. III D. Results for the sts currents, sts charge density, and sts magnetization are provided in Sec. IV.

II. MODEL
We consider the effect of electron-electron interactions on the electric transport through a benzenelike aromatic molecular ring. A simple starting point is provided by a tight-binding approach for the molecule coupled to the left and right metallic leads [see Fig. 1 (top)]: The benzene ring is modelled by considering one atomic p z orbital φ iσ per carbon atom, yielding a six-atomic-orbital Pariser-Parr-Pople [74,75] [(extended) Hubbard [79]] Hamiltonian, where i,j ∈ [1,6] enumerates the six ring orbitals in a clockwise fashion. The six nearest-neighbor bonds are denoted ij . Elementary fermionic operators d iσ /d FIG. 1. (Color online) (Top) Illustration of the device setup: A planar molecular ring is connected to two metallic reservoirs. The setup is placed in a perpendicular magnetic field B. A bias voltage V B is applied between the left and right lead. (Bottom) Bare singleparticle energy levels of the noninteracting, disconnected molecule. Bias voltages required for the specific levels to contribute to transport are indicated on the right. The HOMO as well as the LUMO are doubly degenerate (d = 2) for B = 0. the particle number operator is defined in the usual way aŝ The on-site energy of the ring orbitals is comprised of the bare on-site energy d = −1.5 eV [80] (with respect to the leads), a Zeeman term σ B 2 , and a correction canceling the mean-field contribution of the on-site interaction U and the nearest-neighbor density-density interaction W : − U 2 − 2W . Literature [68,69] provides an optimal parametrization for Eq. (2) when not connected to leads. Fitting to excitation spectra [81], the authors of Ref. [68] find that benzene is best described by a nearest-neighbor overlap integral of t d ≈ −2.5 eV and an on-site interaction of U ≈ 10 eV [82]. Such a model is frequently used in literature [26,27], along with approaches where the interaction parameters are determined in an ab initio way from first-principles calculations [70,71]. We expect the values of U and W to be substantially reduced when the molecule is connected to the leads due to screening from the band electrons [72,73]. Therefore we set in the following t d = −2.5 eV and discuss values of U 9 eV as well as W 3 eV.
We consider the effects of a magnetic field B which is applied perpendicular to the plane spanned by the molecular ring. This is described by a Peierls phase [83] (B) as well as by a Zeeman term.
For simplicity, we model the right (α = R) and left (α = L) leads by semi-infinite tight-binding chains, where the fermionic operators in the lead orbitals c iασ /c † iασ are defined in a standard way. The resulting semicircular electronic density of states (DOS) is centered around the leads' on-site energy α , which is fixed to zero for zero bias. Application of a bias voltage V B is done by shifting the lead's chemical potential (strictly speaking, at infinity) and on-site energies to α = μ α = ∓ V B 2 , so that the leads are kept half-filled. We use a large lead nearest-neighbor hopping t α = −6.0 eV (independent of α), which implies a quite large bandwidth D = 24 eV, so that most of our results are comparable to the wide-band limit [80,84]. We will also discuss and compare to differently shaped lead DOS when effects on the sts properties are to be expected.
The coupling between the ring and the leads is described by a single-particle hopping t to the left and right lead. The leads are connected to the molecule in the so-called para(1,4), meta(1,5), or ortho(1,6) configuration. Here the two numbers in the braces label the ring position at which the (left, right) lead is attached, denoted (1, x ∈ {4,5,6}) in the following: This molecular-lead setup provides a good description of, e.g., two widely used experimental techniques: scanning tunneling microscopy (STM) [6] or MCBJ [8]. Then the STM tip/substrate material or break junction material enters inĤ lead through the DOS, whileĤ ring-leads describes the tunneling from the molecule into the experimental device. In most of the present work, we focus on a typically small value for the molecule-lead coupling of t = −0.05 eV [16], which enhances both the effects of electronic correlations and of the magnetic field. This leads to a level broadening of the order of ∝ πt 2 DOS leads (ω = 0) ≈ t 2 |t| ≈ 10 −4 eV. To study the effects of this lead-induced broadening, we also present results for larger t .
For the sake of simplicity we neglect the effects of mechanical vibrations. The coupling of the electronic degrees of freedom to vibrations can lead to a renormalization of conductance thresholds as well as current plateaus as discussed in Ref. [85] and will be addressed in more detail in future work [86]. Furthermore, we do not discuss charging effects due to the connection of the leads.

III. METHOD
We aim at calculating sts properties of an interacting molecular device. Here we focus on the current between each bond as well as the charge distribution and magnetization.
We consider two currents of special interest: (i) the total transmission current j t and (ii) the circular current j c . In the presence of a bias-induced transmission current, the two ring directions carry different currents. This leads, in principle, to an ambiguity in the definitions of the circular current. However, as discussed by Rai et al. [80], the most natural expression is the one which is directly related to the current-induced magnetic flux through the ring. According to Biot-Savart's formula [87], this is given by the average current obtained by weighting the current flowing through each segment by its length. In our case, in which there are two contacts dividing two ring segments i = 1,2 of lengths L i (in units of the lattice constant), in which currentsj i flow (say, clockwise), we have In this work, we consider couplings to the leads in para: L 1 = 3,L 2 = 3, meta: L 1 = 4,L 2 = 2, and ortho: The total transmission current is given by which equals the inflow as well as the outflow at the bonds connecting the leads to the ring by virtue of the continuity equation.
The sts charge distribution n and magnetization m can be obtained from the sts single-particle density matrix D σ ij , which also encodes all information about the sts current.

A. Nonequilibrium Green's functions
One way of obtaining the sts density matrix is by making use of nonequilibrium Green's functions in the Keldysh-Schwinger [88][89][90] formalism, where G R/A/K are matrices in orbital/spin space and functions of two time coordinates τ [91]. R denotes the retarded, A the advanced, and K the Keldysh component of the single-particle Green's function in Keldysh space G. Since we focus on the sts, time translation invariance applies and we can express the Green's functions in frequency space ω.
The sts density matrix D σ ij is obtained from the sts Green's function G(ω) in matrix notation: The sts current through orbital i in the presence of a magnetic field B can be obtained by the time derivative of the total particle number of orbitals on either side of i [92]. In terms of equal time [93] correlation functions, one can express the current by Keldysh Green's functions ,i (τ,τ )). Generalizing the notation to arbitrary indices and keeping in mind that a definition of current only makes sense for nearestneighbor orbitals, an expression in terms of the sts density matrix D becomes available: where h σ ij denotes the single-particle part ofĤ [Eq.
(1)] in the orbital/spin basis. The expression is purely real because D and h are Hermitian.
Since we are concerned with electron-electron interactions, evaluating the needed Green's functions G is not possible in general.
Within this approach, the thermodynamically large system Eq. (1) is split into individually, exactly solvable parts at time τ → − ∞. The single-particle Green's function g(ω) for each of these parts (clusters) is obtained by analytic or numeric means. The coupling between these parts is switched on at a later time τ 0 using the intercluster (perturbation) matrix T , which holds the couplings between the disconnected parts. The sts Green's function of the full system in the stsCPT [94,95] approximation is given by where 1 1 is the identity in Keldysh space. We use lower case g for the single-particle Green's function of the initially decoupled equilibrium system, while upper case G denote the sts Green's functions of the fully coupled system. The stsCPT approximation made here is to replace the self-energy G of the full system by the self-energy g of the cluster. This amounts to a first-order strong coupling perturbation theory in intercluster terms T . The appealing aspect of this approximation is that it becomes exact in each one of three different limits: (i) for T = 0, (ii) for U = 0, or, in principle, (iii) for an infinite cluster.
In our case, the system at τ → − ∞ thus consists of the molecular ringĤ ring [Eq. (2)], disconnected from the two leadsĤ α lead [Eq. (3)], i.e., t = 0. One can compute the retarded single-particle Green's function g R of the interacting ring by a standard Lanczos approach, [96,97] and that of the noninteracting leads analytically [98,99]. The advanced component is available by the identity g A = (g R ) † . The Keldysh component g K of the initial systems, which are separately in equilibrium, can by obtained by the relation [100] where p FD (ω,μ,β) = 1 e β(ω−μ) +1 is the Fermi-Dirac distribution function, β denotes the inverse temperature, and μ the chemical potential.
The state of the molecular device at τ → − ∞ is half-filled and unpolarized for all parameters under discussion in this work. While, in principle, the sts results should not depend on the initial state of the finite-size central part, the approximate nature of the calculation leads to some dependencies on the initial temperature and chemical potential of the central region in the presence of interactions. In general, we take the zero temperature ground state as a starting (τ → − ∞) point, which leads to very good results in most cases. There are some exceptions, in particular when degenerate states are involved. Therefore blocking effects as discussed in [26,27] cannot be observed within our approach. However, these are expected to occur for B = 0 up to very small fields only [27], since B breaks the degeneracy responsible for these effects.
It is sufficient to calculate two 6 × 6 matrix Green's functions for the central molecule (six orbitals times spin) and four scalar Green's functions representing the contacting orbital of the two leads (times spin). Making use of Eq. (7), these three initial systems are perturbatively connected by the coupling Hamiltonian [Eq. (4)] in the single-particle basis to obtain the sts single-particle Green's function in Keldysh space G(ω).

C. Evaluation of steady-state observables
Using the G(ω) obtained via stsCPT, we rewrite the expression for the sts density matrix Eq. (6) in the corresponding language, where [A,B] − denotes the standard commutator, and Einstein's summation convention is implicit. The occupation matrix P σ (ω) corresponds to the second part of Eq. (8) and is diagonal: P σ ij (ω) = 2δ ij p FD (ω,μ i ,β i ), where μ i , β i are the chemical potential and inverse temperature of site i at τ → − ∞. This expression is used to evaluate all sts observables as defined in Sec. III A.
The integrals in Eq. (9) are nonzero only between μ min and μ max , which is due to the P dependence and renders the numerical evaluation of Eq. (9) much more favorable than directly evaluating Eq. (6). For a numerical evaluation, a high-precision adaptive integration scheme is necessary [101], especially when it comes to resolving small differences in bond currents as induced by magnetic fields [see Eq. (5)].

D. Quality of the stsCPT approximation and systematic improvements
In the following we assess the quality of the stsCPT approximation and discuss systematic improvements for the treatment of correlations. Upon testing we found these improvements to be unimportant for this particular setup. They can, however, be be of considerable importance in systems exhibiting a more complicated quantum mechanical ground state [78].
As discussed above, one can, in principle, systematically improve results by enlarging the central region to include a certain number of orbitals of the leads in addition to the molecular ring. This means that at τ → − ∞ the system is disconnected at some bond t of the lead's chain. This sets much higher computational demands because the Hilbert space of the interacting part of the system grows exponentially in system size. In order to assess the accuracy of our results, we have investigated the convergence by adding one to three sites of each lead to the molecule and found negligible deviations in the sts results. This also allows us to conclude that data presented in this work are accurate and no serious error is being made due to the above-mentioned approximations. Other expansions in the framework of Keldysh perturbation theory are discussed in Ref. [28] in this context.
Another way of improving the approximation for the self-energy is to include a self-consistent feedback within the nonequilibrium variational cluster approach [77,78]. We find this to be a major improvement for "gapless"systems. In this work we consider a system which has a large highest-occupied to lowest-unoccupied molecular orbital gap (HOMO-LUMO, respectively), and we found the improvements introduced by the self-consistency to be less important. Therefore, in this work we are restricted to the nonvariational stsCPT. Moreover, we limit ourselves to zero temperatures, although finite temperatures are easily accessed.

IV. RESULTS AND DISCUSSION
We present results for the sts properties of para-, meta-, and ortho-connected benzene. We start by discussing lead-induced broadening effects on the interacting ring within stsCPT for a broad range of molecule lead coupling strengths. Then we focus on the interesting case of small molecular lead coupling, where the effects of a magnetic field B and electronelectron interactions become important. The current-voltage characteristics of the total transmission current j t and the circular current j c are presented, including a discussion of effects induced by the shape of the leads' DOS. We identify the behavior of threshold voltages V T and maximum reachable currents. A detailed view on the highly interesting bias region of high variation in the current (around V T ) will be provided, identifying the effects of electron-electron interactions on the current signal shapes. Finally, we examine the sts particle density n and magnetization m in the device.
We note that the presented method becomes exact in the noninteracting limit (U,W ) → 0 and therefore has to agree with results for the noninteracting device presented by Rai et al. [16,17], obtained by calculations based on the Landauer formula [102] (not shown). We go beyond their detailed description in terms of transmission coefficients and QI effects by discussing effects of electron-electron interaction and by adding a Zeeman term to the discussion of the magnetic field. As discussed above, our approach may become unreliable in the presence of degeneracies in transport states, which occur in this model for B = 0 (see Sec. III B). Therefore we do not observe the current blocking state which has been reported [26,27] in asymmetrically connected junctions at finite bias voltage for U > 0 and B = 0.

A. Lead-induced broadening
Lead-induced broadening is obviously important for nottoo-small values of the molecule-lead coupling t , as we show in Fig. 2 (top). On the other hand, there are important broadening effects even for small t , such as the ledge in the 155139-5 transmission current visible in Fig. 2 (bottom) for the meta and ortho setup, which is absent in the para case [103]. It is clear that an accurate treatment of broadening effects is highly desirable, even more so when electron-electron interactions come into play. We find that the total transmission current scales with |t | 2 , while the circular current does not show a clear scaling as a function of t . A technique commonly employed in the field of molecular electronics is based on generalized master equations [26,27]. Lead-induced broadening is in general difficult to take into account accurately within such methods. stsCPT is able to capture broadening effects even in the interacting molecule, thus rendering stsCPT an interesting, complementary approach to the generalized master equation technique.
Note that we did not discuss the effects of a magnetic field here since they are fully washed away by the large t (for reasonable fields on the order of a few Tesla). Effects of electronic correlations become increasingly important for smaller t . For these two reasons, in the rest of this work we focus on t = −0.05 eV because we aim at studying the interplay of a correlated device in a magnetic field.

B. Field-dependent current-voltage characteristics
To ease navigation in the rest of the text, we provide a bird's eye view on the current-voltage characteristics in a large voltage regime V B = [0,9] V in Fig. 3. Here, effects of the finite bandwidth on the transmission current are still small. Since the magnetic field induces a very small energy scale in comparison with the on-site, hopping, and interaction energies in the Hamiltonian [Eq. (1)], we will be zooming into regions V B = ±0.001 V around selected voltage points. We note that the absolute position on the bias axis strongly depends on the bare molecule-lead potential difference ( d = −1.5 eV) which governs the first threshold voltage.
We start out by comparing the total transmission current j t and the circular current j c of para-, meta-, and ortho-connected benzene for the noninteracting case, as well as for U = 9 eV, with and without magnetic field B. j t shows multiple threshold voltages V T in a plateaulike structure for all three device setups. At exactly these thresholds, a signal spike in j c is observed. Those signals in j c are of considerable magnitude, which are tens of microamperes, as compared to j t , which is on the order of tens of nanoamperes, the latter being FIG. 3. (Color online) Current-voltage characteristics of para-(top), meta-(center), and ortho-(bottom) connected benzene. We plot the total transmission current j t /nA in the left and the ring current j c /μA in the right column. For every device setup, data for the noninteracting case (blue) is shown in comparison to results for an on-site interaction strength of U = 9 eV. We compare data for the zero-field case (solid lines) to curves obtained for a magnetic field of B = 2 T (dash-dotted lines). mainly determined by the molecule-lead coupling t . Such thresholds in j t and large signals in j c seem to be a generic feature of molecular ring devices [104][105][106] and can be associated with nearly degenerate pairs of states with opposite angular momentum whose contribution to the net current is rendered small by destructive interference [80]. We discuss the interaction dependence of the threshold voltages V T in detail in Sec. IV D. Here we note that one effect of electron-electron interaction is to shift the conduction thresholds V T (see A in Fig. 3). Notice that this is a pure correlation effect, since the mean-field contribution has already been subtracted.
The magnitude of the plateaus in j t is fixed by the amount of current which can be carried by the molecular level starting to participate in transport at the corresponding V T . Due to symmetry, j t in the para device is quite a bit higher than in the corresponding meta or ortho devices, the latter two being comparable. In particular, we find that the plateau current magnitudes are the same in the noninteracting meta and ortho devices but not when interactions are present (see also Sec. IV E). The transmission current of the para setup is independent of the magnetic field, while in the meta and ortho device, j t grows until it seems to saturate, at least for 155139-7 the magnetic field magnitudes discussed in this work (see Sec. IV E). Increasing on-site interaction U or density-density nearest-neighbor interaction W (not shown) always decreases the total transmission current (see B in Fig. 3, discussion in Sec. IV E).
The spikes in j c at V T (which are split by the Zeeman field, see Sec. IV F) are of considerable magnitude, which is tens of microamperes. Similarly to j t , increasing on-site interaction U or nearest-neighbor density-density interaction W reduces the magnitude of these signals. In the para device the ring currents exactly cancel for B = 0 T due to device symmetry, while in the other two setups circular currents are also present in the field-free system. Upon increasing the magnitude of the magnetic field, all currents show saturation (see Sec. IV E). Note that for half-filled systems in the large U limit [107] the current is unaffected by magnetic flux because electron motion becomes severely hindered [22].
We find that effects of electron-electron interactions (shifting of conduction thresholds and decreasing plateau currents in j t and signal currents in j c ) are smooth and do not depend on the specific type of interaction (on-site or nearest-neighbor density-density).

C. Effects originating from the electronic structure of the leads
Instead of the often imposed wide-band DOS of the leads, we use a semicircular lead DOS with a large bandwidth (t = −6 eV). As expected, our setup mimics a wide-band limit regarding j t for low bias voltages where only a small bending down of the current can be observed due to the negative curvature at ω = 0 eV of the leads' DOS [see Fig. 3 (left)].
However, it turns out that j c is highly sensitive to even slight variations in the leads' DOS. Let us first discuss the behavior of j c in detail. In the para setup, the magnetic field induces peaks in j c at V T , while in between it falls back to plateaus in the nanoampere range. In the meta and ortho setup, also at zero field, a finite ring current exists which can be amplified or suppressed by the magnetic flux, depending on its direction according to Lenz's law [87]. These two setups show a behavior of j c which is different from the para case. Without magnetic field, j c shows plateaus (∝ V 2 B ) in the microampere range between the peaks. When a magnetic field is imposed, the two conducting states of the molecule at k = ± 2π 3 and energy ω ≈ 1 eV split up [see Fig. 1 (bottom)]. In our setup j c shows a linear growth between the signals starting at the first threshold voltage. This large effect can be traced back to different occupation of the k = ± 2π 3 states. Calculations for a constant lead DOS do not show this linear increase in circular current. Instead, constant plateaus (and therefore a moderate imbalance in the population of the k = ± 2π 3 states) are obtained. These plateaus have a magnitude which j c acquires for the semicircular DOS right after V T . The crucial parameter here is the curvature of the lead DOS, which in the end amplifies the population imbalance and renders it bias dependent. For a semicircular DOS, the coupling of the two levels at ω ≈ 1 eV to the high-bias and low-bias lead is of different magnitude (effectively becomes different for the two leads), and this difference grows linearly with increasing bias voltage. To our knowledge, this effect has not been discussed before and is probably difficult to observe experimentally, especially since it strongly depends on the features of the DOS of the leads. Moreover, electron-electron interaction does not play a role here. However, if a suitable system can be constructed, bias voltage could be used to linearly tune the circular current over a significant current magnitude.

D. Current signal position
The transmission current-voltage characteristics consist of plateaus with steep jumps between them [see Fig. 4 (left)]. These signal positions at threshold voltages V T are independent of the setup used (para, meta, or ortho connection) and furthermore, are independent of the magnetic field B [108]. We plot data for meta-connected benzene for B = 1 T as a representative in Fig. 4 (left). The width of the current signals is ≈ 10 −3 eV (see Sec. IV F) due to the small lead-ring coupling t .
In the interacting device one can still interpret our data in terms of the single-particle excitations of the molecular system obtained from exact diagonalization of the molecular Hamiltonian [Eq. (2)]. Comparing the HOMO-LUMO gap (U,W ) to the data for the sts current (see Fig. 4), one finds a linear relation between the gap and the threshold voltage V I T = − 3 eV. Note that the constant depends on the position of the molecular on-site energy with respect to the leads, but the effects of interactions in are universal for weakly coupled charge-neutral devices. From our data we find a change in the HOMO-LUMO gap due to interactions beyond mean field (U )[eV] = 5[eV] + (0.02 ± 0.01) e (0.43±0.03)U [eV] . The jump to plateau III exhibits the same dependence on U as the jump to plateau I but located at a different position: V III T = + 3 eV. Interestingly enough, these two thresholds are monotonic in interaction strength but the threshold II in between is not. It shows a strongly nonmonotonic behavior [see Fig. 4 (right)].
The functional dependence of the thresholds on on-site U is the same as for nearest-neighbor density-density interactions W (not shown). It is clear that the HOMO-LUMO gap grows much quicker with increasing W than with increasing U .
We conclude that the voltage thresholds in the parameter regime of the device under study, which consists of (i) small molecule-lead coupling t t d t and (ii) a charge-neutral setting, are perturbatively accessible and can even be inferred to a high degree of accuracy from the single-particle spectrum of the isolated device. However, our results show that they can be subjected to nonmonotonic behavior as a function of interaction strength.

E. Current signal magnitude
We now turn to the analysis of the maximum current j max t in the vicinity (±0.1 V) of the first threshold voltage V I T . As noted before, the para setup does not show a magnetic field dependence in this channel for fixed interaction parameters due to device symmetry [thick black line, stars in Fig. 5 (left)]. In the meta and ortho devices, we find that the maximum current increases as a function of magnetic field |B| until saturation for the considered range of magnetic field strength for all bias voltages [see Fig. 5 (right)].
Note that the eigenstates of the disconnected, noninteracting molecule in a magnetic field are n The inherent circular current driven by the magnetic field is given by j inh c = − e 2π m e 6 2 n occ,σ n σ . From this expression it is clear that in this case the circular current magnitude is bounded from above because the magnetic field just redistributes occupied momenta of compact states. This analysis holds for all finitesize quantum ring devices, but not for the corresponding one-dimensional field theory [110].
Depending on interaction strength U , we find that the maximum transmission current is monotonically decreasing (roughly ∝ const. − U 2 ) in all device setups [see Fig. 5 (left)]. The same effect is observed for nearest-neighbor density-density interactions W (not shown). The maximum transmission currents of the meta and ortho device are identical in the noninteracting device but start to differ when interactions are turned on.

F. Current signal
A zoom-in to the signal in the sts currents at the position of their respective first threshold voltages (I) is provided in Fig. 6. We compare the currents for the noninteracting system with those of on-site interaction strength U = 9 eV. We observe similar signal shapes for all values of interaction strength (U,W ). These signals are, however, shifted in bias depending on the electron-electron interaction parameters (U,W ) and the signal magnitude is decreased. In Fig. 6 we compare the signals of the noninteracting and interacting setups whereby the horizontal axis has been shifted in order for the threshold voltages to coincide.
The total transmission current j t does not depend on the direction of the magnetic field in the para setup. The magnitude of the plateaus is, however, increased with increasing |B| in the meta and ortho setup, while staying on the field-free value in the para setup. The effects of the Zeeman term are visible in j t due to a splitting ∝ |B| of the transition to the next plateau into two subplateaus around the transition points V T . The circular currents j c in the para setup depends on B in a symmetric way, i.e., reversing the direction of B reverses just the sign of j c . The meta and ortho setups show a different magnitude in the circular current for the same |B| due to inherent circular currents which are either amplified by the magnetic field or suppressed, depending on its alignment. The Zeeman splitting introduces a double-peak structure around V T which is proportional to the magnetic field magnitude |B|. This splitting is introduced in a symmetric fashion around V T . The circular currents j c grow with increasing magnetic field until a saturation is reached. Besides the decrease in circular current magnitude, we observe a more pronounced Zeeman splitting with increasing interaction strength.
We again find a generic dependence on interactions beyond mean field, i.e., using on-site or nearest-neighbor densitydensity interaction, which might even permit extracting interaction parameters for model Hamiltonians by carrying out transport measurements and fitting the position and height of sts currents at threshold voltages.

G. Steady-state local charge density
The sts charge density is shown in Fig. 7 for all three device setups for an interaction strength of U = 9 eV and a magnetic field of B = 2 T as an example for its generic behavior. It exhibits features at the same threshold voltages as the sts current, where it either decreases or increases in a plateaulike fashion. In the para setup [ Fig. 7 (left)] we observe an increased charge density at the orbitals at which the leads are connected and a strongly reduced charge density in the other orbitals of the molecule due to interactions, i.e., for the noninteracting molecule all densities would start at one. The sts charge density of the meta and ortho setups [ Fig. 7 (right)] are symmetry related. They show the same sts charge density, except for the fact that two orbitals exchange their respective charge density each.
The suppression of occupation is most pronounced in the orbitals in the "interior" of the molecule, respecting symmetry. These are orbitals 2,3,5,6 in the para device, orbitals 3,6 in the meta device, and orbitals 2,5 in the ortho device.
For smaller interaction strength, the overall shape of the curves does not change. However, we find that in all setups, for on-site interaction strengths U 5 eV, all local densities are 1 in the low-bias regime up to the first threshold voltage, n = 0, and depend on the interaction strength beyond that point. For interaction strength U 6 eV, n depends on interaction also in the low-bias region.
We note that for increasing molecule-lead coupling t , the magnitude of the sts local charge density stays the same but the transitions become washed out (not shown). Furthermore, with increasing bias voltage as well as in the vicinity of the threshold voltages, a splitting of the degenerate orbitals occurs.
The main effect of interactions beyond mean field is to suppress the sts charge density respecting the symmetry of the isolated molecule. Such substantial charge redistribution in sts transport due to interactions is an important effect to be considered when discussing results from self-consistent DFT+NEGF calculations.

H. Steady-state magnetization
In our calculations a Zeeman term is included. Therefore the magnetic field induces a sts magnetization in the vicinity of the threshold voltages V T . We consider just the paramagnetic (spin) contribution, since it is the dominant term. In Fig. 8, we show data for all three device setups in the vicinity of the first voltage threshold (I) for the interacting case. Increasing the magnetic field turns the resonancelike structure in j c around V T into a structure which consists of two maxima on the positive FIG. 7. (Color online) Ss local charge density of the para-connected device (left) as well as the meta-, and ortho-connected devices (right). We show data for a device with on-site interaction of U = 9 eV in a magnetic field of B = 2 T. The drawings of the device at the top serve as a legend for the line colors to identify the respective orbitals. n denotes the shift of the local density away from n i = 1 in the low-bias region. and one on the negative axis in between. The effect of changing the direction of the magnetic field while keeping |B| constant is to mirror the sts magnetization curves around the x axis (not shown). Without interactions the effect is quite similar and only of slightly different magnitude.

〈 〉
We find that the sts magnetization is a highly sensitive quantity as a function of bias voltage (in the submillivolt regime) in the region around V T where it is possible to alter the direction as well as the magnitude of the magnetization by ≈ ±10 −2 /μ B .

V. CONCLUSIONS
We study effects of electron-electron interactions in a benzenelike ring-shaped molecule, subject to a finite bias voltage induced by two metallic leads, as a function of an applied magnetic field. We make use of a Hubbard-typemodel-based description of such a device in the chargeneutral regime. Using steady-state cluster perturbation theory, observables have been computed for para, meta, as well as ortho setups.
Results for the total transmission current and circular current as well as the steady-state charge distribution and magnetization have been presented. By studying physically relevant regimes of electron-electron interactions in addition to an applied magnetic field, we describe the effects of electron-electron interactions on the steady state beyond the mean-field level. We found that these are to shift voltage thresholds and to decrease the magnitude of currents. Additionally, interactions lead to deviating currents in the meta and ortho setup which were comparable in the noninteracting system. The steady-state charge distribution becomes strongly renormalized by interactions respecting the symmetry of the isolated molecule. Due to the Zeeman effect, we obtain a steady-state magnetization which is highly sensitive to bias voltage.
Our results may help to validate model calculations at fixed interaction parameters and contribute to the understanding of sophisticated ab initio-based transport calculations. They might even contribute to designing empirical formulas for junction engineering. Our results indicate that the main effect of interactions is to renormalize voltage thresholds and current magnitudes. Care has to be taken in discussing symmetry relations of meta-and ortho-connected devices. Furthermore, we showed that the charge density is sensitive to electronelectron interactions and becomes strongly renormalized with every additional electronic level contributing to transport. This fact has to be accounted for in self-consistent approaches. We find that insights into the relevant interaction parameters, as well as their magnitude, could be extracted from transport measurements, provided the on-site potential of the molecule with respect to the leads is known.
We presented general expressions for steady-state observables in the language of steady-state cluster perturbation theory. The presented formalism is flexible and simple to apply to a broad range of molecular junctions. Its approximate nature is systematically improvable, and the method does work in all parameter regimes. The method is especially interesting in parameter regimes which exhibit extremely high time scales, rendering approaches based on time evolution not feasible, and when lead-induced broadening effects become important, which renders generalized master equation techniques difficult. The presented approach is, however, unable to capture interaction-mediated interference effects, such as current blocking in asymmetrically connected junctions, when degeneracies in transport channels are present. This can be overcome by a more elaborate construction of the "starting" (τ → − ∞) cluster state. Work along these lines is in progress.