A Simple Model of Ballistic Conduction in Multi-Lead Molecular Devices

: A fully analytical model is presented for ballistic conduction in a multi-lead device that is based on a π -conjugated carbon framework attached to a single source lead and several sink leads. This source-and-multiple-sink potential (SMSP) model is rooted in the Ernzerhof source-and-sink potential (SSP) approach and speciﬁes transmission in terms of combinations of structural polynomials based on the molecular graph. The simplicity of the model allows insight into many-lead devices in terms of constituent two-lead devices, description of conduction in the multi-lead device in terms of structural polynomials, molecular orbital channels, and selection rules for active and inert leads and orbitals. In the wide-band limit, transmission can be expressed entirely in terms of characteristic polynomials of vertex-deleted graphs. As limiting cases of maximum connection, complete symmetric devices (CSD) and complete bipartite symmetric devices (CBSD) are deﬁned and solved analytically. These devices have vanishing lead-lead interference effects. Illustrative calculations of transmission curves for model small-molecule systems are presented and selection rules are identiﬁed.


Introduction
The study of molecular conduction and molecule-scale devices is and has been an active area at the boundary between chemistry and physics for at least half a century [1], and by now, it has accumulated a substantial number of textbooks in the literature, e.g., [2][3][4][5][6]. Over this period, different theoretical approaches have emerged. Techniques for calculation have mainly favoured Green's Function (GF) methods, both equilibrium and non-equilibrium GFs [7][8][9][10][11][12][13][14][15][16][17][18][19]. Especially in chemistry, there has been a sustained development of scattering models for the ballistic conduction of electrons though molecular systems, with an emphasis on conduction through channels governed by energies and symmetry characteristics of orbitals [20][21][22][23][24][25][26][27][28][29][30]. An advantage of such pictures is that they lend themselves to even simpler qualitative modelling of the effects of electron energy, relative interaction strengths, contact placement, and qualitative features of molecular electronic structures, as treated in e.g., [31][32][33][34] and as illustrated by many examples in our own work, as cited below.
A special feature of the SSP model in its graph theoretical incarnation [55] is that it gives a useful qualitative account of the two-lead device, giving selection rules for conduction, rationalisation of the sensitivity of current to placing of leads; models for

Two-Lead Devices
In the standard two-lead SSP model, the device consists of a central molecule that is attached by single-atom contacts to the leads (Figure 1). In the steady state, the source lead supports an incoming electron wave that undergoes partial reflection at the interface with the molecule, and the transmitted part of the wave scatters through the molecule and exits into the sink lead. The essential observation behind the model [40] is that the two semi-infinite leads may be replaced by source and sink pseudo-atoms equipped with appropriate complex potentials, which play the roles of delivering the transmitted fraction of the electron to the molecule and removing it downstream, respectively. In the purely graph theoretical version of the model, this step allows replacement of a doubly infinite system by one in which the n × n Hermitian (in fact, real symmetric) molecular Hamiltonian matrix is expanded to an (n + 2) × (n + 2) non-Hermitian matrix. It is convenient to recast this generalised eigenvalue problem in the form of an inhomogeneous linear equation [63]: where P is the SSP matrix, the solution vector c gives the device wavefunction specified by complex entries on the n + 2 vertices, and b is the inhomogeneity term that embodies the boundary conditions. The solution to this equation for each possible value of the energy of the incoming electron yields a compact expression for the zero-voltage transmission as a function of energy, T(E), and as a bonus, gives a useful analysis of current in terms of bond or molecular-orbital (MO) contributions. The transmission function can be cast in terms of the characteristic polynomials of four graphs, G, G − L , G − R , and G − L − R , where L and R are the centres within the graph that are in direct contact with the leads. These polynomials are denoted s, t, u, and v in [55]. Details of the algebra are given in the reference, but the result is that the device transmission has the form [63]: where B(q L , q R ) = (2β L sin q L )(2β R sin q R )β 2 LM β 2

RM
(3) is a band-pass function ensuring that the electron energy is within the conduction bands of the leads, and the denominator is where the polynomials s, t, u, and v are to be defined below in (7), The incoming and outgoing wavevectors q L and q R satisfy the Hückel dispersion relations for the leads: and we take the symmetric case where α L = α R = α, and β L = β R , and refer the zero of all energies to α. The polynomial j 2 satisfies [74]: where the four structural polynomials (SPs), s, t, u, and v, which are each the characteristic polynomial of a graph, are defined with respect to the characteristic matrix E1 − A, where A is the adjacency matrix of the molecular graph: and j can be calculated directly [55,74] as In these equations, superscripts indicate the sets of rows and columns to be deleted from the characteristic matrix.
The transmission T(E) is therefore a ratio of polynomial functions of E, taking values between zero and one and depending implicitly on the molecular graph and the placement of connections. Questions such as transparency, opacity, selection rules, and linkage of transmission to MO channels can all be explored either analytically or numerically using this function [55,[60][61][62][63][64][65][66][67][68][69].
In going to the more general device with multiple leads, the conceptual framework remains intact (Section 3), but the notation needs to be overhauled. The first point of difference concerns the naming of leads. We replace the L and R nomenclature for leads with numerical labelling. For the construction of derivations and proofs, it is convenient to adopt the convention that the source is lead 1, and the sink leads have labels p = 2 to l (l ≤ n). Centres in the molecule are then numbered 1 to n such that, for the first l centres, p in the molecule is connected to lead p, and the numbers l + 1 to n are then used for any centres that are not connected to leads. Formulas for transmission use this convention so that T 1→M describes the transmission from the source to the molecule and T M→p is the transmission from the molecule to sink lead, p. Contiguous numbering of this type is advantageous for proofs that use the block structure of the characteristic matrix of bipartite molecular graphs, for example. (In this paper, we consider only devices where all leads are connected to distinct vertices of the molecular graphs. Degenerate 'ipso' [60] devices would need a modified numbering convention.) A second adjustment is to the nomenclature for SPs. For larger numbers of leads, the labelling of every SP with a different letter becomes unwieldy, so we adopt a single symbol that can be extended to arbitrarily large cases. It is convenient to define where the SPs of increasing order are now defined as cofactors (cf. Equation (8)) of the characteristic matrix and where the indices indicate the rows (left of the comma) and columns (right of the comma) to be struck out. We also chose to scale each SP by multiplying by the connection parameters between each deleted vertex and its connecting lead, a feature that makes the ensuing derivations simpler. Diagonal SPs, those for which the two deletion lists are equal, are now scaled versions of the four structural polynomials from the earlier two-lead papers, e.g., We can also define some related quantities by dividing by s, i.e., k pq···r, st···u = k pq···r, st···u /s. (11) These 'hatted' SPs are important in two respects. First, they emerge naturally in the solution of the SMSP equations, as we shall see in the next section. Second, they are intimately connected with the molecular GF, which in the Hückel tight-binding formalism is the inverse An extension of a theorem by Jacobi (cf. Equation (12), p. 774, [75]) establishes that all our hatted polynomials are determinants: It is also convenient to use an extra piece of notation for the traces of principal minors of SPs. We define the traced quantities, where the sums are over the set of molecular vertices attached to leads. The transmission for the two-lead device, expressed in this revised formalism, is The reader may object that this expression could diverge when the value of E is equal to a molecule eigenvalue, Mi . The divergences in (15) are easily removed. The zero-rank polynomial, s, is expressed in terms of molecular eigenvalues through and if we multiply numerator and denominator of (15) by s 2 , we obtain which with the equivalences k (1) = tβ 2 LM + uβ 2 RM and k (2) = vβ 2 LM β 2 RM is identical with the expression for T(E) in [55].
Expression (17) has no divergences at molecular eigenvalues and has been used in our previous work [60,63,70] to derive selection rules for transmission. The numerical computation of SSP transmission does not use expressions such as (15) directly, relying instead upon a standard linear numerical solution of the SSP matrix equations that does not suffer from numerical instabilities. The purpose of equations such as (15) is for analysis to reveal the physics involved in the transmission process. Our aim in this paper is to generalise these expressions for multi-lead devices.

Multi-Lead Devices
We consider a device comprising a molecule with n atoms and a set of l leads ( Figure 2). The leads are simple 'unstructured' infinite atomic chains described in a Hückel tightbinding model in which a given lead, p, has edge weights β L p and has a simple connection to atom p (using an appropriate molecular numbering scheme) in the molecule with edge weight β LM p . Without loss of generality, the first lead is the single source lead, and the remaining (l − 1) leads are sinks. Each lead has a dispersion relation E = 2β L p cos q L p (18) that relates the energy, E, of a stream of electrons in the lead to the parameter β L p and the associated wavevector q L p (0 ≤ q L p ≤ 2π) in an infinite wire. The description of the molecule may, in principle, go beyond the Hückel model in the manner we have described previously [70], but, for simplicity, here, we derive the Hückel version of the theory.

The SMSP Equations
The Source-Multiple-Sink potential (SMSP) equations are modelled on our previous work with SSP [63], so that we can write them in the form (1) as where the l-dimensional (diagonal) lead block of the S(M)SP matrix, P, is Each nonzero element in P L is the matrix element for the terminal atom p of the inverse Green's function of the semi-infinite linear chain, g −1 L . The molecule block is the characteristic matrix constructed from the molecular adjacency matrix, A M , The diagonal matrix, M , is the matrix of n molecular-orbital eigenvalues, and U is the associated n-dimensional matrix of molecular-orbital eigenvectors. The unitarity of the eigenvectors implies that The lead-molecule block, P LM = P † LM , is an l × n-dimensional matrix with elements On the right-hand side of Equation (19), the inhomogeneity vector, b L , has a single non-zero element which corresponds to the fact that lead 1 is the single source lead. The vector b L ensures that the system of equations satisfies the correct boundary conditions for a semi-infinite source lead with current travelling in the forward direction with pseudo-momentum, q L 1 . The factor is the normalising factor (cf. [63]), ensuring that the forward current is equal to unity. The SMSP solution vector c in Equation (19) is the (complex) device wavefunction. Consequently, we can write the transmission from the molecule to lead p using the standard formula for the electron current [2] as This formula has the correct sign for transmission into sink leads (p = 1), but the sign must of course be reversed to obtain the transmission from the source lead (p = 1).

The Formal Solution to the SMSP Equations
The SMSP equations in Equation (19) can be written out in terms of the separate lead and molecule blocks as We can solve for c M in terms of c L using Equation (21): In terms of components We can now substitute our formal solution for the molecular coefficients into the lead equation to give where the reduced lead matrix, P L , is Equation (30) can be solved by inversion to give where the cofactor matrix and the determinant are defined as (D L ) pq = cof (P L ) pq , D L = det P L (33) and from Equation (29) c Substitution of Equations (32) and (34) into (26) gives where we moved real quantities to the left of the imaginary part operator. Summing Equation (35) over all leads, p, and using the molecular GF for real energies gives since the terms under the summation on the right-hand side are all real, thus verifying Kirchhoff's law for the conservation of current in the device. From the definition of the reduced lead matrix in Equation (31), We can now use the Laplace expansion of the determinant D L to show that For a sink lead (p = 1), we deduce that It follows from Kirchhoff's conservation rule that Equations (35) and (40) give the general SMSP expressions for the currents that flow through the single-source molecular device. These expressions can be used directly for the calculation of ballistic currents in general SMSP devices. However, some far-reaching simplifications are available when we restrict consideration to the usual case where a device consists of a molecule connected to chemically identical leads, i.e., the symmetric device.

Symmetric Multi-Lead Devices
The leads in symmetric devices have identical lead parameters, β L , and identical connection parameters, β LM . This enables a relatively easy simplification of the transmission formula (39) into a cluster expansion, involving pairs, triples, etc. of leads. All leads satisfy a single dispersion relation and furthermore, the reduced lead matrix of Equation (31) is now where the bar over the matrix U indicates that we require rows of the matrix only for those atoms connected to leads. The rectangular matrix U is obviously not unitary and indeed gives different products but the first of these and the diagonal nature of the first term in the l-dimensional matrix P L , allow us to write where the diagonal n-dimensional matrix, Q M , which involves all eigenstates, is Equation (44) expresses the reduced lead matrix, P L , as a triple product involving rectangular matrices. The Cauchy-Binet theorem [76,77] can be used to expand the determinant of a product of two rectangular matrices. Hence, we obtain a further reduction in the quantities in Equation (33) as a sum over a set of products of l-dimensional principal minors where K = {k 1 , k 2 , · · · , k l } is an ordered subset of 1, 2, · · · , n, with |K| = l, and 1 ≤ k 1 < k 2 · · · < k l ≤ n. By deriving this equation, we are able to separate the roles of eigenvectors, U, and the state energies contributing to Q M . Further manipulation of Equation (46) uses algebra familiar in the reduction of matrix elements over Slater determinants in quantum chemistry. We can replace the restricted summation over the ordered sets, K, by unrestricted summations from 1 to n over the state indices (i.e., Hückel MO or GF pole indices), k 1 , · · · , k l . Equalities between the k i indices are allowed because they give rise to zero contributions in Equation (46), since they produce determinants, det U K , with identical columns. We can also use the antisymmetry of the determinants to replace the first determinant by a simple product. Hence, we can conclude that The next step is to expand the product of the (Q M ) k p k p using Equation (45) to generate a cluster expansion in powers of β L and β 2 LM . Whilst it is easy to pursue this generally, it is probably more enlightening to consider the example of a three-lead case explicitly.

The Symmetric Three-Lead Device
In the case of the three-lead device, the denominator (47) is The central triple product can be expanded exactly in the form where r and s are in the set K = {i, j, k}, and we have defined energy denominators Our aim is to obtain an expression in terms of generalised SPs. We proceed by substitution of Equation (49) into Equation (48) and then by examining the resultant coefficient of each power of β L e −iq L . The first term in (49) generates a coefficient where π is a permutation of parity σ π belonging to the permutational group S 3 of the lead labels 1, 2 and 3, and we used the orthonormality of the U matrix rows.
Using the same methodology, the β 2 L e −2iq L term in Equation (49) gives three contributions to the coefficient, the first of which is Likewise, the other two denominators give contributions (g M ) 22 and (g M ) 33 . Reduction in the remaining terms proceeds in a similar manner, but the increased number of denominators is analogous to going from one-electron to two-or three-electron matrix elements over Slater determinants, and it is necessary to include more of the permutational symmetry of the included determinant. Collecting terms, we reach an explicit formula for D L 33 ), (53) or simply where we have used the Jacobi relation from Equation (13) and the generalised SP notation in Equations (11) and (14). We derived this theory in terms of the Hückel tight-binding approximation, but it also holds for correlated molecular GFs, where the MO coefficients, U, in the definitions of SPs are replaced by the GF Dyson orbital coefficients [70]. Furthermore, the cluster expansion in Equation (53) holds for devices with any number of leads. The alternating sign expansion contains terms with determinants of GF matrix elements with increasing dimension up to l, the total number of leads in the device. Now that we have the denominator for the sink lead transmission, we also need to derive an expression for the cofactor matrix element in Equation (39). As an example, we take the 2, 1 matrix element for the three-lead case. We have the same structure as before, but lead 2 matrix elements are omitted from the left-hand part of the expression and lead 1 matrix elements are omitted from the right so that the determinants are of order two: We write out the determinant in full. The simplification of this expression is carried out exactly as before. The first term, β 2 L e −2iq L , in the expansion of the Q factors, however, gives no contribution because the off-diagonal nature of the cofactor creates matches only between orthogonal U matrix elements. The final result for the cofactor is The difference between the cluster expansion of the determinant in Equation (53) and the cofactor in Equation (55) is significant. The determinant expression has traces over principal minors. These satisfy the interlacing theorems we previously utilised in our work on two-lead devices to deduce selection rules [60].
It is straightforward to extend (53) and (56) to more leads. For example, the equivalent formulas for four-lead SMSP devices are as follows:

The Wide-Band Limit
The wide-band limit (WBL) is a commonly used approximation for calculation of electron transport [41]. It is assumed that the wave-vector of the ballistic electron can be taken to be q = π/2, and all other quantities in the expression for the transmission tend to their values for E = 0. In the SSP model of the symmetric two-lead device, the wide-band limiting transmission, T WBL , for non-ipso devices reduces to a formula in terms of the tail coefficients of the characteristic polynomials of the four subgraphs induced by deletion of 0, 1, and 2 distinct connection vertices of the molecular graph G [55]. This expression has been used to develop selection rules for Fermi-level conduction of two-lead devices [60,61].
A similar approach can be used for many-lead devices. As an example, the WBL version of the sink transmission formula for lead 2 in the symmetric three-lead device is This can be expressed entirely in terms of the characteristic polynomials of the seven subgraphs induced by deletion of 0, 1, 2, and 3 distinct connection vertices as follows: whereβ L = β 2 LM /β L , and we understand the lower-case symbols s, t i , v ij , w ijk as the Fermi limits of the characteristic polynomials of the graphs G, and the upper-case symbols T, V, W as their traces. This expression reduces to the corresponding two-lead equation (Equation (21) in [55]) when all SPs involving lead 3 are set to zero. It is straightforward to extend Equation (59) to arbitrary numbers of leads. As the seven polynomials in the three-lead case are linked by an identity for the principal minors of the 2 × 2 × 2 hyperdeterminant (see, e.g., Equation (2) in [78]), T WBL can be calculated from all seven or just six of them. In either case, the WBL approximation can be applied to this instance of a many-lead device to give the transmission from a knowledge of the graph itself and its connections. Only the tail coefficients of the characteristic polynomials are required; this will have particular significance for future development of multi-case selection rules based on the numbers of zero roots of the various polynomials.

Analysis of Interference in Multi-Lead Transmission
In order to appreciate the nature of the transmission for multi-lead devices, we need to examine the numerator and denominator quantities |D Lp1 | 2 and |D L | 2 . Again taking the three-lead case as the example, we have The first term in the expansion involves the same rank 1 off-diagonal SP that appears in the two-lead device constructed from leads 1 and p, albeit with a different |D L | 2 denominator. We shall refer to this as the direct term for the constituent device (1, p). The remaining two terms represent the effect of the third lead (i.e., q) on the conduction through the second lead (p). We shall call these mixed interference and pure interference terms, respectively. Each term has a different dependence on the lead parameter, β L , and the lead-molecule parameter, β LM (cf. Equation (9)). The direct term varies as β 2 L β 4 LM , whilst mixed and pure interference components vary as β L β 6 LM and β 8 LM , respectively. The different scaling implies that the relative importance of direct and interference effects may change markedly for different choices of the relative sizes of these parameters.
The denominator term in the transmission is expressed using our trace notation as Equation (61) is derived for the symmetric three-lead device with β LM1 = β LM2 = β LM3 = β LM , and so the two-lead device cannot be recovered simply by taking β LM3 → 0, but in the limit thatk (3) → 0, Equation (61) tends to a product of a factor β 2 L multiplied by the denominator in Equation (15). The transmission formula for sink lead p in multi-lead symmetric devices has the pre-factor N 2 1 N 2 p = (4β 2 L − E 2 ), which gives the quadratic rise and fall in transmission as a function of energy near the lower and upper band edges, as in the two-lead device.
The analysis of interference terms can be extended to the case where the device has more leads. In the case of l leads, the leading term in the numerator is of the form It has the same SP as the cognate constituent two-lead device, with a dependence on β 4 LM , and this can be expected to be the most important term in the transmission spectrum. The other terms in |(D L ) p1 | 2 represent sink-lead interference terms. They form a cluster expansion with the inclusion of higher rank SPs up to the number of leads, l, and traces over all sink leads other than p. An analysis in terms of a hierarchy of constituent two-lead, three-lead, etc. devices could be envisaged.
The equations for symmetric devices, exemplified by (60) and (61) can be used to derive general features of transmission spectra, selection rules for Fermi transmission, and systematic trends in interference effects. Even further simplifications are possible with additional assumptions about the choice of leads, in particular when the number of leads is large.

Symmetric Devices with Complete Structure
There are two limiting cases of symmetric multi-lead devices that have enough structure to give special algebraic properties in the Hückel tight-binding approximation. The complete symmetric device (CSD) is a device with each atom of the central molecule attached to a lead, with all leads being characterised by the same lead-parameter, β L , and the same connection-parameter, β LM . A bipartite molecule has a graph with vertices (atoms) that can be divided into two sets, such that there are no edges (bonds) between vertices belonging to the same set. A complete bipartite symmetric device (CBSD) is a bipartite molecule where all the vertices of one set connect to distinct leads, each lead having identical β L and β LM parameters, whilst vertices from the other set are unconnected to leads. Both CSDs and CBSDs turn out to have interesting properties with respect to the conduction and (the lack of) interference between leads.

The Complete Symmetric Device
We begin our derivation with Equation (35), written in the form The CSD is distinguished by the fact that Equation (44) can be written using the complete matrix U, so that we can exploit its unitarity to express the reduced lead matrix inverse as Using the definition of the molecular Green's Function and the orthonormality of the MOs, we find that We now observe that the eigenfunctions of a symmetric matrix can always be expressed in terms of real quantities and that we can multiply the numerator and denominator by (Q M ) ii (Q M ) * jj to limit to the numerator the operation of taking the imaginary part: From the definition of Q M in Equation (45), and hence The final step is to convert the unrestricted sums over the indices i and j into restricted sums with i < j. This is performed by taking notice of which parts of the sum are symmetric in the indices (those in square brackets) and those (in round brackets) which are not. It is also convenient to define a denominator with the singularity removed: The final expression for the sink transmission in a CSD is Using Equation (40) and orthonormality, This last formula shows the implicit dependence of the current transmitted through the molecule on the position of the source lead, on the role of the attached vertex in all the molecular orbitals, and on the shell degeneracies.
A notable feature of Equations (70) and (71) is that the coefficients, U, appear only in the numerator, and therefore, the only energy dependence in the numerators arises from the band-pass factor N 4 1 . The denominators are all finite and positive within the band window, and therefore, the transmission in the device vanishes only at the band edges. As later examples show, however, individual sink transmissions may still be small at particular energies.
The orbital sum in Equations (70) and (71) is actually a sum over shells by virtue of the energy difference factor. Non-zero contributions to the total transmission in (71) arise if and only if the source connection vertex carries non-zero density in both i and j. Such pairs exist for all molecular graphs G, as all vertices have a non-zero entry in the LOMO (the eigenvector with the largest eigenvalue, the Perron eigenvector), and orthonormality requires that every vertex has a non-zero contribution to at least one other shell. Equation (70) is more restrictive in that a non-zero contribution for lead p requires both 1 and p to have a non-zero scalar product of entries over both shells of the pair.
Since the numerator for the transmission to sink p in Equation (70) makes no mention of any other sink lead and the denominator is symmetric in all leads, the predicted transmission for a CSD is free of interference terms: the current in lead p has exactly the same numerator as for the current in the constituent device, and differs only in the denominator, where all leads are democratically treated. Thus, the other sink leads influence the current in a given lead, but only in this mean-field sense.

Analysis of the Complete Symmetric Device
In our previous work on two-lead devices, we were able to describe 11 specific cases for conduction with leads attached to distinct vertices [60,63,70]. In particular, we identified inert MOs (or Dyson orbitals for correlated molecular calculations) that had no effect on the transmission at any electron energy, E. On the other hand, active orbitals-more correctly active shells, allowing for possible degeneracy-each produce a maximum, or a shoulder, in the transmission spectrum. This feature is produced by the presence of a pole in the complex energy plane in the expression for the transmission. The width of the feature associated with the pole depends upon the size of the imaginary part of the pole associated with the eigenstate. Inert states do not have poles in the transmission formula, and hence, they produce no such features.
In that treatment, we were also able to identify states i where the transmission vanished at E = Mi . Such pointwise insulation effects depend upon non-cancellation of factors ∆ i between numerators and denominators in the expressions for the transmission, Equation (39). The separation into the 11 cases depended solely on the properties (specifically, the numbers of zero roots) of the SPs in the expression for the two-lead transmission. In the multi-lead case, the quantities D L and (D L ) p1 contain many more SPs of higher rank. The number of different cases allowed by the Interlacing Theorem would be very large. For the present, therefore, we restrict ourselves to a more limited analysis.
For CSD it is straightforward to identify two inertness classes. Strong inertness implies that a specific MO (or Dyson orbital) does not affect the transmission in any lead at any energy, E. From Equation (70), we can see that strong inertness for pole i obtains when the vertex in contact with the source lead 1 is at a node in all orbitals of the shell i, as the denominator L i does not then appear in the transmission formula. We can also have inertness restricted to a specific sink lead, p. This we call weak inertness, and it occurs when the connection vertex p is at a node for all orbitals of shell i: As noted earlier, the transmission is strictly positive across the energy range allowed by the band-pass factor. The maxima in the CSD transmission spectra are close to the poles of the device GF, which in this case are the zeroes of the denominators,

The Complete Bipartite Symmetric Device
We consider a bipartite molecule with n 1 vertices in one partite set (the black set), and n 2 in the other (the white set). If we number the black and white sets contiguously, the adjacency matrix takes the form The (n 1 × n 2 )-dimensional matrix, B, can be reduced to diagonal form using the singular value decomposition where the (n 1 × n 2 ) dimensional matrix, Σ, is with σ 1 ≥ σ 2 ≥ · · · ≥ σ r > 0 = σ r+1 = · · · = σ min(n 1 ,n 2 ) .
Here, r = Rank(B), and V and W, are square n 1 -and n 2 -dimensional unitary matrices. This approach leads to the well-known Coulson-Rushbrooke Pairing Theorem [79,80], which states that the eigenvalues of a bipartite system are paired, with values ±σ i . The eigenvectors corresponding to those eigenvalues are specified mixtures of the V and W. In the case where n 1 = n 2 , there is also a set of supernumerary zero eigenvalues to make up the total number, n 1 + n 2 , for the molecule, and the eigenvectors for these lie entirely in one or other of the unitary spaces, depending on which of n 1 or n 2 is the larger.
In the present context, we assume, without loss of generality, that the black set of vertices are all connected to identical leads with the source being lead 1, but we make no assumption about the relative sizes of the two sets. The SVD decomposition shown in Equation (76), the unitarity of the V and W matrices, and the identical nature of the leads enable us to write the SMSP matrix as by analogy with our CSD treatment in the last subsection. We can now introduce transformed equations in terms of transformed quantities where c M1 , and c M2 refer to the black and white sets of atoms, respectively. The SMSP matrix is composed of a set of diagonal blocks, which implies that the equations separate into 3-dimensional matrix equations one for each singular value, σ i . The solutions are where The original quantities, c L and c M in the atomic orbital basis can be obtained by back transformation, and then, the remainder of the derivation for the transmission follows that of the CSD subsection.
The final expression for the CBSD sink transmission is The denominators are derived from the determinant of the 3 × 3 matrix in (81). The expression should be compared with Equation (69). We can also use Equation (84), with sign reversed, for the source transmission The argument used above for CSD to show that there are no lead-lead interference effects also applies here to CBSD.

Analysis of the Complete Bipartite Symmetric Device
In contrast to the CSD, the CBSD has a pre-factor of E 2 in Equation (84). If uncancelled, this would imply insulation at the Fermi level. Cancellation can occur for graphs with one or more zero eigenvalues. For a non-zero contribution to total transmission, this requires that the source connection vertex is not a node for every orbital of the non-bonding shell. Conditions for individual lead transmissions involve both source and sink having non-zero coefficient products in the non-bonding shell. In these circumstances, the device is a conductor at the Fermi level. Two cases of conducting CBSDs are those based on 4N-cycles, and on odd chains where the leads are attached to the larger partite set. We see examples of these in Section 6 below.

Internal Conduction Channels
We now turn to the question of how electrons are conducted through the central molecular part of the device. We considered this in detail in our previous work [63] in terms of two different models. The simpler of these models treats current as the resultant of the set of atom-to-atom bond currents (see Section V B of [63]). Whilst this approach is chemically appealing, it is closely associated to the Hückel tight-binding Hamiltonian. A more general approach is to use the states of the (isolated) molecule as a vehicle for understanding the conduction process. In the context of the tight-binding approximation, the wavefunctions for these internal states are the MOs obtained from the diagonalisation the Hückel molecular adjacency matrix. However, we also showed in [70] that the Dyson orbitals (DOs) coming from a correlated GF treatment of the molecular ionisation and attachments states can be understood to play the same role with respect to hole and particle conduction through the SMSP device. In this light, the Hückel MOs are just the DOs from the uncorrelated GF defined with respect to the Hückel Hamiltonian. DO currents offer a suitable SMSP model for understanding conduction through molecular devices that is not limited to low-level semi-empirical Hamiltonians and which can range from a one-electron treatment, through Hartree-Fock, to a sophisticated correlated GF method.
We first consider the case of l non-equivalent leads and then proceed to expressions for the special cases of CSD and CBSD devices.

The General Case
The SMSP equations, Equation (19), were defined using the AO basis for the molecular block, P M . An alternative choice is to use the MO vectors, U, to describe the molecule in terms of its internal states, so that in the MO representation, the molecule block is diagonal. In this representation, the connection block is transformed to so that each lead, p, connects with all the MOs through the (p, i) matrix elements of the coefficient matrix, U. The SMSP solution vector is unchanged in the lead block, but the molecular block becomes The inhomogeneity vector is unchanged because its non-zero component is within the lead block.
The SMSP equations can be set up in the molecular basis and solved in the usual manner, using the linear solution method on the SMSP equations, Equation (1). Having obtained the solution vector, c MO M , we can compute the transmission from a given component, c MO Mi , to lead p, using the current formula [2], With these definitions, it is possible to carry out a derivation that parallels the treatment in Section 3.2, to obtain We can see immediately that so that individual MO currents satisfy a sum rule. It is evident that we have achieved a formula for internal channel currents, which is equivalent to a Mulliken-type current analysis of the transmission expression given by Equation (35).

The CSD and CBSD Cases
The derivation of the CSD transmission starts with Equation (63), which is exactly the form in which we can use the current analysis, since one merely removes the sum over the index i present in the g M matrix element. The derivation follows exactly the same steps as in Section 4.1, except that we have a fixed index, i, and a summation over j. The final expression is The lack of a double summation implies that the conversion of the denominator (Q M ) jj to L j as shown in in Equation (69), leaves an uncancelled factor, ∆ j .
The derivation of the CBSD internal currents follows exactly the same pattern. The resulting expression is As they must, the CSD and CBSD internal currents satisfy the sum rules

Symmetry Considerations
The presence of molecular symmetry has significant consequences for the transmission. We investigate this using the three-lead symmetrical device in the tight-binding approximation. Using the cluster expansion of D L in Equation (53), we can explore the consequences of symmetry by constraining the connection vertices for leads 2 and 3 to be equivalent. This can be achieved by insisting that the molecular GF has elements in the transmission formula, Equation (39). This is the direct term in Equation (60), so that the mixed and pure interference terms in the numerator have vanished. The denominator, |D L | 2 , on the other hand, is a modified version of the two-lead analogue, with three-lead terms giving shifts in pole positions that alter the positions of the features. Symmetric three-lead devices, nevertheless, have sink-lead transmissions that are similar to those of the corresponding two-lead device. Furthermore, leads 2 and 3 in this particular example have identical transmission spectra, with each carrying exactly half the transmission at any energy, E.
We now consider what happens when all three connection vertices in a three lead device are equivalent. Then, 33 , and (g M ) 12 = (g M ) 13 = (g M ) 23 . (100) Hence, and (D L ) p1 factors in the identical manner to Equation (98), because its underlying properties reflect the permutational symmetry of the sink vertices. Once again, this ensures vanishing of interference terms. The cluster expansion form of D L , on the other hand, does not distinguish between source and sink vertices and reflects higher permutational symmetry. Devices with more leads exhibit even richer symmetry effects.

Results
The SMSP formalism described in Section 3 has been implemented in our suite of routines that use the Maple 2019 package [81] to provide analytical and numerical calculations of the zero-voltage transmission curves T(E) for multi-lead systems.
The Hückel tight-binding calculations that are used for the central molecule use a single 2p z -orbital basis function on each atom with a hopping parameter, β, to represent the interaction between π-bonded carbon centres. Hartree-Fock (HF), and second-order Green's function (GF2) calculations require two-electron ππ-interactions. We used a Hamiltonian with a Hubbard single-centre interaction parameter, U = 1β, for all centres and a single two-centre interaction parameter, W = 0.5β, between atoms that are π-bonded. The methodology is identical to that in our earlier work [70] for the two-lead SSP model, where it is described in more detail. No two-electron interactions were used to describe the source and sink atoms that represent the leads in the SMSP method. For most calculations, we used β L = 2β for the lead hopping parameter and β LM = 1β for the lead molecule connection parameter. All energies are shown in units of the negative quantity β.
The labelling schemes for the molecules under study are shown in Figure 3. Figures 4-7 show transmission curves plotted against electron energy. In each case, the black outer envelope is the total transmission from the source lead into the molecule, and coloured curves show (symmetry distinct) transmissions from the molecule into sink leads. A colour code is used to distinguish the leads in clockwise order from the source. The plots also indicate poles of the molecular Green's function, using an encoding of blue for attachment and red for ionisation poles.

Devices Based on the Benzene Ring
As a first example, we take the six-membered ring of carbon atoms, which can be connected as a two-lead device in three ways ('ortho', 'meta', 'para') and as a three-lead device in six ways, leading to the nine transmission curves shown in Figure 4. All curves are symmetric about the Fermi energy in the Hückel model, as the six-ring is a bipartite graph. Transmission curves for the three two-lead devices are qualitatively different. The ortho and para devices have non-zero transmission at the Fermi energy, with different patterns in the wings, as predicted by the analytical expression for T(E) of a cycle [55] and selection rules [61,63] for conduction at degenerate and non-degenerate eigenvalues (here, α ± β). In contrast, the meta device is insulating at the Fermi level.
Plots for the various three-lead devices show a variety of patterns, but it is striking that the overall shapes of the curves can be interpreted in terms of constituent two-lead components: to a first approximation, the curve T(E) for the three-lead device with source 1, and sinks a and b has lead contributions that follow the patterns for devices (1, a) and (1, b) in the region of the Fermi level but are damped in the higher and lower energy wings. The implication is that the interference terms (see Equation (60)) have smooth, predictable effects on T(E) for the composite device. This pattern is especially clear for the symmetric (1, 3, 5) device (Figure 4i), where the source is in a meta relations to both sinks, and the Fermi-level transmission is zero. In fact, this device is a CBSD: interference terms between leads do not appear in the numerators of the lead currents, and the device is interference-free. In other devices such as (1, 2, 3) and (1, 3, 4), a vanishing Fermi-level contribution from the meta constituent device is masked by the contribution from the other constituent. The counting of ortho, meta and para constituent devices works well as a rough guide to T(E), though we may expect this simplicity to be diluted in larger systems and devices with more leads.
The final row of Figure 4 shows the transmission curves calculated at the GF2 level for (1, 2), (1, 4) constituent devices, flanking the plot for the composite three-lead device (1, 2, 4). The plots show strong resemblances to those calculated at the Hückel level, and interpretation of the lead contributions in terms of constituent devices survives intact in the correlated treatment.

Devices Based on Anthracene
This section presents results on devices based on anthracene, which, similar to benzene, has a bipartite molecular graph. The devices considered here all have three leads, with the source on atom 4 in the numbering of Figure 3. The devices for which transmission curves are shown in Figure 5a-c are (4,8,14), (4,9,13) and (4,8,11), which span the types (•, •, •), (•, •, •), and (•, •, •), where the colours of the circles represent the partite sets ( Figure 3b). All three necessarily have single-vertex deleted and triple-vertex deleted graphs that are non-Kekulean and have at least one non-Kekulean constituent two-lead device. Devices (4,8,14) and (4,9,13) are C 2v symmetric, and their total tramsission is simply twice the contribution of each sink lead. The second row of the Figure shows the three two-lead constituent devices (4,8), (4,9), and (4,11). In each case, there is a clear qualitative correspondence between the three-lead device and its constituent devices, i.e., (a) with (d); (b) with (e); and in the mirror-symmetric (c), the distinct lead contributions strongly resemble (f) and (d), generally as expected. Exceptions to this agreement are in the features where the transmission falls sharply to zero at eigenvalues ±2 and ±1 in the two-lead (4,8)and (4,9) constituent devices; these are suppressed in the mirror symmetric three-wire devices as an effect of specific cancellations with the denominator.
The final row of panels Figure 5g-i refers to interference terms in lead currents in the three-lead devices. Explicit calculation of interference contributions to the numerators of the expressions for the lead currents shows that the mixed and pure interference contributions are identically zero for the mirror-symmetric devices (4,8,14) and (4,9,13) but not for the non-symmetric one (4,8,11). As shown in Figure 5i, interference terms of both mixed and pure types appear. It is risky to generalise too much from one example, but it is at least interesting to see that these terms exert little effect in the central region of the spectrum (where they are both small) but are more influential in the wings (where they are more intense and have regions of opposite sign). Figure 6 shows a cascade of calculated curves T(E) for non-bipartite (i.e., non-alternant) pentalene. We take one example of a four-lead device and compare it with all of the constituent two-and three-lead devices. The chosen four-lead case has connections at vertices 1,3,5,7 (Figure 6a), which are the four core vertices for the non-degenerate nonbonding LUMO of pentalene (i.e., they carry the non-zero charge/spin densities for this orbital).

A Non-Bipartite Case: Devices Based on Pentalene
The plots are no longer symmetric, as the spectrum of a non-bipartite graph is not paired, but there are still discernible qualitative features of the transmission curve that are associated with eigenvalues: broad or compound maxima at most eigenvalues, and a sharp spike in the transmission at the HOMO eigenvalue (α + 0.47068β in Hückel theory). The plot also shows a broad dip in transmission at antibonding energies above the LUMO. The origins of the lead contributions to devices in the middle row are plausibly traced from the appropriate two-lead devices ((e) and (f) to (b); (e) and (g) to (d); (f) and (g) to (d)), and the constituent device curves for (e), (f), and (g) appear in the plot for the full device (a), with some blurring of features by interference effects. The family relationships of devices with progressively increasing numbers of leads are yet again helpful as a rough guide to the features of the transmission spectrum of a complex device.  Figure 7 shows calculated curves T(E) for molecular devices of CSD and CBSD types based on cycles and linear polyenes. The first three panels show the transmission curves for CSD based on 6-, 8-, and 10-membered rings (Figure 7a-c). They show qualitatively similar behaviour, with broad maxima, shoulders, and minima associated with the eigenvalues of the respective graphs. The pattern of a central minimum for (4N + 2)-rings and a central maximum for 4N-rings persists to larger ring sizes. Traces of the constituent devices are apparent in the lead contributions, e.g., compare ortho, meta, and para benzene devices (Figure 4a-c) with the CSD in Figure 7a. As noted in Section 4.1, transmission is non-zero across the whole band window. However, the summed current can include some very low sink-lead contributions. In the case of the 10-ring, for example, the lead transmissions at the Fermi level for leads 2, 3, 4, 5, and 6 are 0.1581, 0.0445, 0.0853, 0.0042, and 0.0674, respectively. The next three panels (Figure 7d-f) show the transmission curves for CBSD based on 8-, 10-, and 12-membered rings (Figure 7a-c). Again, a pattern of alternating behaviour at the Fermi level is apparent and is confirmed by extension of the calculations to larger rings: CBSD based on 4N-rings show a central maximum in total transmission, consistent with superposition of sink-lead curves with central maxima, whereas (4N + 2)-rings have an insulating minimum at the Fermi level, indicating vanishing of all lead contributions. As noted in Section 4.3, this is part of a general pattern based on the number of nonbonding orbitals (NBOs), or in other words, the nullity of the graph. In the 4N-ring, the nullity is 2 and all centres carry non-zero density arising from any occupation of the nonbonding shell (i.e., a 4N-ring is a core graph [82]). Cancellation of the E 2 factor between numerator and denominator in the expression for lead currents therefore yields conduction at the Fermi level for these rings. The (4N + 2)-rings have no such cancellation. This convex/concave pattern for antiaromatic/aromatic rings is one example of a nullity-based selection rule for a family of devices. There will be many more to be found.

Complete Devices of Types CSD and CBSD
Finally, the third row (Figure 7g-i) deals with complete bipartite devices based on the heptratriene linear polyene. As shown in Figure 3, the molecular graph has two partite sets of vertices. The graph has nullity one and the unique Hückel NBO has non-zero entries only on the larger (black) partite set. This leads to a qualitative distinction between devices of CBSD type based on the two sets. For the (1,3,5,7) device, there is non-zero density in the NBO at all lead positions, and cancellation of the E 2 factor in the numerator leads to Fermi conduction. On the other hand, the two devices based on the smaller (white) partite set show insulation at the Fermi level, even though the underlying graph is still of course singular, as all lead positions are at nodes in the NBO.
A comparison of Figure 7h-i shows that this central feature of the transmission spectrum is independent of the choice of source lead position within the same partite set. There is, however, a significant difference between the (2,4,6)-and (4,2,6)-devices. The molecule has a mirror plane through the central vertex, so that MOs that are antisymmetric with respect to the mirror have U i4 = 0. The (2,4,6)-device in panel (h) displays weak inertness for lead 2 attached to vertex 4 (in red) but not for lead 3 attached to vertex 6 (the green curve). Hence, no maximum is visible in the red curve near the second eigenvalue from the right or the second from the left. The (4,2,6)-device, in panel (i), exhibits strong inertness for these same eigenvalues because the source lead is attached to vertex 4. These observations are again indicative of the wealth of underlying selection rules.

Connections with the Meir-Wingreen Formula
The Meir-Wingreen (MW) formula gives an exact expression for the current through a molecular device described using a correlated molecule but with leads and lead-molecule interactions treated in the one-electron approximation. In previous work [70], we showed that the SSP method is consistent with the MW formula, in the elastic-scattering limit when the electron-phonon interactions are neglected (i.e., in the Born-Oppenheimer approximation). SSP produces identical formulae for the transmission, including the case in which the description of the central molecule includes electron correlation.
The transmission in the SMSP approach is determined by boundary conditions on a device wave function with conceptual advantages in the chemical interpretability of the internal molecular channels for conduction. These channels are determined by the attached and ionised states appearing in the Lehmann representation of the equilibrium molecular GF. The properties of these channels are defined in terms of the characteristics of the DOs associated with each of the GF poles. The sets of orbitals and poles define spectral expansions of structural polynomials used to formulate the expressions we derived in this paper. The relevant formal properties of the polynomials are retained regardless of the level of theory, from Hückel, through HF, to sophisticated GF formalisms. The molecular part of the problem may therefore be treated with empirical, semi-empirical, or ab initio methods, with or without the inclusion of correlation, thus giving a seamless single formalism to describe conduction through molecular devices.
At any of these levels of approximation, the inverse (P) −1 of the SSP/SMSP device matrix is the device GF (Equation (85) of [70]). The poles of this matrix are complex-valued and represent device resonances corresponding to the presence of an extra electron (or hole) in the sea of (possibly correlated) electrons. The imaginary parts of these poles are inversely related to the lifetimes of the resonances. These lifetimes represent transit times for the particles or holes crossing the device, times that are particularly short for Koopmans poles but long for shake poles. In the SMSP context, it can be seen that the poles of the device GF can be obtained by solving the equation |D L (E)| 2 = 0. Each complex conjugate pair of poles is then associated with a Dyson orbital and (real) pole energy of the molecular GF. Hence, the SMSP model gives a direct route to calculation and interpretation of all the properties that can be accessed through GF methods.

Conclusions
Working equations have been derived for an extended version of the SSP model for calculation of ballistic currents flowing through molecular conductors under a potential difference. They constitute the new Source-and-Multiple-Sink Potential model. The new model describes systems with an incoming electron travelling along a single source lead that drains out through the molecule into multiple sinks. As with its SSP predecessor, SMSP lends itself to a graph-theoretical approach, which leads to a formulation that is compatible with, but not confined to the Hückel/tight-binding description of electronic structure of π systems. SMSP therefore retains the advantages of SSP. It yields closed-form expressions for device transmission as a function of electron energy: for general lead parameters (Section 3.2); for chemically symmetric leads (Section 3.3); and for leads with coverage of all π-centres (Section 4.1), or of all starred or unstarred centres of an alternate π-system (Section 4.3).
As its SSP predecessor did, the SMSP model gives predictions for qualitative selection rules and generic patterns of conduction for families of molecular systems, such as the opposite behaviour of aromatic and anti-aromatic rings at the Fermi level (Section 6.4). In addition to a calculus for the prediction of total transmission with energy, the model allows for an analysis of molecular current at various levels of detail: by bond current, by orbital-based molecular channels, and by destination lead (Section 5). Patterns of current in multi-lead devices can be interpreted as built up from notional constituent two-lead devices, where SMSP gives a natural breakdown into direct contributions, interference terms (again governed by selection rules), and mean-field corrections (Section 3.5). The wide-band limit yields particularly simple expressions for transmission involving only the characteristic polynomials of the set of vertex-deleted graphs (Section 3.4). In the hypothetical limit of full or half coverage of a molecule by leads, the specific interference effects are predicted to vanish (Sections 4.1 and 4.3). The symmetry of the device can also cause the vanishing of interference effects (Section 5.3) when the connection point for the sink leads are equivalent.
SMSP, similar to SSP, is a parameterised model, with two parameters β L , and β LM (defined in units of β, the resonance integral for the molecule). Adjusting the parameters affects the detailed appearance of the transmission spectrum and the relative importance of interference effects. The effective lead parameter β L is likely to be larger in magnitude than β, even in all-carbon devices since the co-ordination number is higher with structured leads. The molecular energy levels, therefore, lie inside the band as assumed in our sample calculations. The role of β LM is to describe the perturbation of the molecule by the leads, the effect of which is to move poles of the molecular GF off the real axis, and hence, lower β LM values give sharper transmission peaks. Interference effects become more important as the ratio β 2 LM /β L increases, as shown in (Section 3.5). Perhaps the main advantage of models such as SMSP is that many interpretative features persist when more realistic models of electronic structure are used. In particular, the analysis by molecular-orbital-based channels survives the introduction of electron-electron interaction into the theory.
Thus far, we are not aware of experimental work on molecular systems with dense attachment of leads to single molecules. Mechanically controlled break-junction (MCB) techniques [83] provide IV characteristics that give insight into averaged conduction for two-lead molecular devices rather than their detailed transmission spectra. Even so, an impressive range of data on the factors influencing molecular conductivity can be inferred from painstaking experimentation based on careful molecular design. For example, encouraging agreement with experimentally observed regiospecificity of conduction [84] can already be achieved using non-empirical and graph-theoretical SSP calculations [43,60]. The long-term aim of the present work is to carry forward this modelling of trends to multi-lead devices, in the hope that it may be of use in the design of future experiments on these intriguing systems.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: