Ab initio complex band structure of conjugated polymers: Effects of hydrid DFT and GW schemes

The non-resonant tunneling regime for charge transfer across nanojunctions is critically dependent on the so-called \beta{} parameter, governing the exponential decay of the current as the length of the junction increases. For periodic materials, this parameter can be theoretically evaluated by computing the complex band structure (CBS) -- or evanescent states -- of the material forming the tunneling junction. In this work we present the calculation of the CBS for organic polymers using a variety of computational schemes, including standard local, semilocal, and hybrid-exchange density functionals, and many-body perturbation theory within the GW approximation. We compare the description of localization and \beta{} parameters among the adopted methods and with experimental data. We show that local and semilocal density functionals systematically underestimate the \beta{} parameter, while hybrid-exchange schemes partially correct for this discrepancy, resulting in a much better agreement with GW calculations and experiments. Self-consistency effects and self-energy representation issues of the GW corrections are discussed together with the use of Wannier functions to interpolate the electronic band-structure.

The non-resonant tunneling regime for charge transfer across nanojunctions is critically dependent on the so-called β parameter, governing the exponential decay of the current as the length of the junction increases. For periodic materials, this parameter can be theoretically evaluated by computing the complex band structure (CBS) -or evanescent states -of the material forming the tunneling junction. In this work we present the calculation of the CBS for organic polymers using a variety of computational schemes, including standard local, semilocal, and hybrid-exchange density functionals, and many-body perturbation theory within the GW approximation. We compare the description of localization and β parameters among the adopted methods and with experimental data. We show that local and semilocal density functionals systematically underestimate the β parameter, while hybrid-exchange schemes partially correct for this discrepancy, resulting in a much better agreement with GW calculations and experiments. Self-consistency effects and self-energy representation issues of the GW corrections are discussed together with the use of Wannier functions to interpolate the electronic band-structure.

I. INTRODUCTION
The fields of molecular electronics and charge transport through nanojunctions have been deeply investigated in the past fifteen years. 1-3 At the experimental level many different techniques have been developed, including those based on break junctions, nanostructured and scanning probe layouts, or self-assembled monolayers. 3,4 Significant improvements in the accuracy with which these junctions are characterized have been achieved over the years, eg to address the I-V characteristics of single molecular junctions. Moreover, a large experimental literature exists 5-7 on non-resonant tunneling experiments, where it is possible to determine the exponential decay (β 0 ) of the current I = I 0 exp(−β 0 L) as a function of the length L of the tunneling layer [such as eg a layer of organic self-assembled monolayer (SAM) connected to metallic electrodes]. Even though the β 0 parameter depends 6,8-10 also on the detailed nature of the interface, it carries mostly information about the properties of the tunneling layer itself, which makes β 0 an interesting analysis and characterization tool.
Experimentally, these measurements are performed using different setups, ranging from metal-insulator-metal (MIM) junctions, as mentioned above, to the evaluation of kinetic constants of electro-transfer reactions (optically or electrochemically induced) in donor-bridge-acceptor molecular complexes. [11][12][13][14][15] In terms of systems, measurements have been performed on a number of cases ranging from saturated olephins (alkanes) 5,12,13 to biological molecules (such as DNA [16][17][18] ). Theoretically, the electronic mechanism underlying these experiments has been analyzed and understood. 1,5,8,10 As stated in Ref. [8], the key parameter β can be expressed (eg in MIM junctions) in terms of (i) the band gap E g and the (frontier) band widths (or hopping parameter) t of the insulating layer, and (ii) the alignment of the Fermi level in the metals with the energy gap of the insulator. Indeed, the effect of the electronic structure of the insulating layer can be singled out by evaluating 8 the complex band structure (CBS), or evanescent states, in the limit of an infinitely long insulating region. The CBS approach is also particularly interesting for an ab initio evaluation of β, where the calculations can be performed either using wavefunction- 19 or Green's function-based 8 approaches. A cartoon describing the relation between the electronic structure of the MIM junction and the evaluation of the β-decay factor is given in Fig. 1.
Nevertheless, since the β parameter can be 8 directly related to the ratio between the energy gap and the band width of the insulator layer, the accuracy of standard electronic-structure simulations based on the Kohn-Sham (KS) framework of the density functional theory (DFT) can be questioned. In fact, using eigenvalues computed from KS-DFT it is well known that the fundamental band gap is badly underestimated and (when using local and semilocal approximations) the delocalization of electronic states is typically overestimated. Moreover, the description of this class of experiments in terms of single particle energies would require them to be interpreted as quasiparticle energies, in order to address the electronic dynamics of the system. This is valid for advanced MBPT methods, 20 such as Hedin's GW approximation, 21,22 and, at least in a perturbative sense, also for Hartree-Fock calculations. However, DFT Kohn-Sham states are fictitious orbitals with no direct physical interpretation, and their use in this context can only be justified by the assumption that the exchange-correlation kernel is an approximation for the quasi-particle Hamiltonian. Model self-energies have also been recently 23 used to correct the electronic structure of metal-molecule-metal junctions and found to be important for the evaluation of β. To date there has been no systematic investigation of the performance of different ab initio schemes in the calculation of the β decay factors.
In this paper we study the effects of hybrid exchangecorrelation functionals 24 and the GW 21,22 approximation on the calculation of the β decay-factors, according to a scheme based on the the complex band structure (CBS) formalism. 8 A detailed comparison with local and semilocal functionals is also provided. The theoretical background of CBS, and GW and hybrid functionals are described in Sec. II and III respectively. We discuss a general application of the Wannier functions interpolation to the case of GW electronic structure (Sec. III B). Our approach is applied to a number of polymers as reported in Fig. 2: We compute the CBS for poly-ethylene (PE) and poly-acetylene (PA) as references for saturated and conjugated chains respectively. We then consider poly-(para phenylene-vinylene) (PPV) 25 and poly-(phenyleneimine) 7,26 (PPI), which are relevant from a technological point of view, where we compare also with recent experimental data. 7 All chains are studied as isolated, see App. B for full numerical details.

II. METHOD: TRANSPORT
To simulate the decay factors of non-resonant tunneling experiments we adopt the CBS algorithm proposed by Tomfohr and Sankey (TS). 8 For a recent discussion of the connection of transport properties with the CBS theory see also the work by Prodan and Car. 10 Within the TS approach, we need to evaluate the CBS in the limit of an infinitely thick insulating region (β is in fact an asymptotic behavior). The outcome of this procedure is a set of β(E) curves. The value β 0 which has to be compared with the experiments is the smallest one (i.e. the most penetrating) aligned with the Fermi level of the junction: Since the electrodes are not considered in the calculation, together with the proper metal-insulator interface, E F is not known a priori and must be either estimated or calculated separately. This issue is discussed in detail in Ref. [8]. In the present work we will not compute the Fermi level alignment explicitly, and will rely on the estimation method proposed in the above reference, 8 which consists in evaluating β(E) at the energy

FIG. 2:
Polymers studied: (a) PE (poly-ethylene), with each Carbon atom in the polymer chain being fully saturated by two Hydrogen atoms; (b) PA (poly-acetylene), where sp 2 hybridization of Carbon atoms in the chain (each one also bonded to one Hydrogen atom) implies π-conjugation, i.e. alternation of single and double C-C bonds along the chain; (c) PPV [poly-(para phenylene-vinylene)], constituted by benzene rings connected through vinyl groups, also presenting conjugation; (d) PPI [poly-(phenylene-imine)], differing from PPV for the substitution of one Carbon atom in each vinyl group by a N atom, which still preserves the polymer conjugation.
where dβ/dE = 0 (branch point). This method is originally due to Tersoff 27 and based on the pinning of E F by metal induced gap states (MIGS). We discuss this approximation in Sec. V where we compare with experimental results. Note that the value of β at the branch point is also connected (for non-metallic 1D systems with a local potential) with the degree of localization of the density matrix, since β determines 28,29 its spatial decay. Pictorially, a better description of β(E) means then an improved description of the electronic localization.
Scanning the energy spectrum, the CBS procedure searches for evanescent solutions to a given effective single-particle Hamiltonian. By definition, states with (real) Bloch symmetry k satisfy the relation where R is any direct lattice vector,T (R) a translation operator, and λ = e ik·R . In the same way it is possible to define a complex Bloch symmetry κ = k+iβ/2 setting λ = e iκ·R = e ik·R e −β·R/2 , which is thus no longer a pure phase. The imaginary part of κ implies a real space exponential decay of the wavefunctions and it is customary to define β(E) = 2 |Im[κ(E) ·ê]| (ê transport direction). The energy dependence of κ comes from the fact that for a fixed energy E, the solutions are searched in terms of κ, as it is usually done in scattering theory. By adopting a localized basis set {|φ iR } (i and R orbital and lattice indexes, respectively), it is possible to define the Hamiltonian H and overlap S operators through their matrix elements H ij (R) = φ i0 |H|φ jR and S ij (R) = φ i0 |φ jR , and the wavefunctions as |ψ l = iR C il (R) |φ iR . Setting Z = H − E S, the eigenvalue equation for H can be written as: where matrix multiplication is implicit, and N is the number defining the last non-zero matrix Z(R N ). Here we are assuming a real space decay of the Hamiltonian and overlap matrices, which is typically physical even if the range can be strongly dependent on the scheme used to define the Hamiltonian. We will comment later on this point when discussing the use of GW or HF methods. It is then possible to derive 8 the following system of 2N matrix equations: where Eq. (5) is a direct consequence of Eq. (2). Such matrix equations become an eigenvalue problem for λ assuming we can invert the matrix Z(R N ). As pointed out in Ref. [8], this is very often not the case as the matrix is singular, but the singluarities can be avoided. The reader is referred to the original work for the details. The full algorithm proposed in Ref. [8] has been implemented in the WanT code 30,31 and used for the present work. A simple tight-binding analytical model (a generalized version of the one presented in Ref. [8]) is discussed in App. A. This model will be used in Sec. IV to fit and interpret the real and complex band-structures of the polymers we have investigated here.

III. METHOD: ELECTRONIC STRUCTURE
According to the above discussion, in order to simulate the decay coefficient of a MIM junction, we need to compute the electronic structure of the insulating layer (considered as infinitely extended). The underlying reason for this simplification is that we interpret the computed single-particle energies of the system as quasi-particle (QP) energies, which in turn determine the dynamics of singly charged excitations. In general the transport problem for interacting systems is more complicated than that and requires a more sophisticated treatment. [32][33][34][35][36][37][38][39][40] For DFT, the understanding of how the exact KS-DFT Hamiltonian performs to compute transport has been recently the subject of several investigations. [41][42][43][44] Apart from the properties of the exact functional, currently available DFT approximations like LDA or GGA have been demonstrated to systematically overestimate the conductance, especially for off-resonant junctions. 32,45,46 Pragmatically, this suggests that corrections beyond local and semilocal KS-DFT approaches are needed.
Using QP-corrected electronic structure to compute charge transport (eg by means of the Landauer formula 47 ) through interacting systems seems to be a reasonable approximation when finite-lifetime effects are weak. 37,48 Indeed, a number of works computing QP energies by means of the GW approximation 35,[49][50][51][52] or model self-energies 23,53-55 have been reported in the literature. On the other hand, hybrid-exchange functional methods like B3LYP 56 or PBE0 57 are widely used and lead to band gaps and band widths which are usually closer 58 to the experimental values than simple semi-local KS approaches, for both molecules and solids. Recent works [59][60][61][62] have further investigated the accuracy of hybrid exchange functionals, also in comparison with GW calculations. In this work we compare GW and hybrid functionals for the calculation of electronic structure and transport properties of selected organic polymers. In the following we summarize the GW approximation and underline some formal similarities with hybrid functionals.

A. The GW approximation and hybrid DFT
A many-body theoretical formulation of the electronic structure problem can be obtained by using the Green's function formalism. The one-particle excitation energies of an interacting system are the poles of its interacting Green's function G(E), 20,63 which can be written as: where h 0 is an effective single particle Hamiltonian and Σ(E) is the non-local, non-hermitean, frequency dependent self-energy operator. In general, Σ is not known a priori and must be approximated. In this work the selfenergy is computed within the GW approximation: 21,22 where W (ω) is the screened Coulomb interaction evaluated in the random phase approximation (RPA). For more details see e.g. Refs. [63,64]. In the simplest implementation of the GW approximation, the self-energy is computed non self-consistently, i.e. by evaluating G and W according to the eigenvalues and eigenvectors of a reference non-interacting Hamiltonian (typically the KS Hamilotonian at the LDA or GGA level). Such a procedure is known as G 0 W 0 and gives reasonable results for the quasiparticle energies in a number of cases. 63,65,66 In the present paper we exploit the G 0 W 0 approximation, and evaluate the frequency integrals in Eq. (7) by using a plasmon pole model according to Godby and Needs. 67 In this work, the main quantities we are interested in are the QP energies. If we neglect finite lifetime effects and take (or symmetrize) the self-energy to be hermitean, QP energies are given by first order perturbation theory as: As customary, 64 in order to solve for QP m , the self-energy in the above equation is expanded to first order as a function of E.
In order to get more physical insight we also refer to the (static) COHSEX 21,64 approximation of GW, where the self-energy is written as Σ = Σ COH + Σ SEX : Here γ(r 1 , r 2 ) is the one particle density matrix, and W p = W − v is the dynamical contribution to W (being v the bare Coulomb interaction). The COHSEX selfenergy is thus the sum of a statically screened exchange term and a local potential. This partition is particularly useful for discussing the connection between hybrid exchange functionals and GW. In the former case, the potential can be written as: where v NL x is the non-local exchange potential, while v L x and v L c are local potentials. It is then straightforward to interpret α as an inverse effective screening to stress the formal analogy of Eqs. (9-10) and (11). Similar considerations are of course valid for more complex forms of hybrid functionals, like range separated or local formulations. 24,68,69 This formal analogy is well-known in the literature, 24 and it has also been further investigated recently. 70 Besides their accuracy for thermochemistry, this analysis highlights that also the electronic structure computed by non-local hybrids can benefit from the inclusion of some screened exchange term. Indeed, improved description of the electronic structure for finite and extended systems are typically found, 58,70 even though the accuracy may vary significantly depending on the system.

B. Interpolation of GW electronic structure using Wannier functions
Dealing with periodic systems, interpolation over the first Brillouin zone (1BZ) is a long standing issue. Calculations are typically performed by discretizing k-points in the 1BZ, and some post-processing schemes [such as eg those to compute density-of-states (DOS), Fermi surface, band structure, or phonons] might need a better discretization of 1BZ than some of the previous steps (often aimed at computing total energy, forces, and charge density). This is particularly critical when the computational requirements of the adopted methods limit the k-point discretization, as is the case for GW calculations. Schemes able to refine or interpolate 71 over the 1BZ are particularly useful for this purpose. One of these is the Wannier interpolation, 31,[72][73][74] where the localization of the Wannier function (WF) basis together with the finite range in real space of the Hamiltonian are used to perform a Fourier interpolation of the eigenvalues, and eventually eigenvectors. The use of this scheme to interpolate GW results has been also reported elsewhere by Hamann and Vanderbilt. 75 The procedure can be applied not only to the Hamiltonian, but in principle to any operator A(r, r ) with the translational symmetry of the Hamiltonian. First, we define the projector P over a subspace of interest: in terms of the eigenvectors of H. When the subspace P is complete, we can represent A as: Here, A is diagonal with respect to the k-index because it commutes with the translation operators of the direct lattice (as assumed). In practice, limiting the number of eigenstates of H included in P is equivalent to considering the projection of A on the P subspace, namely A P = P AP instead of A. At this point we can use the definition of maximally localized Wannier functions (MLWF's), according to Ref. [73]: to obtain an expression for the matrix elements of A on the Wannier basis, A P ij (R) = w i0 |A|w jR : Note that when the original Marzari-Vanderbilt procedure 73 is applied without any disentanglement, 74 the U k matrices are a unitary mapping of N Bloch states (usually occupied, but not necessarily) into N WF's. Instead, when the disentanglement is performed, the resulting WF set does not span the whole P subspace (U k are then rectangular matrices). This means that in the general case the final representation of A is actually not projected on P but on the smaller subspace spanned by the WFs.
Assuming that A P is decaying fast enough in real space to have ||A P (R)|| 0 for |R| > |R 0 | (where R 0 is within the finite set compatible with the initial k-point grid), we can perform the following Fourier interpolation to obtain the matrix elements of A for any k point: When interpolating GW results, we want to represent the operator Σ(E) = Σ GW (E) − v xc which is in general non-local, non-hermitean and frequency dependent. For the sake of the Wannier interpolation, we are mainly interested to check that the intrinsic non-locality of P ΣP is compatible with the selected k-point grid (or, in other terms, that the GW calculation is converged with respect to the number of k-points used). The localization of the GW self-energy is further discussed in Sec. V A, especially in connection with the usual approximation that neglects off-diagonal Σ mn (kE) matrix elements.

C. Numerical approach
In this work, DFT and hybrid-DFT calculations have been performed using the CRYSTAL09 package. 76 The code implements all-electron electronic structure methods within periodic boundary conditions and adopts an atomic basis set expanded in Gaussian functions (further details in App. B). Once the Hamiltonian matrix elements are obtained, 77 the real and complex band structures are interpolated (as discussed in the previous Sections) using the WanT 30,31 package.
GW results have been obtained using the plane-wave and pseudopotentials implementation of SaX 78 which is interfaced to Quantum-ESPRESSO 79 (QE) for DFT calculations. In this case, once the Kohn-Sham electronic structure is evaluated, we first compute MLWF's 73,74 using WanT and then apply the CBS technique. In order to assess any systematic error in comparing GW and hybrid-DFT results (which have been obtained using different basis sets such as plane waves and local orbitals), we have also performed hybrid-DFT and HF calculations using QE and SaX. In this case WFs are computed on top of the already corrected electronic structure. Results are shown for the case of polyacetylene (see Tab. III and Fig. 8(b) in particular). The excellent agreement between the two sets of data suggests that the pseudopotential approximation, the basis set, and the numerical thresholds are sufficiently well converged to have negligible influence on the results presented. Full computational details and parameters are reported in App B.

A. Structural properties
Before focussing on the electronic and transport properties of the isolated polymer chains from Fig. 2, we investigate their structure by fully relaxing both the atomic positions and the cell parameters using different exchange-correlation schemes. All systems are treated with one-dimensional periodicity, the details of the calculations (performed using CRYSTAL09) have been given in App. B. The results for the lattice parameters are reported in Tab. I.
In the case of PA, electronic properties such as the band gap (as well as the evanescent states) are strongly dependent on the dimerization of the C-C bond lengths (Peierls distortion). Such bond alternation is not easily captured by local and semilocal DFT schemes leading to band gaps that are far too small, i.e. to an overestimation of the metallicity. Since these parameters are critical to our purpose, we report also the bond length alternation (BLA) for PA in Tab. I. Our HF results are in excellent agreement with previously published re-sults. 80,81 In the following we will label as PA HF the calculations performed using the geometry from Ref. [80], where c = 2.469Å and BLA is 1.339/1.451Å. We also decided to consider two frozen geometries of PA (namely PA 1 and PA 2 ) according to other data in the literature. 8,82 In the case of PA 1 82 we set c = 2.451Å and BLA to 1.370/1.460Å, while for PA 2 8 c = 2.496Å and BLA to 1.340/1.540Å. In passing we note that the PA 1 geometry is also very similar to the one adopted in Refs. [83][84][85], where c = 2.451Å and BLA is set to 1.360/1.440Å, according to experimental data. 86,87 Since these theoretical studies report GW results, in Sec. IV B we will compare with PA 1 .
Data from multiple geometries of PA are useful to decouple the electronic and structural effects of the adopted XC schemes on the complex band structure. Since the extent of this goes beyond PA, in the following we decided to look at the effect of XC treatment first using a frozen geometry, independently of the adopted method, and then to compare also with the same quantities obtained using fully relaxed geometries. Moreover, since GW corrections are usually computed without performing a further structural relaxation, working at fixed geometry allows us to compare GW corrections and results from hybrid functionals for identical geometries. For each polymer except PA we adopted the geometry obtained after full relaxation at the PBE level using Quantum-ESPRESSO. Lattice parameters (2.564Å for PE, 6.702 A for PPV and 13.004Å for PPI) are in very good agreement with those in Tab. I obtained using PBE in CRYS-TAL09.
B. Electronic structure: poly-ethylene and poly-acetylene In this section we investigate the effect of different XC schemes (local and semilocal DFT, hybrid functionals, HF, and GW) on the electronic properties of two prototype polymers (PE and PA), while results for PPV and PPI will be reported in Sec. IV C. We compute both the real and the complex band structures. Figures 3-4 refer to the case of frozen geometries (i.e. geometry is not changing according to the adopted scheme, see also Sec. IV A). Details about the electronic structure of the three PA geometries studied are reported in Tab. III. In the case of PA (as well as PPV), we can also compare with previously published GW calculations. [83][84][85]88 Note that the GW results are interpolated using MLWF's, which results in filtering out some of the states above the vacuum level.
In Fig. 3, we report 77 the real and complex band structure for PE. Our results are in reasonably good agreement with previously published theoretical data. 8,10,19,89,90 As can be seen from panels (a,b), the maximum value of β from the arc across the fundamental gap is not changing much when passing from PBE to PBE0 (β max 0.8 A −1 ). In agreement with Ref. [10], we believe this to be one of the reasons why the β computed at the LDA/GGA levels have been found in agreement with the experiment.
The G 0 W 0 -corrected band structure is reported in panel (d). It compares favorably with existing literature data. 19 The main qualitative difference with Fig. 2(a) of Ref. [19] comes from the filtering of some vacuum-like states, due to the use of the localized basis (or WF interpolation) in our approach. Our results are more similar to those presented in Ref. [8], obtained using a localized basis set. This has little influence on the part of the CBS spectrum that is physically relevant for the non-resonant tunneling experiment. The G 0 W 0 CBS is not computed here because the missing self-consistency is critical for CBS, as it is discussed in Sec. V A. Instead, we have computed CBS from the self-consistent 91 COHSEX electronic structure. All results for β max are collected in Tab. II, where we compare data obtained by using CRYSTAL and Quantum-ESPRESSO or SaX after Wannierization. Results and trends compare reasonably well. In Fig. 4, we show the results for real and complex band structure in the case of PA 1 . The dependence of the β parameter on the exchnage-correlation scheme is displayed. exchange and correlation for the β parameter is displayed. As expected, when increasing the percentage of non-local exchange from PBE (0%), to PBE0 (25%) up to HF (100%), the band gap opens, and the "degree of localization" increases as indicated by the increasing values of β max . This trend is general and also found for the other conjugated polymers studied in this work. A detailed description of the computed band structure for PA in the PA HF , PA 1 , and PA 2 geometries is reported in Tab. III for all the methods.
In order to gain more physical insight from our calculations, we have also fitted the data by using the generalized nearest-neighbors (NN) model presented in Sec. II and App. A. This model has three parameters, which approximately correspond to the band gap E g , and the band widths of the HOMO and LUMO bands (related to the parameters t 1 and t 2 ). The proper relation between the band widths and t 1,2 is given in Eqs. (A11-A12).
The main difference between this model and the one introduced in Ref. [8] is that the one used here allows for different band widths for the HOMO and LUMO. While this difference is found to be small (but generally not negligible) in the cases studied, the generalized model allows for a more accurate fitting of the electronic structure of polymers. As can be seen in Fig. 4(a-c), the model fitting (dashed, red lines) is very accurate for all the local, semilocal and hybrid functionals. In general the HF data shows the largest deviation from the model. We believe the non-local nature of the exchange potential to be the origin of this behavior. fundamental gap is 0.8 eV at the PBE KS-DFT level, while it increases to 2.05 eV when G 0 W 0 is applied. This is in very good agreement with previous GW data for PA. [83][84][85]92 As already mentioned, GW calculations are performed by using plane-waves and pseudopotentials, while hybrid-DFT is evaluated on a localized basis set. In order to assess a possible systematic error due to this procedure, we have also computed the electronic structure of PA HF at the LDA, PBE and PBE0 levels, by using the plane-wave implementation of Quantum-ESPRESSO(QE). Results are reported in Tab. III at the lines X-QE (X=LDA,PBE,PBE0). The two sets of data are found to be in excellent agreement, allowing for a direct comparison of GW and hybrid-DFT data. See also Sec. V A and Fig. 8(b) for further details. As a last remark, according to the standard approach for GW calculations, 63,64 only the Kohn-Sham eigenvalues are corrected, without modifying the wavefunctions. In other words, only the diagonal matrix-elements of the self-energy (on the original DFT Bloch eigenvectors) are considered. The effects of this is analyzed in detail in Sec. V A, when it will be demonstrated that this procedure has a sizable effect on the calculation of the complex band structure. For this reason, the GW-CBS directly computed is reported in a lighter color in the right panel of Fig. 4(d), while β max is put in parentheses in Tab. III.

C. Electronic structure: PPV and PPI
In this Section we consider two further polymers, namely poly-para phenylene-vinylene (PPV) and poly phenylene-imine (PPI). PPV has been largely investi-  gated for its role in organic (opto-)electronics, 25 while oligo phenylene-imine molecules attached to gold leads have been recently considered and the β decay coefficents measured. 7,26 This makes these two polymers particularly appealing for our analysis. Results for PPV are reported in Fig. 5 and Tab. IV. The behavior of the real and complex band structures are qualitatively in agreement with what we have found for PA. PBE0 results (β max = 0.26Å −1 ) are those among the hybrid functionals that best compare with the interpolated GW data (β max = 0.28Å −1 ), though slightly underestimating β max and to a larger extent the band gap. In the case of PPI, results are reported in Fig. 6 and Tab. V. For this polymer, the model fitting has been applied on a cell with half length and then bands have been folded leading to four interpolated bands instead of two. Despite the reduction of translational symmetry, this procedure leads to a better fit because it describes a larger part of the frontier electronic structure.
In agreement with previous cases, while we do not have GW results for PPI, we consider PBE0 data (β relax max =0.29 A −1 ) as our best estimate. In the next Section we will also discuss the comparison of these computed data with recent experimental results. 7,26 For the PPV and PPI cases, we have also studied separately the effect of the geometrical relaxation induced by the different functionals on the CBS. Such effect is consistent with the trends already observed at fixed geometry. In the last column of Tabs. IV, V we report the β max value for the relaxed polymer geometries (labelled as β relax max ). The β max parameters increase further with increasing fraction of non-local exchange, when the geometries are relaxed according to the adopted functional. While the coupling of the electronic structure with the structural properties is particularly evident and critical in the case of PA (since the opening of the gap is due to Peierls distortion of the C-C bonds), it is much less pronounced for PPV and PPI where it accounts for a correction term only, the leading contribution being the description of the electronic levels.

A. Analysis of the GW data
In this Section we discuss the effect of some of the approximations involved in the evaluation of the GW selfenergy. In particular we address issues related to the representation of Σ when computing the CBS, as well as the effect of the self-consistency. The GW self-consistency is investigated within the static COHSEX approximation. In doing so, we discuss the localization properties of the resulting Hamiltonians together with the quality of the NN model fitting of the CBS. We focus on the case of PA, which is a good prototype for this class of one-dimensional systems.
Let us begin with the representation problem. As already recalled in Sec. III B, assuming a large enough subset of Bloch vectors (in principles, all of them), the selfenergy operator can be represented as: [see also Eqs. (13,14)]. In the usual GW practice, besides evaluating the self-energy by using the underlying KS-DFT electronic structure for G and W (G 0 W 0 approximation, i.e. no self-consistency), it is also customary to neglect the off-diagonal band indexes m = n in Eq. (18) when computing QP energies. This approximation forces the self-energy to be diagonal on the KS-DFT Bloch states, and thus allows us to modify the quasiparticle energies without changing the DFT wavefunctions. If we assume that the correctly represented selfenergy [Eq. (18)] is physically short-ranged in real space (further comments follow), as it happens eg for HF and COHSEX, when the representation is taken to be diagonal on the Bloch basis spurious long range components of the self-energy may (and typically do) arise, as evident from simple Fourier transform arguments. The diagonal approximation has been found 64 to have little effect on the quasi-particle energy spectrum (at least when LDA wavefunctions are a reasonable starting point). Our investigations confirm this picture at both the HF and COHSEX level. As discussed in the following, the case of the complex band structure is more critical. First we focus on a tight-binding model. In Fig. 7 we report the real and complex band structures for such a model, according to App. A. In order to simulate the effect of a diagonal self-energy, we refer to a picture where the Σ correction can be modelled as a stretching of the bands (which may be different for valence and conduction states) plus a scissor operator applied to the HOMO-LUMO gap. A simple scissor and a scissor+stretching corrections are applied to the model in Fig. 7(a,b) (solid black lines for the original model, thick light gray lines for the corrected ones). The new Hamiltonian including the corrections is now longer ranged than the original nearest neighbors Hamiltonian h 0 , because of the non-local projectors used to express the scissor and scis-sor+stretching corrections. The different spatial decay of the pristine and corrected Hamiltonians is shown in the lower panels of Fig. 7. We then fit the corrected Hamiltonian by using again the NN model (dashed, red lines). This fitting emulates the real band structure, but using a short ranged (nearest-neighbors) Hamiltonian. The effect on the CBS is evident and sizeable. The simply shifted and stretched electronic structures leads to little corrections to the CBS, while much bigger corrections are obtained when considering the fitted short-ranged Hamiltonians. Because CBS measures the decay of the evanescent states in real space, it is not surprising that a method which does not update the wavefunctions (as the diagonal self-energy corrections) is not able to capture all of the physics involved in the change of the electronic structure.
In order to further numerically support this interpretation and to investigate the effect of the self-consistency, we have evaluated the HF and COHSEX self-energies for PA HF . At first we have done so non-self-consistently (self-energies evaluated on the LDA wavefunctions), with and without the diagonal approximation. Then selfconsistent 91 COHSEX results are provided. HF results are plotted in Fig. 8(a,b), COHSEX data in panel (c). LDA β max is shown with a dashed line in the CBS panels as a reference. Regarding the real band structure, we found that the inclusion of off-diagonal matrix elements is not very relevant for PA, the bands being almost overlapping for both HF and COHSEX [diagonal data shown in panels (a,c) by dashed gray lines]. The situation is different for the CBS, as it is highlighted in Fig. 8(a). The HF 0 correction of β max from the full off-diagonal representation is almost twice as big as the diagonal correction. This confirms the behavior observed with the models in Fig. 7.
This observation also correlates with the decay of the HF 0 (LDA) Hamiltonian reported in the lower part of panel (a). As for the models [Fig. 7], the diagonal representation induces a much longer (and unphysical) de- . Real and complex band structure for PAHF at the HF and COHSEX levels. Panel (a) shows nonself-consistent HF results with a diagonal (gray lines) and fully off-diagonal (black lines) representation of the correction. Panel (b) reports the same data for the self-consistent HF solution, as computed from SaX (solid black lines) and CRYSTAL (green triangles). Panels (c) shows COHSEX data: Diagonal non-self-consistent (gray lines) and fully selfconsistent (black lines) data are reported. In both cases, a NN tight-binding fit of the real and complex band structures is performed (dashed red lines). Lower panels show a measure of the spatial decay of the COHSEX and HF Hamiltonian matrices on the WF basis.
cay. The same situation is found for COHSEX (the offdiagonal results at the first SCF iteration are not shown). The proper Hamiltonian decay (black line, circles) is clearly longer ranged than the LDA results, because of the non-local contribution of the exchange operator. The decay of the exchange potential is driven by that of the density matrix, which in turn is related 29 to β at the branch point. Being β max typically underestimated at the LDA level, so is the decay of the HF 0 Hamiltonian. The effect of self-consistency of HF (as well as COHSEX) is then to reduce such over-delocalization and to produce shorter ranged self-energies. This is shown in the decay plot of panels (b,c).
This behaviour has strong consequencies regarding the quality of the NN model fit. In the case of local and selfconsistent hybrid functional calculations from CRYS-TAL, the model fit (dashed red lines in Figs. 4, 5) works very well when compared with the full calculations, despite its semplicity. This is not the case when comparing with the diagonal corrections of Fig. 8(a,c). The offdiagonal representation improves the situation but does not solve the problem. The failure of the model fit in Fig. 8(a) is in fact mostly related to the decay of the exchange operator. A much better fit is then found when full self-consistency is included, as in Fig. 8(b,c). Finally, we have also compared the HF results from SaX (planewaves and pseudopotentials) and CRYSTAL (localized basis, full electron). Results are reported in Fig. 8(b) by the black solid lines and the green triangles, respectively. We find excellent agreement for both the real and complex band structures, confirming that the results from the two codes are well comparable and free of any systematic error.
In the light of the discussion above, when presenting the diagonal G 0 W 0 data for PA and PPV [Figs. 4(d),5(d)], the CBS is better described by the interpolated data (dashed lines) instead of that directly calculated from the diagonal GW corrections (solid thin lines). Similar conclusions about the importance of describing changes to wavefunctions when applying GW to transport calculations are reported also in Refs. [49,50,93].

B. Electronic structure
Now we can turn to discussing the accuracy of the electronic structure calculations for conjugated polymers. First of all we note that our results for PA 1 and PPV are in good agreement with previously published [83][84][85]92 G 0 W 0 results. In particular, for PA 1 we obtain a GW gap E g = 2.05 eV, to be compared with 2.1 eV (Ref. [85,83]) and 2.13 eV (Ref. [92]). These results are obtained for isolated chains of PA. In the case of a crystal, the gap shrinks 84 to 1.8 eV due to interchain interactions. For the isolated chain of PPV, Rohlfing et al. 83,85 found E g = 3.3 eV (G 0 W 0 ), which compares reasonably well with our G 0 W 0 result of 3.09 eV [see Tab. IV]. In general, the overall shape (including the band widths) of the GWcorrected band structures for PA and PPV computed in this work is in excellent agreement with that of Ref. [85].
From a methodological point of view, the accuracy of G 0 W 0 corrections (based on the LDA electronic structure) for organic molecules has been recently widely addressed, 59,62,[94][95][96] comparing with different implementations of self-consistent GW and experimental results. G 0 W 0 (LDA) is found to underestimate ionization potentials more than in extended systems, suggesting that a certain degree of self-consistency tends to improve on the results. Moreover, self-consistency is found to further lower the HOMO level and increase the fundamental gap. This leads to larger estimates of β max .
In terms of a direct quantitative comparison with ex-periments, some issues have to be taken into account. First, electronic structure (and optical) measurements for polymers typically distinguish between crystalline grains and amorphous regions. Isolated chains are considered to resemble more (and to be used as rough models for) the amorphous regions. Clearly, a direct theory-experiment comparison may suffer from systematic errors (interchain interactions, medium polarization, electrostatic effects). These features generally tend to reduce the fundamental gap wrt that of the ideally isolated chain. With this in mind, we can compare data calculated here with experimental data from photoemission (PES) or scanning tunneling (STS) spectroscopies. Rinaldi et al., measured 97 the electronic gap of PPV films (on a GaAs substrate) by means of STS. They were able to estimate E g ∼ 3 eV. Kemerik et al. 98 used also STS and found the fundamental gap of PPV [film deposited on Au (111)] to be around 2.8 eV. All of these results are to be reasonably considered as lower bounds of the theoretical gap for the isolated PPV chain.

C. Complex band structure: trends
As can be directly inferred from the model described in App. A, as well as from Ref. [8], the important parameters that determine the behavior of β(E) (and β max ) are the band gap E g , and the effective band widths of the states around the gap, given eg in terms of the hopping t 1 . While our model (see App. A) includes a second parameter t 2 to describe the difference in the band widths of the frontier bands, considering that the ratio t 2 /t 1 ranges from 0.1 to 0.01 or less, corrections to the model due to t 2 are not particularly relevant for the cases studied here. According to Eqs. (A4,A18), β max is mostly determined by the E g /t 1 ratio. Even though this is just a simplified NN model, our numerical investigations suggest that the model is widely applicable (using a folding technique as in the case of PPI when needed).
In general the band gaps are expected to increase with the fraction of non-local exchange included in the hybrid functional. For covalently bonded systems, the same trend is expected for the band widths. Numerically, in all our examples, while increasing the energy gap, HF increases also the band widths. The same trend is found for all the hybrid functionals we have investigated. Since both E g and t 1 increase, it is not trivial to understand a priori which mechanisms would dominate. Indeed, it is very clear that the band gap opens more than the band width, leading to a clear trend of β max increasing when a larger fraction of exchange is included in the calculation. The same trend is also found for the GW results, even though we had to extrapolate the CBS from the real band structure by using the model fitting (as described in details in Sec. V A).
In Fig. 9 we report a synthetic view of all the computed values of β max (times a/2, a being the polymer lattice parameter), including all electronic structure methods for PA (PA 1 , PA 2 are shown; PA HF is not shown because almost overimposed to the other PA geometries), PPV and PPI. Those data are plotted against the E g /t 1 value. The ideal curve from the tight-binding NN model is reported (dashed red line). The agreement between the computed and the modelled data is remarkable for all the cases studied. HF data are reported with empty symbols, showing in general a slightly worse agreement (as already discussed in Sec. IV C). On the basis of the above relations, and according to the results reported in Fig. 9, we suggest the use of the model fitting to extrapolate information about β max from experiments able to investigate the electronic structure. This would allow for an indirect measure of the CBS and the related parameters (as β max ).

D. Comparison with transport data
As discussed in the introduction, in order to compute the β decay of the current (or conductance) in a metalinsulator-metal (MIM) junction, we need two different ingredients: (i) the knowledge of the CBS for the infinitely long insulating region, and (ii) the position of the Fermi level of the MIM junction wrt the band structure of the insulator. This is depicted in Fig. 1. Assuming that at low bias the current is carried by the states close to the Fermi level of the junction, we basically need to know how these states decay into the insulating region (the polymer, in the present case). Once the CBS is known, either the Fermi level alignment is computed explicitly or it is estimated on the basis of physical considerations. While the direct calculation is feasible (but demanding) and in some cases necessary, it is also possible to give a first estimate of the Fermi level according to the so-called MIGS (metal-induced gap states) theory. 8,27,99 As proposed by Tersoff,27 if the metal DOS around the Fermi level is sufficiently featureless and the MIGS penetrate deep enough in the insulating region, the Fermi level of the metal-insulator junction is approximately pinned at the charge-neutrality level (often close to the midgap point) in order to avoid charge imbalance at the interface. The charge neutrality level can be easily identified from the CBS, as the energy where β(E) reaches its maximum inside the gap. From this perspective, β max is a first estimate of the experimental β decay. In general, the Fermi level will move from the charge neutrality level, and the β decay will change accordingly. The physical reasons leading the Fermi level to shift are mostly related to the charge transfer (and dipole formation) at the interface. This explains why different chemical linking groups on the same molecule may lead to different values of β. In general, the β max value computed from the CBS may be regarded as an upper bound for β and thus as a lower bound for the ability of the insulating layer to allow the current to tunnel through the junction.
The above discussion stresses the fact that the experimentally measured β does not in general depend 8 only on the electronic structure of the insulator (especially the CBS), but also on the details of the interface (which determines the position of the Fermi level). For instance, recent measurements on alkanes (oligomers of PE) determined a β decay length of 0.71-0.76Å −1 for for -NH 2 terminations, 100-102 while it has been found in the range 0.8-0.9Å −1 (and more disperse) for thiolated molecules. 5,102-104 Recent calculations confirm this picture also for conjugated polymers connected to gold leads through different chemical groups. 9 In that work, the Authors have studied a number of oligomers including oligo-phenylenes (whose infinite polymer is poly-paraphenylene, PPP). In the case of the molecule connected to gold leads via thiol groups, after investigating the interface-DOS projected on the molecule, they find that the Fermi level aligns close to the midgap point (reasonably the charge neutrality level considered here), off by few tenths of eV (∼0.2-0.3 eV, the molecular gap being about 2 eV) at the LDA level. Since the basic unit of PPP (a phenyl ring) is the same as part of the monomers forming PPV and PPI, we assume the charge neutrality condition to be almost fulfilled if we were considering Au-PPV-Au and Au-PPI-Au junctions with thiol anchoring. We can thus estimate β max to be a good estimate of β, keeping in mind that the a small deviation from midgap would slightly decrease β.
Recent experiments 7,26 reported β = 0.3Å −1 for oligomers of PPI connected to gold through thiols. This number is in very good agreement with the value β max = 0.29Å −1 we found for PPI using PBE0 (see Tab. V). According to our findings for PA and PPV (see Tabs. III,IV), GW should give results for β comparable to PBE0, but slightly larger. In order to compare with other existing (experimental and theoretical) results for oligo-phenylenes (PPP in the infinite limit) 23,103,105,106 we would need to address separately the issue of the Fermi level alignment (tending to decrease β wrt the ideal β max ) and that of the phenyl twist-angle (which goes in the direction of increasing β). This will be the subject of future work. Coming to the case of PA, we note that recent experimental results 107,108 report β-values of 0.22 A −1 for molecules similar to oligo-acetylenes. While this number is in very good agreement with our GW (and PBE0) results for PA 1 (and consistent with PA HF , see Tab. III), the large variability of β with the structural parameters of PA does not allow us to be conclusive on the assessment of the theory vs the experiment for this specific case. Nevertheless, our findings suggest that the use of PBE0 or GW results, together with a proper determination of the Fermi level alignment, will provide a reasonable approximation. We also note that gap underestimation and a priori assumption of the validity of the MIGS theory tend to partly cancel each other, and fortuitous agreement of experimetally measured β values with LDA calculations may occur.

VI. CONCLUSIONS
In this work we have computed from first principles the real and complex band-structure of prototype alkylic and conjugated polymer-chains using a number of theoretical schemes, ranging from local and semilocal to hybrid-DFT, and GW corrections. The accuracy of these different methods has been evaluated and compared with existing theoretical and experimental data, both in terms of the electronic structure and transport properties. From the CBS the β decay parameter, which governs nonresonant tunneling experiments through metal-insulatormetal junctions, can be computed.
In doing so we have stressed the formal analogy of hybrid-DFT and GW (especially in the COHSEX formulation), and the interpretation of the hybrid-DFT electronic structure as an approximation to the proper quasiparticle spectrum. We have also described in detail how to interpolate GW results by using a Wannier function scheme. In this case we have found that while the real band structure is always well interpolated, the CBS needs the self-energy real-space decay to be properly treated (off-diagonal representation and self-consistency of the wavefunctions).
We have numericaly investigated four polymers, namely poly-ethylene (PE), poly-acetylene (PA), polypara-phenylene-vinylene (PPV), and poly-phenyleneimine (PPI). Our results compare well with the existing theoretical and experimental literature. Among the hybrid functionals studied, PBE0 results compare best with the G 0 W 0 electronic structure. While the band gaps may still have a non-negligible deviation from GW, the agreement is remarkable on the CBS and β coefficient. The comparison with transport data (when available) is also very promising. This suggests PBE0 as an efficient and reliable alternative to GW for these class systems, at least for transport properties. More generally, a systematic application of hybrid functionals to improve the accuracy of DFT-based electronic structure results is appealing, while further developments along the lines of Ref. [70] are probably needed.

VII. ACKNOWLEDGEMENTS
We thank D. Varsano In the case of one dimensional systems like conjugated polymers, numerical results for β can be rationalized in terms of a simple tight binding model as presented in Ref. [8]. In their work, Tomfohr and Sankey presented a two-band model which provides an analytical expression for the complex band structure (CBS) within the fundamental energy gap. This model can also be used to fit the real band structure of realistic systems in order to evaluate the CBS analytically. The model is described in terms of two inequivalent sites ( a,b ) with nearestneighbors (NN) hopping t 1 . These three parameters can be recast into E g (the fundamental gap), t 1 and a further shift of the energy levels, which has no physical meaning. Moreover, it is shown that in this case β max (the maximum value of β(E) within the fundamental gap) depends only on E g /t 1 , according to: . (A4) All the details are given in Ref. [8]. We note however that this model is unable to reproduce any difference in the widths of the two bands (thus resulting in the CBS maximum being located at mid-gap), while in our realistic simulations we typically find the LUMO bandwidth to be somewhat greater than that of the HOMO. We have then generalized the above to the following three-parameters model, which is more suitable to determine the physical relevant quantities from our simulations. Extending the previous model, we add second NN interactions (with strength t 2 ) between equivalent sites, t 1 being the hopping between inequivalent sites as in the previous model. The model Hamiltonian is the following: where a, b indicate inequivalent sites and R is a cell index. The model is pictorially described in Fig. 10. The parameters t 1 , t 2 are in general complex numbers. Taking t 1 , t 2 to be real for semplicity, the analytical expressions of the energy bands read: E 1,2 (k) = Σ + 2t 2 x(k) ± ∆ 2 + 2t 2 1 x(k) where, assuming a > b , we have set: (A7) x(k) = 1 + cos(ka); (A8) In this picture E c,v are the onsets of the valence and conduction bands, a is the lattice parameter, k runs from − π a to π a . Under the condition |t 2 | |t 1 |, the band gap of the model is still at the Brillouin zone edge. A differenc choice of the relative phases of the model parameters would be needed to have the band gap at Γ (which is the case for PPI).
It is also possible to express the results of the model in terms of more physical parameters such as the energy gap and the HOMO and LUMO band widths; W c = ∆ 2 + 4t 2 1 1/2 − ∆ + 4t 2 (A11) Because in the realistic calculations the bands of interest may cross other bands far from the k corresponding to the gap, we consider partial ( W c , W v ) and not full band widths. These quantities are defined by the amplitude of the bands in a limited range of the BZ around the fundamental gap. This yields a much better agreement of the model bands with the calculated bands close to E g , which is the energy range of interest. In order to extract the parameters E g , t 1 , t 2 from our calculations we used the following relations (under the restriction that the band gap is direct and located at the Brillouin zone edge k = π/a): Following Ref. [8], once we have parametrized the model, we can give an analytical expression for the CBS, which in our case reads: where Eq. (A1) connecting β to γ holds unchanged. The maximum of γ(E) can be found analytically, and the resulting expression can be further simplified under the assumption |t 2 | |t 1 |; γ(E βmax ) ∼ 1 + 1 8 The CRYSTAL09 software package 76 perform calculations based on the expansion of the crystalline orbitals as a linear combination of a local basis set consisting of atom centered Gaussian orbitals. A 6-31G* contraction double valence (one s, two sp and one d shells) quality basis sets have been selected to describe carbon and nitrogen atoms; the most diffuse sp (d) exponents are α C = 0.1687 (0.8) and α N = 0.2120(0.8) Bohr. −2 The hydrogen atom basis set consists of a 31G* contraction (two s, one p shells): the most diffuse s and p exponents are 0.1613 and 1.1 Bohr. −2 The self consistent field procedure was converged to a tolerance in the total energy of ∆E = 2 · 10 −7 Ry per unit cell. 109 Reciprocal space sampling was performed on Monkhorst-Pack grid with 12 kpoints. The thresholds for the maximum and RMS forces (the maximum and the RMS atomic displacements) have been set to 0.00090 and 0.00060 Ry/Bohr (0.00180 and 0.00120 Bohr).
The calculations performed with Quantum-ESPRESSO (QE) adopt a plane waves basis set and norm-conserving pseudopotentials to describe the ion-electron interaction.
The kinetic energy cutoff has been set to 45 Ry for wavefunctions. For ionic relaxation, total energy LDA and GGA calculations use a Monkhorst-Pack grid of 8, 6, and 6 k-points for PE, PPV and PPI respectively (PA geometries are taken from the literature and not relaxed with QE). The convergence threshold on the atomic forces has been set to 10 −3 Ry/Bohr. A minimum distance of 20 Bohr between chain replica is used.
When performing GW calculations using SaX, the kpoints grids have been made finer by using 50, 50 and 20 k-points for PA, PE and PPV, respectively. The longrange divergence of exchange-like Coulomb integrals is treated using a generalized version 110 of the approach given by Massidda et al. 111 The same approach has also been used when performing hybrid-DFT calculations with QE. Note that other schemes to treat exchange in one-dimensional systems have been proposed. [112][113][114][115] The Godby-Needs plasmon-pole model 67 has been used, setting the fitting energies at 0.0 and 2.0 Ry along the immaginary axis, for all the cases. A kinetic energy cutoff of 6 Ry has been used to represent the polarizability P and the dynamic part of the screened Coulomb interaction W on a plane-wave basis, while a cutoff of 45 Ry has been used for the exchange operator. In order to converge the sums over empty states for the polarizability (self-energy), a total number of 288, 288, 288 (288, 288, 608) states has been used for PA, PE, and PPV respectively. This corresponds to an equivalent transitionenergy cutoffs of 53, 53, 35 eV (53, 53, 44 eV). Interchain distance has been increased to 30 Bohr to control spurious interactions of periodic replica. The QP corrections are computed by evaluating the diagonal matrix elements of the self-energy operator nk|Σ xc |nk , unless explicitly stated.