Efficient simulations with electronic open boundaries

We present a reformulation of the Hairy Probe method for introducing electronic open boundaries that is appropriate for steady state calculations involving non-orthogonal atomic basis sets. As a check on the correctness of the method we investigate a perfect atomic wire of Cu atoms, and a perfect non-orthogonal chain of H atoms. For both atom chains we find that the conductance has a value of exactly one quantum unit, and that this is rather insensitive to the strength of coupling of the probes to the system, provided values of the coupling are of the same order as the mean inter-level spacing of the system without probes. For the Cu atom chain we find in addition that away from the regions with probes attached, the potential in the wire is uniform, while within them it follows a predicted exponential variation with position. We then apply the method to an initial investigation of the suitability of graphene as a contact material for molecular electronics. We perform calculations on a carbon nanoribbon to determine the correct coupling strength of the probes to the graphene, and obtain a conductance of about two quantum units corresponding to two bands crossing the Fermi surface. We then compute the current through a benzene molecule attached to two graphene contacts and find only a very weak current because of the disruption of the $\pi$-conjugation by the covalent bond between the benzene and the graphene. In all cases we find that very strong or weak probe couplings suppress the current.


I. INTRODUCTION
Atomic scale computer simulations of nanoscale systems of necessity have to approximate the environment that the system finds itself in as it is of unlimited size. One way to incorporate a model environment is through the boundary conditions of the system being treated explicitly. Here we focus on the boundary conditions for the electrons. Traditional choices include: free boundaries, where the system is treated as an isolated cluster in vacuum; periodic boundaries, where the system plus its near environment are repeated periodically to make an effectively infinite system; and open boundaries, where the system is finite but electrons can enter and leave as though connected to an external reservoir. The correct choice of boundary conditions is determined by the problem being addressed.
There exist very efficient algorithms for free and periodic boundary atomistic simulations [1], and these will not be considered further here. Open boundaries are important for a number of problems [2], and mature open boundary codes also exist [3][4][5][6]. However, relative to free and periodic boundaries, they tend to be more computationally expensive to implement, and simulations can require more human effort to set up. These technical con- * a.horsfield@imperial.ac.uk siderations tend to limit the range of problems addressed, often to molecular conduction, whereas if they could be overcome new areas would become accessible, such as electrochemistry. Our purpose here is to map out a possible way forward by extending a light weight scattering theory scheme known as Hairy Probes [7] to systems more general than those to which it was originally applied, and to show that simulations can be made computationally efficient and easy to set up. Hairy Probes originally was designed to address time dependent problems; here we only consider the case of steady state current and static atoms.
Using an Empirical Tight Binding (ETB) model we investigate a Cu atomic wire, and then using a Density Functional Tight Binding (DFTB) model [8], we apply the method to the study of a chain of H atoms. These two simple, but well understood systems, allow us to investigate the correctness of the method. For both the Cu and H wires we get ideal ballistic conductance provided the strength of the coupling to the probes is neither too large nor too small: extreme couplings suppress the current. We then look at current flow through a benzene molecule between two graphene contacts as a way to investigate the properties of graphene as a contact for molecular electronics [9,10]. We find that the presence of a covalent bond between the benzene and the graphene suppresses the current as it disrupts the π-conjugation. As preparation for this calculation, current flow through a carbon nanoribbon is studied to find the correct cou-pling strength for the probes to the graphene. We obtain a conductance of slightly less than two quantum units corresponding to two bands crossing the Fermi level.

II. FORMALISM
The Hairy Probes formalism was originally introduced for orthogonal tight binding models, and covered both static and time dependent simulations [7]. Here we generalize the static limit to the non-orthogonal case [11], summarizing the key steps in the theory. The expressions are derived using the Lippmann-Schwinger formalism [12,13], which is equivalent to using non-equilibrium Green's functions (NEGF) for non-interacting or mean field Hamiltonians [14].
We note that this method has a number of similarities with the sink-source potential method [15][16][17]; however, additional simplifications allow for arbitrary bias, any number of terminals, and full self-consistency. Similar simplifications have also been achieved previously by applying the wide band limit directly to the leads [18]. However, we note that Hairy Probes can deliver accurate results for low dimensional systems, and charge self-consistency can be introduced straightforwardly, as demonstrated below.
The starting point is to imagine that our system is connected to one or more particle reservoirs by a set of atomically thin leads (which we call probes) that each attach to just one atomic orbital in the system of interest. The reservoir of electrons from which a probe emerges is characterized by a chemical potential and a temperature for the electrons. Each probe, then, is a bit like a wire attached at one end to a terminal of a battery, and at the other end attached by a kind of alligator clip to an atomic orbital. Each probe thus corresponds to both a source of incoming electrons of given chemical potential and temperature, and a channel for outgoing electrons.
In practice, the probes are attached to contact regions in much the same way that leads are attached to contacts in many NEGF calculations: see Fig. 1. However, because the probes are not system specific, we can define them in a manner that is computationally convenient. Thus there is no need to compute surface Green's functions, the embedding self-energy can be made energy independent while avoiding imposing the wide band limit directly to the leads, and the mean field self-consistent potential profile is taken care of automatically. These attributes are what enable the Hairy Probe formalism to be easy to use (you just need to specify where the probes are to be attached, and how strongly, but do not have to build Green's functions for the leads), and computationally very efficient (an effective Hamiltonian is produced that can be diagonalized, and all subsequent integrals can then be performed analytically). We note that it is shown in [7] that in the limit of long electrodes and small coupling the Hairy Probes steady state reduces to the conventional 2-terminal Landauer picture. The argument we present here that leads to the Hairy Probes equations is based on the Lippmann-Schwinger formulation of scattering theory. That is, we treat each probe as transporting independent electrons from a reservoir, with the electron wavefunctions being viewed as scattering states that travel down the probes and scatter off the system of interest, being partially transmitted (producing a current) and partially reflected.
As we employ ETB and DFTB models, the basis set used to expand the single particle wave functions is composed of atomic orbitals. Let our atomic basis set be denoted by |α where α is a combined index that spans both atomic sites and orbitals. We can now define the Hamiltonian and overlap matrices by H αα = α Ĥ α and S αα = α | α , respectively. Note that the Hamiltonian includes all the terms associated with self-consistent charge redistribution [8,19]. We partition these orbitals between the system (|β ) and the probes (|pγ p ), where p is the index of the probe and γ p indexes an orbital in probe p.
Consider a state with index n p in probe p with energy E pnp that is stationary before the probe is connected to the system. Let us denote this state by is an expansion coefficient. A scattered wave forms from this state after the probe is attached to the system, which we denote by ψ (pnp) = β ψ is an expansion coefficient for orbitals in the system, and ψ (pnp) p γ p is an expansion coefficient for orbitals in probe p . The scattered wave is related to the initial state by the Lippmann-Schwinger equation, giving where G R is the retarded Green's function matrix for the whole system, including all probes, and W β ,pγp (E) = H β,pγp − ES β,pγp is an effective coupling matrix element between the system and probe p.
If we now define the retarded self energy Σ γp,γ p is the retarded Green's function matrix for isolated probe p, we get the following central results where η is a positive infinitesimal, G A β β and Σ (p)A β β are the advanced Green's function and self energy respectively, ρ ββ is the single particle electronic density matrix, and f (p) (E) is the occupancy of the levels inside the isolated probe p, and hence the energy distribution with which electrons are injected into the system by probe p.
We now introduce the Hairy Probe anzatz for the retarded self energy. We note that we want the simplest possible form that still possesses the properties required by a self energy. Making it (almost) energy independent allows us to reduce the problem of finding the scattering states to a simple diagonalization, and by making it local to one orbital we minimize the parameters we have to set. We then end up with the following form where β p is the index of the orbital in the system to which the probe p is attached, and E p,c is the bottom of the band for the electronic states in probe p, taken to be well below any energy levels in the system. As the self energies are imaginary, they have the effect of allowing electrons to be added to, and removed from, the system [15]. The quantities Γ p set the broadening of the states in the system, and define the rates at which electrons can enter or leave. Provided E p,c lies below all levels in the system, then we can substitute Eq. 4 into Eqs. 2 and 3 to give Note that Eq. 6 offers an alternative, albeit unphysical, interpretation of E p,c : it is the energy of the lowest occupied state in probe p, with states with E < E p,c being unoccupied. A discussion of the implications of the choice of E pc is presented in the Appendix C. We write the retarded Green's function as where ζ (r) β and χ (r) β are left and right eigenstates, and (r) the corresponding complex eigenvalue. These satisfy In principle setting χ should satisfy Eq. 9 as the Hamiltonian matrix is symmetric, but we have found that better results are found by solving Eq. 9 explicitly, especially in the presence of degeneracies.
To solve these equations, the numerical procedure we have adopted is as follows. We first transform Eq. 8 from a generalised eigenvalue problem to an ordinary one in the usual way. First we carry out a Cholesky decomposition of the overlap matrix and use the resulting triangular matrices to express the Hamiltonian (including the self-energies) in an orthogonal representation. We then diagonalise the Hamiltonian matrix using a general complex eigensolver as the problem is complex and symmetric, rather than Hermitian, and obtain the complex eigenvalues and right eigenvectors. The left eigenvectors are then obtained by inverting the square matrix of right eigenvectors, and then all eigenvectors are transformed back to the original representation using the triangular matrices from the Cholesky decomposition.
Substituting Eq. 7 into Eq. 6 gives where and can be thought of as a generalized occupancy. We note that in the limit of very weakly coupled leads (Γ p → 0) all having the same coupling strength, the occupancy simplifies to This limiting form of the occupancy matrix is real and diagonal, so the system carries no current, and is a weighted sum of the contributions from each probe. We include it in the spirit of moving open boundaries to problems outside the usual range, as it might be relevant to the case of an electrode in an electrochemical cell. Finally we note that if the populations f (p) ( ) are independent of the probes, then we get back the usual equilibrium expression for the density matrix. We compute the current through the bond between orbitals β and β from where E ββ = rs g rs χ (r) β χ (s) * β and g rs = 1 2π and a factor of 2 for spin degeneracy has been included. A derivation of this expression is given in the Appendix B. We use finite temperature occupations for the electrons in the probes. To enable analytic and efficient evaluation of the integrals involving the occupancies, we use the following piecewise linear approximation to the Fermi-Dirac distribution function: where µ p and T p are the chemical potential and temperature for the electrons in probe p. The integrals are given in the Appendix A.
Finally we note that the transmission between two probes p 1 and p 2 is given by III. RESULTS

A. Atomic wires
The Hairy Probe algorithm has been implemented in the tight binding program Plato [20]. To test the method we first investigated an atomic wire made from 300 Cu atoms; probes were attached to the first 100 and last 100 atoms. We used the orthogonal TB parameterization of Sutton et al. [21] that assigns just one s orbital to each atom. That is, there are 200 probes in all, one per orbital on each of the 200 lead atoms. The probes all have the same coupling strength Γ p , and the same temperature k B T p = 0.001 Ry. Open boundary calculations are carried out in two stages. First, every probe is assigned the same chemical potential, and its value is adjusted until the system as a whole is charge neutral; this we term the reference chemical potential. Each atom individually is allowed to acquire a net charge, described by a monopole with a gaussian charge distribution [22], and charge selfconsistency is imposed. Second, a bias is applied with the chemical potential on the left probes being raised by half the bias relative to the reference chemical potential, and the chemical potential on the right probes being lowered by half the bias. This allows the wire to acquire a net charge, though this is typically less than 1 electron for the whole system for biases up to 3.9 V. The first step is necessary because the probes do not correspond to a known physical system, so an anzatz is needed to give them sensible characteristics.
We computed the current as a function of applied voltage for a range of coupling strengths of the probes; the results are shown in Fig. 2 a). We see that the current varies close to linearly with bias for all probe coupling strengths, and that the slope (conductance) is nearly independent of that coupling for values in the range 0.01 Ry to 0.10 Ry, and in this range the slope is equal to the quantum unit of conductance (G 0 = 7.748 × 10 −5 S). The current is reduced for both larger and smaller couplings. With small couplings the current is restricted by the rate at which charge can be injected and removed by the probes. At very large couplings the hopping matrix elements between atoms in the wire become a weak perturbation on the interaction between the atoms and the probes; in this limit incoming electrons are reflected back into the probes before they can contribute to the current in the wire. The lower and upper bounds for reasonable couplings are roughly the mean spacing between levels (to ensure we have a continuous density of states) and the bandwidth (to ensure the probes do not overwhelm the system).  In Fig. 3 a) is shown the density of states (DOS) projected onto each atom (the atom index is on the x axis) as a function of the electron energy (y axis) for a wire with a bias of 1 V applied, and a coupling of 0.01 Ry for each probe. We see that in the middle of the wire (atom position 150) we have a DOS that is sharply peaked at the band edges. This is consistent with the cosine band structure associated with an infinite chain of atoms with one s orbital per atom. At the ends of the wire there is considerable weight towards the middle of the energy range, consistent with the square root type DOS associated with the end atom of a semi-infinite chain of s orbitals. We note that the states in the lead regions (atoms 1 to 100, and 201 to 300) are significantly broadened by the probes. Finally, the potential in the wire region is essentially independent of position ( Fig. 4 b)). This is to be contrasted with the interface regions where the probes end and begin; here there is a clear variation of potential with position suggesting that this is where the potential drop occurs.
The variation of potential with position can be understood in the following way. The potential in the probe free wire is uniform as it is metallic and the electrons can move to screen out any charge accumulation; current in a perfect conductor requires no field, locally [23]. That leaves the regions with probes. Consider electrons arriving at the left region with probes from the middle region, with energies within the conduction window. In this region, the lifetime of electrons before being absorbed into a probe is τ = /Γ and λ = v group τ is the mean free path, with v group being the group velocity of the electrons at the Fermi energy. For a cosine band with band filling ξ we have v group = 2av sin(πξ)/ , where a is the interatomic spacing, and v is the hopping integral between neighbouring sites. The fraction of electrons that make it to position x (measured from the junction between the perfect wire and the region with probes) dies out as exp(−x/λ). To keep the metal neutral, the bandbottom has to adopt the same shape, to compensate. We thus have the following form for the potential at position x where V is the applied voltage. The functional form of φ(x) clearly has a shape corresponding to that seen in Fig. 3 b). From Eq. 16 we get φ(x)/φ(∞) = 1 − exp(−x/λ). If we let x 1/2 be the point where φ(x 1/2 )/φ(∞) = 1 2 then we get λ = x 1/2 / ln 2. From Fig. 3 b) we see that x 1/2 ≈ 20a and hence λ ≈ 29a. As the hopping integral is v = 0.212 Ry, the band filling is 0.243 [21], and Γ p = 0.01 Ry, we would expect λ/a = 2v sin(πξ)/Γ p ≈ 29; this is in full agreement with the measured value.
We have repeated the above calculations using a nonorthogonal DFTB model for hydrogen [24]: an atomic wire made from 300 H atoms with probes attached to the first 100 and last 100 atoms. The resulting current against bias plot is shown in Fig. 2 b). We see that it has the same structure as for the orthogonal Cu wire (see Fig. 2 a)), and that the maximum conductance is again one quantum unit. This suggests that the method for including overlap into the formalism is correct.
We note that, for the case of orthogonal tight binding, agreement with the two terminal Landauer solution was demonstrated previously for a non-uniform wire, provided a sufficiently large number of probes was employed [7].

B. Graphene contacts
Having studied simple one dimensional atomic wires, and found good agreement with the expected conductance, we now consider electron transport through a more complex system: a benzene ring attached to two graphene contacts by means of covalent bonds. We have a) b) The current through the nanoribbon as a function of the bias voltage for a range of probe coupling strengths (indicated by the symbol G).
selected this system because graphene's electrical properties [9] suggest it might make a good contact material for molecular electronics [10]. As we shall see below, care will have to be taken with how connection to the contacts is made. We note that this system has some similarities to the well studied benzene-dithiol between two gold contacts [25].
To estimate the correct coupling strength of the probes to the graphene contacts we first perform calculations of current through a carbon nanoribbon. To compute the current through a small carbon nanoribbon, whose edges have been terminated with hydrogen atoms (see Fig. 4), we again use a non-orthogonal DFTB model [24]. The probes all have the same coupling strength Γ p , and same temperature k B T p = 0.001 Ry. The variation of current with bias is shown in Fig. 4 for a range of coupling strengths. For coupling strengths of 0.1 Ry and below we find that the current increases roughly linearly with coupling strength for a given bias, and is sensitive to details of the electronic structure of the nanoribbon. The current is fairly insensitive to coupling strength for 0.3 Ry ≤ Γ p ≤ 1 Ry. At large coupling strengths the current is again heavily suppressed. From this we conclude that for carbon flakes of this size, setting Γ p = 0.4 Ry is appropriate. At this coupling, the conductance is 1.44 × 10 −4 S which is 1.9 times the quantum of conductance; this can be understood as resulting from two bands crossing the Fermi energy forming two conductance channels.
Our final simulation is now of the current through a benzene ring coupled covalently to a pair of graphene contacts. The contacts are modelled as small flakes, whose edges are terminated with hydrogen (see Fig. 5 a)). The probes are then attached to the atoms around the edges of each flake, with the probes on one flake all having the same electron chemical potential. The difference between the potentials of the two flakes then creates the bias across the benzene molecule. We use the probe coupling strength of Γ p = 0.4 Ry found from our nanoribbon calculations. Comparing the current versus voltage plot from Fig. 5 b) with that from Fig. 4 b), the first thing to notice is that the current has dropped by a factor of over 1000. This can be understood by looking at the transmission function for the the benzene molecule (Fig  5 c)). Here we see that the reference chemical potential sits well within a tunelling gap several eV wide, thus there are very few free carriers. As the bias increases a small number of holes appear in the valence band; the benzene molecule acquires a small positive charge of order 0.03e, which grows between 1.5V and 4V to about 0.04e. The presence of the band gap is a consequence of the covalent bond between the benzene ring and the graphene: at the point of contact, the carbon atom in the graphene adopts sp 3 hybridization, disrupting the π-conjugation. Thus, to form a good contact, a method is required that maintains the conjugation. Finally, we note that the transport is dominated by holes rather than electrons because the reference chemical potential lies about 0.46 eV closer to the valence band than to the conduction band.

C. Graded probes
Above we have applied the simplest implementation of the Hairy Probes battery, where all probes have the same coupling strength to their respective atoms. This implementation has the conceptual advantage of corresponding most closely to the physical interpretation of the Hairy Probes as external particle baths, in which the system is immersed. In Ref. [7] it was shown that when the length of the hairy leads increases, and Γ decreases (while always remaining larger than the lead energy-level spacing), the Hairy Probes steady state tends to the conventional 2-terminal Landauer steady state.
However in practice one would like to keep the leads as short as possible for computational reasons. The rough rule of thumb for the optimal Γ then is that it should be as small as possible, while remaining larger than the level spacing in the leads. The resultant steady states then approximate the conventional 2-terminal limit, but not exactly. This is not right or wrong as the Hairy Probe battery is intended to be a stand-alone transport setup, with its own interpretation (as above). But the need to consider finite-size effects, and the precise choice of Γ, could then be seen as irksome.
To overcome this complication, a simple alternative is to make Γ position-dependent, so that its value rises gradually from zero, as we move along each lead, away from the central region. We refer to this scenario as Graded Probes. Below we compare these two implementations numerically, and then comment.
The comparison uses the simplest case of a perfect linear atomic chain, with 10-atom long leads with probes and a 4-atom central region without probes. For simplicity we use a single-orbital orthogonal model, with a nearest-neighbour hopping integral set to −1, defining the energy unit. The corresponding energy band then lies in the energy interval (−2, 2), and the 2-terminal Landauer solution has unit transmission throughout that interval.
First we consider the earlier implementation of the Hairy Probes, with a position-independent coupling Γ in each 10-atom long lead. The surface plot in Fig. 6 shows the transmission as a function of energy and Γ.
Consider the limit of small Γ first. In that limit, the 10 + 4 + 10 atom system thinks of itself as a 24-atom linear molecule weakly coupled to an envirnment, which just broadens its 24 molecular states into 24 narrow resonances. This is the origin of the 24 sharp transmission peaks at the small-Γ end of the plot. To understand the opposite limit -large Γ -consider first each 10-atom lead coupled to its probes, but not yet to the central piece (corresponding to the isolated leads in the usual Green'sfunction partitioned approach). If Γ is big enough, it dominates all other energy scales in the lead, ultimately making the lead itself a wide-band system, with a density of states (DOS) going down as 1/Γ. If we now couple the components together, then the 4-atom central region just sees low-DOS adjoining leads, with a correspondingly small embedding self-energy. The upshot is that now the 4-atom central region behaves as a resonant system with weakly broadened states, producing the 4 resonances at the large-Γ end of the plot. In between these two extremes, there is an optimal region of Γ-values, as expected, producing a roughly uniform transmission close to 1, but for the given short leads the corrugation always remains visible. The reason is that even at its optimal value, the finite Γ results in an effective interface (between the regions with and without probes), which -like any interface -generates additional scattering. The longer the leads -and the smaller the optimal Γ -the weaker the disruption.
The Graded Probes provide an alternative way to suppress this boundary scattering, without having to make the leads long. The plot in Fig. 7 shows the Graded Probes transmission, with Γ rising linearly from zero to 1.4 along each 10-atom long lead. It is clear that -at no extra computational expense -we are now much closer to the ideal 2-terminal limit, even for the given modest lead length. The Graded Probes thus provide an alternative, if one wishes to avoid very long leads, or having to consider the precise choice of Γ in the uniform-Γ setup.

IV. CONCLUSIONS
The primary purpose of this paper is to show how to extend the Hairy Probe open boundary method for the steady state to non-orthogonal atomic orbital basis sets. By considering the well understood case of the one dimensional atomic wire (using both orthogonal and nonorthogonal basis sets) we find that we obtain the expected conductance provided the coupling of the wire to the probes has a suitable value. Couplings that are either too large or too small suppress the current: small couplings reduce the rate of charge injection, while large couplings result in high levels of reflection of electrons back into the probes. The optimal value results in the broadening of the system states by the probes to produce a continuous density of states. There is still more work to be done to understand completely the properties of the probes. In addition to studying the dependence of current on applied bias, there are a number of calculations that could be performed, such as the transmission as a function of electron energy for different couplings, or the self-consistent charge distribution.
The method is sufficiently simple that it can be implemented by finding the eigenvalues and eigenvectors of an effective energy independent Hamiltonian, and then performing all the subsequent integrals over energy analytically to produce the single particle density matrix. This results in an efficient algorithm that makes self-consistent open boundary simulations easy to carry out, as it eliminates the need to construct lead self-energies and to perform numerical integrals over energy. The most time consuming part of the calculations is the construction of the density matrices (Eq. 10). For sparse Hamiltonian and overlap matrices, the scaling for building the density matrix is O(N 3 ), which is no worse than the diagonalization step. The absence of numerical integrals also helps keep the prefactor low.
The method was applied to the problem of current flow through a benzene ring attached by covalent bonds to two graphene contacts. It was found that the formation of the contact covalent bonds disrupts the π-conjugation, and thus heavily suppresses the current. We thus conclude that the contacts must either involve physisorption, or a different way to form covalent bonds must be found.
We have also introduced a possible way to accelerate the convergence of the current with respect to lead length by using graded coupling strengths for the probes. The results shown here look very promising, though more work is needed to fully undeerstand them.

ACKNOWLEDGMENTS
We gratefully acknowledge funding from the Leverhulme Trust (RPG-2014-125). M.B. was supported by funding from EOARD (FA8655-12-1-2105) and through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London, funded by EPSRC under Grant No. EP/G036888/1. A.T. was funded by A-star. We also gratefully acknowledge support from the Thomas Young Centre under grant TYC-10. R.D'A. acknowledges support by NAN-OTherm (CSD2010-00044) of the Spanish Ministerio de Economia y Competitividad, and the Grupo Consolidado UPV/EHU del Gobierno Vasco (Grant No. IT578-13). Finally we thank Matas Petreikis for providing improved data for Fig. 4b.
In the limit that |E pc | | |, and for E pc < 0, we have ln sign(E pc ) − |E pc | = ln sign(E pc ) − |E pc | +i arg sign(E pc ) − |E pc | → −iπsign( ) Substituting Eq. C4 into Eq. C3, and noting that sign( (r) ) = −1, we get rs becomes independent of E pc , while I (1) rs varies with E pc as − ln |E pc |, which is independent of r and s. We now recall the expressions for the density matrices The contribution from ∆ρ ββ becomes arbitrarily small for large enough E pc , while the contribution from ∆E ββ is logarithmically divergent. We note that the E pc dependent parts of both matrices are real and symmetric, thus they make no contribution to the electric current, but make a contribution to the atomic forces (not discussed further in this paper). These contributions decrease as the distance from the probes increases, and can be suppressed entirely if we set to zero those overlap matrix elements that link orbitals not attached to probes to those that are attached to probes.
In the main text we offer an alternative interpretation of E pc that allows us to avoid these difficulties, but at the expense of being unphysical: it can be interpreted as the lowest energy for which states in the probes are populated.