On the resilience of magic number theory for conductance ratios of aromatic molecules

If simple guidelines could be established for understanding how quantum interference (QI) can be exploited to control the flow of electricity through single molecules, then new functional molecules, which exploit room-temperature QI could be rapidly identified and subsequently screened. Recently it was demonstrated that conductance ratios of molecules with aromatic cores, with different connectivities to electrodes, can be predicted using a simple and easy-to-use “magic number theory.” In contrast with counting rules and “curly-arrow” descriptions of destructive QI, magic number theory captures the many forms of constructive QI, which can occur in molecular cores. Here we address the question of how conductance ratios are affected by electron-electron interactions. We find that due to cancellations of opposing trends, when Coulomb interactions and screening due to electrodes are switched on, conductance ratios are rather resilient. Consequently, qualitative trends in conductance ratios of molecules with extended pi systems can be predicted using simple ‘non-interacting’ magic number tables, without the need for large-scale computations. On the other hand, for certain connectivities, deviations from non-interacting conductance ratios can be significant and therefore such connectivities are of interest for probing the interplay between Coulomb interactions, connectivity and QI in single-molecule electron transport.


Supplementary Note 1: The model
The system consists of a molecular core, coupled to conducting leads. Non-interacting part of the core is described by the tight binding model Here i and j run over the sites of the molecule, i.e. the p z orbitals centred on each of the carbon atoms. c † i,s and c i,s are the electron creation and annihilation operators for the orbital centred on site i and with spin s. n i,s = c † i,s c i,s is the electron number operator with n i = s n i,s . ε i is the energy of the orbital relative to that at the chiral symmetric point. γ ij are hopping integrals. Taking into account that the next nearest neighbour hoppings in graphene are at least an order of magnitude smaller than the nearest neighbour ones [3], in what follows we retain only the nearest neighbour hopping integrals that we set to γ = 2.4eV.
When the Coulomb electron-electron interaction is taken into account, we include it according to the Parr-Pariser-Pople (PPP) model [1,2] and the whole Hamiltonian is of the It consists of the on-site interaction U ii and the long range interaction U ij . For the latter we 1 use the Ohno interpolation [4] U ij = λ where U = λU 0 is an interaction strength, U = U 0 = 11.13eV giving the physical value of the interaction. d ij is the distance between sites i and j, for nearest-neighbour sites is equal to d 0 = 1.41Å.
The leads are modelled as chains of atoms on sites i connected by nearest neighbour hopping integrals γ 0 and are included in the system with a term − i,s,α γ 0 c † α,i+1,s c α,i,s +H.c.. Here α ∈ {L, R} labels the left (source) and the right (drain) lead, respectively, and i ≥ 1 runs over all atomic sites of one lead. c † α,i,s and c α,i,s are the electron creation and annihilation operator for lead sites, respectively. Eigenstates of an infinite lead are plane waves with wave vector k and eigenenergy E(k) = −2γ 0 cos k. The coupling between lead and the molecule is described by − α,s V c † α,1,s c iα,s + H.c.. Here V is the hopping integral between the lead site closest to the molecule and the molecular site i α to which lead α is attached. We take a wide-band limit, γ 0 = 10γ and V = γ. The appropriate quantity that describes strength of coupling between the molecule and the lead is the spectral width This is the relevant coupling parameter since it includes the density of states of the lead, ρ ∝ 1/γ 0 , and the probability for an electron to jump between the lead and the molecule, V 2 .
It is constant in the energy interval of interest, Γ α = 2V 2 γ 0 = 0.2γ and small, justifying the weak coupling limit. The retarded Green's function of the whole system can be expressed in the elastic cotunneling approximation we adopt for the Lanczos calculation which has a simple form where G(E) is the Green's function for the isolated molecule and the influence of the leads is included via Σ α , the retarded self-energy. Its value due to coupling with lead α is Recently it was shown experimentally [5] that molecular levels shift as a result of electron interaction with image charges in the metal leads, resulting in a HOMO-LUMO gap renormalization. We take the image charge effects into account by analytically solving [6] the Poisson's 2 equation for the electrostatic Green's function in a simplified geometry, namely we assume the leads are two infinite parallel plates. The renormalized interaction values are Here r i = (x i , y i , z i ) is a vector pointing to site i and L is the distance between the leads. It depends on the distance between the connectivities and on d -the distance between the lead and the site, which it is connected to. d is measured in units of d 0 -the lattice constant.

Supplementary Note 2: The Hartree-Fock method
The Hartree-Fock method (HF) [7] is an approximation in which the interacting term in the PPP Hamiltonian is evaluated to the first order, leading to an effective HF Hamiltonian (8) s denotes anti-parallel spin polarization of s. In the restricted Hartree-Fock (HF) approximation the first two terms are zero because n i,↑ = n i,↓ = 1 2 which is due to chiral symmetry that is not broken in the HF approximation (see Supplementary Note 3). The Hamiltonian in the HF approximation simplifies to the tight binding Hamiltonian with effective hoppings New long-range hoppings are introduced for i and j on different sublattices. The conductance can again be calculated with the Landauer-Büttiker formula [8,9] where T (0) is the transmitivity at the gap centre E F = 0 as calculated from the HF Hamiltonian.

Supplementary Note 3: Hartree-Fock approximation preserves chiral symmetry
For bipartite lattice models the chiral symmetry [17] is defined with the operator S which acts on a site creation and annihilation operator as where α ∈ {0, 1} is the sublattice index. A system has a symmetry whenever the equality The HF approximation transforms the PPP model to the tight binding model with long range hopping amplitudes (9) and the chiral symmetry is therefore preserved if c † iα,s c jα,s = 0. The following calculation proves that HF preserves chiral symmetry: The HF Hamiltonian is calculated iteratively and self-consistently and the first iteration starts with the tight binding Hamiltonian, which is chiral symmetric because on the leads and on the molecule only nearest neighbour sites are connected. If one proves that c † iα,s c jα,s = 0 for a chiral symmetric system, then after the first iteration no inter-lattice hoppings are introduced and the Hamiltonian remains chiral symmetric with c † iα,s c jα,s = 0. The system therefore stays chiral symmetric for all later iterations proving that the HF Hamiltonian is chiral symmetric. Single particle tight binding Hamiltonian of the whole system has the block off-diagonal form where the connectivity matrix C contains hopping amplitudes γ ij between sites on different sublattices. A non-interacting system has chiral symmetry if there exists an unitary operator U S that transforms the Hamiltonian as U † S HU S = −H [17]. For this system the corresponding symmetry operator is the Pauli operator s z acting on different sublattice spaces. The single particle eigenstate |ψ k,s with energy eigenvalue E k can be decomposed into two components, one on sublattice α = 0 and another on α = 1, where ψ k,i,s = ψ k,i,−s = ψ k,i and |α represents a state on the sublattice with index α and |0 the vacuum state. In chiral symmetric systems, every |ψ k(E),s with eigenenergy E has a partner s z |ψ k(E),s = |ψ k(−E),s with eigenvalue −E and eigenvector Here we assumed the there is a HOMO-LUMO gap for a half-filled system, so there is no state at E F = 0. The expectation value c † iα,s c jα,s is calculated for a half-filled system in the ground state of the HF Hamiltonian, that has all of the eigenstates with E k < E F filled with For a completely full system this expectation value is zero i.e., Here we used ψ k(E),iα ψ * k(E),jα = ψ k(−E),iα ψ * k(−E),jα that is evident from equation (13). This proves that c † iα,s c jα,s = 0 for a chiral symmetric Hamiltonian, which further proves that HF Hamiltonian stays chiral symmetric. From here also follows that the single-particle HF Hamiltonian is again of off-diagonal block form as the non-interacting one in equation (11) where C HF is the Hartree-Fock connectivity matrix filled with effective hopping amplitudes from equation (9), C HF ij = γ HF ij .

Supplementary Note 4: Validity of the Hartree-Fock method
Lanczos diagonalization produces exact results in the limit of weak coupling between the molecule and metallic leads but processing time grows exponentially with number of 5 atomic sites making the calculation impractical for larger molecules. Calculation with the Hartree-Fock method is much faster, growing polynomially with size. Although it is an approximate mean-field method and is known to accurately describe systems with sufficiently weak Coulomb interactions. It is therefore prudent to test its validity for the interacting PPP model of polyaromatic molecules, for which the Coulomb interaction might be significant.
Hartree-Fock results for transmitivity, density of states and energy gap are compared with results from the Lanczos method. Fig. 3 Supplementary Fig. 15. It is evident that staggered magnetization becomes non-zero at around λ ≈ 1.5. This phenomena occurs also in graphene when described by the Hartree-Fock PPP model [7]. Phase transition happens at slightly lower λ than for anthanthrene, which is consistent with the observation that the larger the molecule is the lower the phase transition point below which the system is paramagnetic.

Supplementary Note 5: Infinite range interaction limit
The PPP model can be considered as a model that is intermediate between two extreme cases, a system with localized interaction and a system with infinite range interaction. The first one is described by the Hubbard model that has U ii = U and U i,j =i = 0. As seen from equation (8) in the HF approximation of the Hubbard model the interaction has no effect.
The model with infinite range interaction has U ij =Ũ for all pairs of i and j whereŨ is the 6 average value of the PPP interaction integrals in a given molecule. In this section the M -table in the infinite range interaction limit is derived. Since the PPP interaction is somewhere in between both limiting cases (Hubbard and infinite-range), its M -table will be somewhere in between M -tables of those two limit cases.
The non-interacting Hamiltonian H of an isolated molecule, that is a bipartite lattice with N 2 sites per sublattice, can be expressed in terms of a N 2 × N 2 connectivity matrix C as in equation (11) and is symmetric under the chiral symmetry operator s z . When wave functions are expressed in terms of wave functions on each sublattice, see equation (12), the Schrödinger equation is Since the system has a gap and no gap states, E k = 0 and the eigenstate sublattice wave functions are related, Therefore, only the problem on one sublattice needs to be solved, where C † C is Hermitian and and its eigenvectors form a unitary matrix, Since the probabilities to find an electron on each sublattice are equal: the normalized eigenstates with E k > 0 (E k < 0) are therefore Let us form a N × N 2 matrix V containing occupied (E k < 0) eigenstates of H in its columns: 7 In terms of V the single-particle correlations are In particular, for j and i in the same sublattice c † j,s c i,s = 1 2 δ ji . For conductance ratios in the limit of weak coupling between the molecule and leads we consider the Green's function at the centre of the HOMO-LUMO gap which is expressed in terms of the M-table which can be calculated from the connectivity matrix asM Note that the M -table as defined here does not contain integers in general, but the non-integer factors are cancelled out when calculating conductance ratios.
The PPP Hamiltonian in the HF approximation leads to an effective Hamiltonian H HF that has chiral symmetry so it can also be expressed in the form of equation (11), that is in terms of the HF connectivity matrix C HF containing effective hoppings. The HF conductance ratios are then given by the HF M-table (In the main text M int is an M table for interaction system for general method of calculation, here we explicitly use the name of HF approximation.) Assuming the interaction is independent of distance (infinite range interaction, U ji =Ũ ), the HF self-consistency equation reads, see Eqs. (8) and (25): 8 Its solution is Note that the HF Hamiltonian see Eq. (17) and (19), has the same eigenvectors as the non-interacting one, with the HOMO-LUMO gap increased by the value of the interactionŨ being the only difference in the eigenvalue spectrum. The infinite range HF Hamiltonian can be expressed with non-interacting one as sgnA is a matrix one gets if one diagonalizes the matrix A = P ΛP −1 and then replaces eigenvalues Λ ii by their sign [sgnΛ] ii = Λ ii /|Λ ii | and then rotates the matrix back to original basis of A, sgnA = P −1 sgnΛP . Equation (33) follows from the fact that H HF has the offdiagonal form as in (11) and the fact that E k /|E k | are eigenvalues of sgnH.
The exact PPP Hamiltonian in the infinite range interaction limit can be written as 9 The exact expression for Green's function of the core for the infinite range interaction is (see ForŨ γ the M -table ratios become independent of the interaction strength and the HF M -table can be rescaled toŨ or, using Eq. (25)Ũ In this limit, all occupied states merge into the HOMO level at E = −Ũ 2 and all the empty states merge into the LUMO level at E = U 2 . In order to compare the results obtained using infinite range interaction and the PPP parametrisation we show, in Supplementary Fig. 16(a), correlations of the Hartree-Fock conductance ratio for a particular pair of connectivities (horizontal axis) with the non-interacting (blue dots) and the infinite range interaction (orange dots) conductance ratio for the same pair of connectivities. Results for all possible pairs of connectivities are shown. The first observation is in remarkable agreement between infinite range and the PPP model results. In is also clearly demonstrated that non-interacting results (i.e., magic ratios) are correlated with the PPP results on average, but exhibiting noticeable deviations in some cases. Maximum deviations of PPP results compared to non-interaction ratios are limited within the scaling range 0.1 ∼ 10. In Supplementary Fig. 16(b) the same type of analysis is shown but in the limit of large interaction strength, U → ∞. Here the correlation between infinite range and PPP model is more dispersed, but still remarkable.
In Supplementary Fig. 16 In Supplementary Fig. 17 are presented pyrene results corresponding to Supplementary   Fig. 16. Qualitatively all conclusions are unchanged.

Supplementary Note 6: Systematic view on the change of conductance ratios
By comparing magic ratios with conductance ratios for a system with Coulomb interaction, with or without screening, one can analyse the changes and try to find a general description for deviations from magic ratios.

Supplementary Note 6.1: Effects of interaction with no screening
For systems with Coulomb interaction and no screening there is a qualitative rule which describes deviations of conductance ratios from magic ratios. To illustrate this correlation graphs for anthanthrene are shown in Supplementary Fig. 18. Correlation between two quantities is defined as This behaviour can be explained by observing the molecule energy spectra at different interaction strengths, shown in Supplementary Fig. 19 and Supplementary Fig. 20. In the limit of strong interaction, energy levels above and below the HOMO-LUMO gap come much closer to each other, as if they were squeezed to a very narrow band. Energy levels are renormalized so the HOMO-LUMO gap stays constant with interaction. By definition, elements of the M -table are not affected by such renormalization factors. In the limit of strong interaction one can therefore make an approximation for the eigenenergy E k of the eigenstate |ψ k,s : where ψ k,i,s = ψ k,i,−s = ψ k,i are wave function coefficients of the eigenstates of HF  Supplementary Fig. 19 and Supplementary Fig. 20 the ψ k,i ψ * k,j are presented at every level k at different U and it is evident that this quantity, directly connected to wave functions, does not change with interaction. The same observation is analytically derived in Supplementary Note 5 for a model with infinite range interaction, from which its follows that this expectation values can be evaluated in the non-interacting ground state. In the limit of strong interaction one can therefore qualitatively estimate conduction ratios according to a simple non-interacting expression As seen from equation (38), the approximation becomes exact for the infinite range interaction model with infinite interaction strength. The fact that above relation qualitatively reproduces conductance ratios can be read from Table 1 in the main text. M HF ij are therefore proportional to s c † i,s c j,s and since it is known [16] that this quantity decreases with distance between site i and site j this explains why M HF ij also decreases from magic integers if i and j are far apart.

Supplementary Note 6.2: Effects of lead screening
In some cases screening does not affect the conductance ratios whilst in others they are changed drastically. As its name implies screening usually decreases effective interaction strength. However, the width of the HOMO-LUMO gap is changed differently for different lead connectivities because screening depends on the distance between the two leads and the distance depends on connectivity. The consequence is that the calculated M-tables are different for different connectivities, for example [M HF 1,2 ] ij = [M HF 3,7 ] ij , in first case the leads are connected to sites 1 and 2 and in second case to 3 and 7.
Change of conductance ratios can be seen as a combined effect of different HOMO-LUMO gap renormalization and different M-tables. This is evident by rewriting conductance ratios in the basis of energy eigenstates |ψ k,s ij and |ψ k,s lm , which are different for different connectivities. Eigenenergies E ij k and coefficients ψ ij k,i correspond to the system with connectivity i − j while E lm k and ψ lm k,l correspond to the system with connectivity l − m.
where E ij g is the HOMO-LUMO gap of a system with connectivity i − j and M-table ratio in case of screening as In Supplementary table 2

Supplementary Note 7: Electron currents
An illustrative representation of conduction processes is to plot electron currents through the molecule in terms of so called bond currents [14]. The bond current I ij is defined as the current that flows along a bond connecting site i and site j. At T = 0 K and infinitesimal difference of lead chemical potentials µ L − µ R = eV and V → 0, it is equal to the expectation value of j ij , the current operator between sites i and j in the single-particle scattering state at the Fermi energy from left (source) lead state |ϕ k F ,s,L = j ϕ k F ,j c † j,s,L |0 . The final expression that was used to calculate bond currents is Currents can be expressed in terms of Green's function for a real Hamiltonian as while the total current through the junction is equal to The current distribution through the bonds is In Supplementary Fig. 21(b) and (c) , bond currents I ij are plotted as arrows between sites. The arrow direction is that of the current and its thickness the magnitude of the bond current. In    The comparison of transmitivity of energy from HF (red) and Lanczos (black) calculation for naphthalene at connectivity 3 − 8 with λ = U/U 0 = 1 and with screening for different lead distances is shown. d is the distance between the lead and nearest atom on the molecule. Supplementary Figure 19: Energy levels of naphthalene Energy levels (grey lines) in naphthalene are shown in dependence to interaction strength λ. They are renormalized so that the HOMO-LUMO gap stays constant. Coloured circles represent ψ k,i ψ * k,j , radius is proportional to magnitude and red/blue colours denote +/− sign.