Signature of the Dirac cone in the properties of linear oligoacenes

Abstract

Organic electronics offers prospects of functionality for science, industry and medicine that are new as compared with silicon technology and available at a very low material cost. Among the plethora of organic molecules available for materials design, polymers and oligomers are very promising, for example, because of their mechanical flexibility. They consist of repeated basic units, such as benzene rings, and the number of these units N determines their excitation gap, a property that is often used in proposals of organic photovoltaics. Here, we show that contrary to a widely held belief the magnitudes of excitation gaps do not always decay monotonously with N, but can oscillate due to the presence of a 'Dirac cone' in the band structure. With an eye on the more fundamental question how a molecular wire becomes metallic with increasing length, our research suggests that the process can exhibit incommensurate oscillations.

ARTICLE

Received 19 Mar 2014 | Accepted 14 Aug 2014 | Published 23 Sep 2014

Organic electronics offers prospects of functionality for science, industry and medicine that are new as compared with silicon technology and available at a very low material cost. Among the plethora of organic molecules available for materials design, polymers and oligomers are very promising, for example, because of their mechanical exibility. They consist of repeated basic units, such as benzene rings, and the number of these units N determines their excitation gap, a property that is often used in proposals of organic photovoltaics. Here, we show that contrary to a widely held belief the magnitudes of excitation gaps do not always decay monotonously with N, but can oscillate due to the presence of a Dirac cone in the band structure. With an eye on the more fundamental question how a molecular wire becomes metallic with increasing length, our research suggests that the process can exhibit incommensurate oscillations.

DOI: 10.1038/ncomms6000

Signature of the Dirac cone in the properties of linear oligoacenes

Richard Korytr1, Dimitra Xenioti1,2, Peter Schmitteckert1,3, Mbarek Alouani2 & Ferdinand Evers1,3,4

1 Institut fr Nanotechnologie, Karlsruhe Institute of Technology (KIT), PO Box 3640, 76021 Karlsruhe, Germany. 2 Institut de Physique et Chimie des Matriaux de Strasbourg, UMR 7504 CNRS-UdS, 23, Rue du Loess, BP 43, 67034 Strasbourg, France. 3 DFG Center for Functional Nanostructures, Karlsruhe Institute of Technology (KIT), Wolfgang-Gaede-Str. 1a, 76131 Karlsruhe, Germany. 4 Institut fr Theorie der Kondensierten Materie, Karlsruhe Instituteof Technology (KIT), Wolfgang-Gaede-Str. 1, 76128 Karlsruhe, Germany. Correspondence and requests for materials should be addressed to R.K.(email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6000

Oligoacenes are molecules that consist of fused benzene rings (Fig. 1a). The shortest examplesbenzene, naphthalene and anthraceneare of fundamental impor

tance in chemistry. While the synthesis of gas-phase oligoacenes with more than six rings (hexacene) remains a challenge1, there are reports that attempts are on the way for on-substrate synthesis2. On the theory side, oligoacenes have been subject of an intensive debate due to pronounced correlation effects3,4. Much of the interest in this class of molecules is driven by organic electronics5 and photovoltaics6,7, where oligoacenes and their derivatives nd frequent applications812. From a broader perspective, a study of systems with linearly fused rings, nanographenes, offers means to learn how properties of a molecular wire converge to the one-dimensional limit.

As the number of fused molecular units increases, the optical gap either approaches a constant or vanishes. In the former case, the innite chain is a band insulator and in the latter it is a metal. Metallic behaviour is usually manifested by a single band that crosses the Fermi energy. This band is partially lled and the gaps of chains of nite length should tend to zero as quickly as 1/N, N being the number of units.

The 1/N behaviour merely reects the well-known level spacing of the particle-inside-a-box model. Hence, the common belief that the decay is generally monotonous (possibly, with decorating even-odd oscillations). However, in this work, we demonstrate that the gap can decay in a non-monotonous way, with strong and even incommensurate oscillations. Our calculations based on density-functional theory (DFT, detailed in Methods) show that the gap D Ng of the rst 10 oligoacenes drops to zero quickly, but it rises up again and drops repeatedly with periodicity of 11 units (rings). By a sophisticated post-DFT analysis, we clarify the conditions for observing the oscillations in oligoacenes. Remarkably, this behaviour has not been reported before in spite of the

intensive research focus on these molecules. As we argue below, the gap oscillations originate from a band crossing at the Fermi level, that is reminiscent of the Dirac points and Dirac cones in certain lms, such as graphene, and bulk materials. A consequence of the oscillations is that molecules with similar excitation gaps can have very different lengths. In practical terms, our research implies that the molecular length can be a parameter that can be optimized in addition to the energy gap.

ResultsSingle-particle spectrum. The segmented (piecewise-smooth) form of the oligoacene gap, shown in Fig. 1b, indicates that an orbital crossing occurs at the extremal points of the gap function. We conrm this by inspecting the KohnSham eigenvalue spectra in Fig. 1c. The orbital energies evolve as almost straight lines as a function of the inverse chain length and create a web-like pattern. For a given N, the highest occupied molecular orbital (HOMO) is given by the blue symbol of the highest energy; the lowest unoccupied molecular orbital (LUMO) corresponds to the red symbol of the lowest energy. Close to certain crossings, the LUMO and HOMO interchangein this situation, the gap has a minimum. The local maxima of the gap manifest another level crossing.

Band crossing. To understand the origin of the oscillations and level crossings, we invoke the band structure of the innite chain, polyacene, in Fig. 2. As one passes from the G point towards the corner of the Brillouin zone, we encounter a crossing of the valence and conduction bands at kD 0.9102p/a. This peculiar band

structure of polyacene wires was discussed rst by Kivelson and Chapman13 based on an analysis of a tight-binding Hamiltonian with rst (t), second (t0) and third nearest neighbour hopping (t00).

As they argued, t00 shifts the conduction and valence bands against each other. If t and t00 share the same sign, the shift is negative and we witness a non-avoided band crossing. From our DFT-calculation, Fig. 2, we estimate t00 0.1t 0.31 eV by matching

to the tight-binding model.

As a consequence of crossing, electrons and holes with energies sufciently close to the Fermi level have a piecewise linear dispersion (see Fig. 2) reminiscent of the Dirac cones in graphene. Just as in graphene, the Dirac point in polyacene is a property of the conjugated p electron system. However, in graphene, the position of the Dirac points in the Brillouin zone is dictated by the symmetry of the honeycomb lattice, while in polyacene it is controlled by t00.

1

3.2

0.9

3.4

0.8

Occ.Empty

0.7

3.6

0.6

g(eV)

Eigen energy (eV)

3.8

0.5

4

4

0.4

2

0.3

4.2

0.2

b E F(eV)

4.4

0

0.1

kD

4.6

2

0 10 0.05 0.1 0.15

Figure 1 | Electronic structure of oligoacenes. (a) Ball-and-stick representation of benzene, anthracene and hexacene with N 1, 3, 6 rings

(blue: hydrogen, orange: carbon). (b) optical gap of oligoacenes as a function of the number of rings N calculated in DFT. (c) Flow of orbital energies as a function of 1/N, showing that a repeating interchange of orbitals occurs as the number of rings grows. Blue symbols are occupied orbitals and red are unoccupied. Close to N 10, N 21 and N 32, there

is an interchange of frontier orbitals. The data points are connected by straight lines for visual guidance.

4

5

15

20

25

30

35

0

N

1/N

0 /a

k

Figure 2 | Band structure of polyacene highlighting the band crossing at a Dirac point. Blue lines are occupied bands, red lines are unoccupied. The gap between the valence and conductance bands decreases as we move from the G point (k 0) to the zone border, where it becomes negative.

a 2.462 is the unit cell length. The band structure was calculated in

DFT, see Methods for details. Inset: ball-and-stick model of polyacene with carbon (yellow) and hydrogen (blue) atoms.

2 NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6000 ARTICLE

Origin of oscillations. At least in the case of sufciently long oligoacenes, one should be able to derive their spectrum from the polyacene band structure by imposing selection rules for wave numbers due to hard-wall boundary conditions (zone-folding). For a given number of rings N, only selected wave numbers are allowed. We choose the following set

kNm

m N 1

N

for

N44. We attribute this to the fact that the relaxation effects compress the length of the rings at molecules edges. Since the total length shrinks, the level spacing and gap increase.

We have seen that on DFT level, band structure effects control the excitation gaps of oligoacenes. It is known that the Kohn Sham excitation gaps obtained from the DFT ground-state calculation agree with optical excitations only up to corrections that are incorporated in time-dependent DFT. These corrections include a Hartree and an exchange-correlation kernel that can renormalize the excitation energies. In ref. 14, these dynamical correction have been calculated for the smaller system sizes N 26 (with a hybrid functional) and turn out to be below 10%

of the bare KohnSham gap. As screening effects generally enhance with decreasing gap, that is, increasing system size, it is reasonable to expect that dynamical corrections do not affect the picture that we present.

Charged excitations. We now argue that the charged excitation energies, that is, the ionization energies and electron afnities, also inherit the oscillatory structure from the Dirac cone. In this section, we provide an analysis of results by DFT and various post-DFT techniques. The latter can exhibit signicant deviations from the DFT estimates and therefore correct for artifacts of approximate exchange-correlation functionals1416.

As it is well known, in KSDFT, the ionization potential I agrees with the HOMO energy of the charge-neutral molecule with Nel electrons: I EHOMO(Nel). The correspond

ing electron afnity A can be obtained from a reference calculation that has been performed with one additional electron: A EHOMO(Nel 1).

In contrast to DFT, in Greens function methods (such as G0W0), the energies of charged excitations are already incorporated in the quasiparticle (QP) energies. The G0W0 QP band

structure is shown in Supplementary Fig. 1. It shows a band pattern very similar to the DFT results. Since computational demands do not allow a ne resolution near the Fermi points, we inspect the symmetry of the four Bloch states in Fig. 4, which provides evidence of a band crossing in G0W0. Therefore, one expects that the charged excitations support oscillations analogous to the ones already observed in the optical gaps. The

pa ; m 1; . . . ; N; 1

which applies for a guitar string of the length a(N 1), a is the

length of one ring. In Supplementary Note 1, we inspect the boundary condition in detail and demonstrate in Supplementary Fig. 2 that the optimal quantization length is the one given in equation 1. Let us denote Eg(k) : Ec(k) Ev(k), where Ec,v are the

conduction and valence band energies. The optical gap D Ng will be given approximately by the energy Eg k N

, where k N

denotes the wave number from the allowed set, that is closest to the Dirac point, kD. Hence, k N

labels the HOMO/LUMO pair.

As we increase the number of rings, k N

will move and eventually cross the Dirac point, where the HOMO and LUMO interchange. This explains that the periodicity of the gap oscillations of the oligoacenes is controlled by the position of the Dirac cone in the Brillouin zone. More precisely, the period is inversely proportional to p/a kD, the distance to the zone

boundary. Similarly, the maxima of the gap are also xed by band structure consideration: it is only when the Dirac point lies in the middle between two consecutive wave numbers from the allowed set. If a local maximum occurs at N0 with value D N0g, its upper bound must be uFp/a/(N0 1), where uF is the Fermi velocity

(numerically uFp/a 6.485 eV).

The above considerations suggest natural energy and length scales. In Fig. 3, we show the gap D Ng expressed in these natural units (blue trace). Remarkably, our plotting exhibits a sawtooth form with well developed periodicity (with the exception of very small molecules). Furthermore, in the same plot we show Eg k N

the gaps reconstructed from folding the polyacene band structure (red traces). Remarkably, both traces lie almost on top of each other. Naively, corrections of the order of 1/N to the zone-folding procedure should be expected, which means a correction of the order of 100% for benzene, N 1. In practice,

the nite-size correction turns out to be much smaller. For benzene, the deviation is only 11% (tetracene 1%, hexacene 8%).

We comment that the molecular gaps D Ng in Fig. 3 already involve relaxation effects that can not be accounted for in the zone folding. We see that D Ng is always higher than Eg k

1.8

1.6

Molecule (DFT)

Zone folding

a(N +1) / ([afii9843]h F) g

1.4

1.2

1

+

0.8

EF

0.6

+

[afii9843]

0.4

kD a

+

0.2

0 0.5 1.5 2 2.5

1

3

(N + 1)(1 akD/[afii9843])

Figure 3 | Comparison of optical gaps of oligoacenes from DFT and zone folding. The gap energies are divided by the factor puF/(N 1)/a to scale

away the leading behaviour. The horizontal axis is rescaled by the reciprocal space distance of the Dirac point to the zone border. Blue squares are calculated for a set of geometry-optimized molecules from DFT. Red circles are calculated from the Bloch bands of the innite chain by Brillouin zone folding. The data points are connected by lines for visual guidance only.

+

Figure 4 | Evidence of band crossing in G0W0. (a) Charge density of QP states of polyacene in the G0W0 approximation. We plot charge density in the plane offset by 0.7 from the molecular plane. The colour changes linearly from blue, corresponding to zero. Shown are four states whose energies lie in the vicinity of the Fermi level, as sketched in b by four circles. The QP band whose energy grows with k has a nodal plane along the CC bonds of the rungs and we indicate it by the sign. The other two states

that pertain to the band of decreasing energy are labelled as and lack the

nodal plane.

NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6000

renormalized parameters of the Dirac cone are compared in Supplementary Table 1.

In what follows, we will represent the number of rings N by a dimensionless number x (N 1)(1 akD/p), where kD is the

Dirac wave number and a is the lattice parameter. Brillouin zone folding predicts gap minima at integer x. The orbital energies Ei (i

stands for HOMO or LUMO) will be represented by a dimensionless quantity Ei : a N 1puF Ei EF

.

First, we discuss the behaviour of the HOMO and LUMO (KohnSham) energies on DFTPerdewBurkeErnzerhof (PBE) level. The LUMO energy enters because we replace EHOMO(Nel 1) by ELUMO(Nel) thus ignoring effects of the

derivative discontinuity. In Fig. 5, we show HOMO (symbols red, lled) and LUMO energies (symbols red, open), along with their counterparts from the Brillouin zone folding (solid lines). As expected, with increasing number of rings the frontier orbitals cross near x 1 in a linear fashion.

The frontier orbitals can be classied according to the symmetry with respect to a mirror plane that is normal to the molecular plane cutting the molecule along the long axis into two identical pieces. Notably, the energy of symmetric (antisymmetric) orbitals decreases (increases) with x. Hence, at x 1, the

band theory predicts a change in the symmetry of the HOMO orbital. The reverse symmetry transitions appear at half-integer steps in x. We have veried that this is indeed the case for the molecules calculated in PBE.

In the following, we explore an interesting consequence of the symmetry transition for the energies of charged excitations. To this end, we present HOMO and LUMO energies calculated in PBE0. PBE0 is a (hybrid) density functional that replaces a certain fraction l of the PBE exchange by an exact (HartreeFock type)

exchange. In conventional PBE0, one takes l 25%. The derivative

discontinuity, that is absent in PBE, is thus partially restored.In Fig. 5, the (dimensionless) PBE0 energies EHOMO,LUMO are

shown (green symbols). They deviate signicantly from the PBE zone-folding estimates. This difference is due to l: it vanishes continuously in the limit l-0. The PBE0 data brings two

qualitatively new features. With hybrid functionals, LUMO and HOMO energies contain a fraction of the charging energy. In PBE0, the latter is entirely given by HartreeFock exchange, which does not account for interaction screening. Hence, for long linear chains, the charging energy will be dominated by a spurious logarithmic length dependence. The slowly increasing deviation between the HOMO/LUMO energies and band lines for the rst 10 molecules can be understood in this way. At N 11,

we witness an abrupt change: a slope inversion and a step in the PBE0 trace. It can be understood in the following way. The electronic structure data for PBE0 and PBE are adiabatically connected via the l parameter. Hence, for every system size N, the frontier energies obtained with PBE0 can be thought of as originating from a PBE reference result corrected by perturbative terms that are controlled by l. The correction terms describe the charging energy and therefore contain Coulomb-matrix elements of the PBE reference calculation. At N 11, the symmetry of the

frontier orbitals of the PBE reference calculation changes. The Coulomb-matrix elements react to that change and consequently the PBE0 traces exhibit a jump. In the presence of screening beyond PBE0, we expect that the spurious logarithm disappears in the limit x-N and the step will be weaker.

The G0W0 QP energies are also shown in Fig. 5. It is easy to understand the differences between all methods presented if we recall that the PBE0 values contain a fraction 14 of the charging energy, while the G0W0 values fully include charging effects. Consequently, for benzene, the difference in orbital energies E GWHOMO E PBEHOMO is approximately four times the

difference E PBE0HOMO E PBEHOMO and the same holds true for LUMO.

We notice that for longer molecules, the relative difference drops below four, because of the screening mechanism in G0W0 that is absent in PBE0. Due to high computational demands, we could not obtain QP energies for systems beyond N44 rings.

Correlation effects. To rigorously establish the validity or our conclusions even in the presence of strong correlations (where ab initio methods tend to fail), we analyse a (parameterized) model system, the oligoacene Hubbard model. We solve it numerically at different values of the interaction strength by employing the density-matrix renormalization group (DMRG), see Methods for details. We rst discuss a model with spin-less fermions. Results for the optical gap are show in Fig. 6 for various values of N and nearest neighbour interaction U1. As we see, the N oscillations of the optical gap survive even up to interactions comparable in strength with the bandwidth. The period can be renormalized to the extent that for large interaction the oscillation period increases. At a certain critical value comparable to the bandwidth, U1E5 eV, the period diverges. At large interactions U1445 eV,

2

PBE PBE0

GoW0

1.5

1

HOMO,LUMO([afii9841])

0.5

0

0.5

1

12

1.5

10

0.5 1

[afii9841] =(N+1)(1akD/)

1.5

(N+1) g(eV)

2

8

6

4

2

4

Figure 5 | Comparison of charged excitation energies obtained with different methods. We plot dimensionless HOMO and LUMO energies

Ei : a N1p hu Ei EF

, over the length of the oligoacene. The parameters a, uF, EF and kD were obtained from a separate PBE band structure calculation.

PBE results are represented by red symbols; red lines are estimates from Brillouin zone folding. Green symbols represent the PBE0 counterparts and magenta are the G0W0 QP energies. Empty (lled) symbols stand for

LUMO (HOMO).

10

8

6

U 1(eV)

2

0 5 10 15 20 N

25 30 35

0

Figure 6 | DMRG results of the optical gap in the spin-less case with varying interaction strength. With growing interactions the oscillation period increases until it diverges at an interaction strength approximately equal to the bandwidth.

4 NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6000 ARTICLE

we nd a (twofold degenerate) nearest neighbour charge density wave ground state.

Our DMRG results for the spin-full calculation at low interaction strength are depicted in Fig. 7a. The data conrms that the oscillation phenomenon is stable against the presence of weak interactions. Figure 7b displays the dependency of the optical gap on the on-site interaction, U0. Three molecules are shown with a length N that is close to the optimal length NminE10 where the gap is expected to exhibit particularly small excitation gap. In the perturbative regime U0oo2t, switching on

U0 even stabilizes the oscillations and makes them more pronounced in the sense that the ratio between the minimal gap near NminE10 and the neighbouring gaps goes through a maximum at nite U0. At interaction strength U0t2t, the excitation energy for N 10 is signicantly smaller than for the

neighbouring molecules N 9,11. Only at U0E6 eV, the two

traces (N 9,10) cross. At larger values, a monotonous order (in

N) of the traces is established.

By performing an independent DMRG calculation in the spin-sector S 1, we have checked that the optical gap coincides with

the spin gap for the values of U considered in this manuscript.

The scaling of the excitation gap with system size is shown in Fig. 8. For smaller molecules, an approximately exponential decay is observed that can be tted by

DEgN D0e N=N ; Nt2N

with D0E7.8 eV and N*E2.3. Since the quasiexponential decay is already seen at weak interactions including U0 0, it is unrelated

to correlation effects, such as ground state magnetism. Interaction effects become relevant only at N\N* and very strong values of

U0 that signicantly exceed the bandwidth E2t. In this regime, the gap oscillations disappear and a monotonous behaviour sets in, see Fig. 8

DiscussionOver the years, there has been intensive theoretical and experimental work that addresses the nature of the ground-state of oligoacenes, for review, see refs 1,4. The topic of interest is whether the ground state could be magnetic (open shell, radical) or at least close to a magnetic (or other) instability. We discuss a selection of earlier work that appears most relevant to us in the light of our new ndings.

In DFT studies employing hybrid functionals, broken-symmetry, spin-polarized ground states have been reported17. Since approximate functionals underestimate quantum uctuations, they tend to overestimate the tendency for symmetry breaking. Therefore, such studies are difcult to interpret. Indeed, our DMRG calculations indicate that the ground state does not exhibit a (simple) broken-symmetry phase. To be consistent with this nding, we employ in our DFT studies non-magnetic, closed-shell calculations.

To address correlation effects, other methods have been applied. Where only relatively small system sizes could be addressed, an extrapolation of excitation energies to the polymer-limit of large N was performed (theory3 Nr12, experiment18

Nr6). For instance, Hachmann, et al.3 performed a DMRG study on a model of polyacene similar to ours. These authors t their data for the spin-triplet gap in the manner DEstEDst dst

exp( N/Nst) with a non-vanishing offset DstE0.1440.017 eV

(orig: 3.330.39 kcal mol 1). The offset indicates the presence of strong correlations that have been interpreted as a singlet ground state with an open shell (di- or polyradical). A nite offset has also been reported in another DMRG study on an oligoacene-like Heisenberg model (up to N 28)19. Since the Heisenberg

Hamiltonian ignores charge uctuations, the latter result carries over to oligoacenes only in the Mott-limit of very large U0. In this

regime, our results are consistent with refs 3,19 in the sense that we conrm the absence of oscillatory behaviour, as shown in Fig. 8 for very large interaction U\2t. In this regime, our results are also in agreement with the study of ref. 20 where in the large U limit an enhanced mottness for an honeycomb like two-leg ladder in comparison to a quadratic two-leg ladder was shown.

In conclusion, we show by a combination of DFT, post-DFT methods (GW theory) and numerically exact Hubbard-type model calculations that oligoacenes exhibit oscillations in their fundamental excitation gaps with increasing length N. It has been shown that these oscillations result from the well-known Dirac-type dispersion of polyacene in combination with a zone-folding argument. This basic origin suggests that the effect is robust, which is an important requirement for prospective technological

12

1

0.8

0.6

U 0(eV)

10

6

(N+1) g (eV)

8

0.4

0.2

0

4

2

5 10 15 20

0

N

(N+1) g(eV)

43.5 3

2.5 2

1.5 1

0.50 0 2 4 6 8 10 12 14 16

U=0.0 U=1.0

U=15.0

Exp. decay

1

N = 9 N = 10 N = 11

g(eV)

0.1

0.01

U0 (eV)

Figure 7 | DMRG results of the optical gap for a spin-full Hubbard model. (a) Oscillations in the optical gaps are stable for small interactions U0 (vertical scale). The interaction has little or no inuence on the period of the oscillations. (b) gaps for increasing on-site interaction U0 for system sizes N 9,10 and 11, which are close to the rst minimum of the optical

gap. (Lines are guides to the eye.).

10

5 10 15 20

N

Figure 8 | Scaling the optical gap with the number of rings N in the spin-full model. Traces with three values of the on-site Hubbard interactionU0 0,1.0,15.0 eV are shown.

NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6000

applications. We have seen that the oscillatory modulation of the gap dependence may be hindered if strong interactions are active. At present, a realistic estimate of the interaction strength is lacking. However, even if the Coulomb repulsion wipes out gap oscillations in the gas phase, our study suggests that these band structure effects will revive when molecules are embedded in screening environments, such as electrochemical solutions21 or synthesized on surfaces, where the effective screened interaction will be several times reduced as compared with the gas-phase value.

Methods

Density-functional calculations. We employ the usual KohnSham density-functional method22,23 implemented in the Fhi-Aims package24. The optical gaps were calculated using the PerdewBurkeErnzernhof25 functional. The spin-restricted (closed shell) KohnSham equations were solved in a non-relativistic form, with relativity included a posteriori by rescaling the Hamiltonian eigenlevels, according to the zeroth order regular approximation26. The basis set for hydrogen (carbon) contained totally 15 (39) functions of angular momenta up to l 4 and

distinct radial dependences (the basis set is called tier2). For periodic calculation of the band structure, we sampled the Brillouin zone using 1,000 points in the reciprocal chain direction. For both periodic and non-periodic calculations, the ground state was reached with strict convergence criteria for the difference in density (10 6), total energy (10 7 eV) and forces (10 5 eV 1). Geometries were considered optimized only after the maximum residual force dropped bellow 10 3 eV 1.

The optical gap, shown in Figs 1 and 2, was approximated as the energy difference between the LUMO and the HOMO. A discussion of effects beyond this approximation is included in Results.

The G0W0 calculations. For the QP band structure and QP spectra of nite chains, we have used Vienna Ab initio Simulation Package (VASP)27,28. In VASP, the one-electron orbitals are expressed using the projector augmented wave basis set29. In this calculation, we used plane waves with 400 eV energy cutoff for the DFT starting calculation, and a cutoff of 300 eV for the response functions. First, standard KohnSham DFT calculations were performed, using the PBE25 exchange-correlation functional. Then, QP corrections were included on top of the KohnSham states at the non-self-consistent, single-shot, G0W0 level. After convergence checks, a frequency grid of 80 points and 600 bands in total was used. The QP energies were represented on a mesh of 11 k-points in the Brillouin zone.

The oligoacene Hubbard model and DMRG. Following Kivelson and Chapman13, we construct a tight-binding model with a single orbital site per carbon atom at half lling and including up to third nearest neighbour hoppings t,t0 and t00. We set t 3.1 eV, t0 0.341 eV, t00 0.31 eV (as estimated from our ab initio

band structure). We allow for on-site (U0) and nearest neighbour (U1) interactions, that we take as a variable parameters.

We perform an analysis based on the DMRG. DMRG is a well-established method suitable for studying many-body eigenvalue spectra of quasi-one-dimensional systems with strong correlations30,31. We study the oligoacene Hubbard model in a spin-less and spin-full fashion. In the spin-less case, we take U1 variable (U0 is meaningless). Discarded entropies during the DMRG sweeps are always below 10 7, typically much lower. In the spin-full case, results are parameterized by U0, while U1 0. We keep up to 8,000 states per block; the

dimension of the target space grows up to 4.3 107, the discarded entropies are

considerably below 10 4 (even for N 21) in each DMRG step, typically

signicantly smaller. We employ at least seven nite lattice sweeps with an optimized innite lattice warm up.

The optical gap is the difference between rst excited and ground-state energies.

References

1. Sun, Z. & Zeng, Z. Benzenoid polycyclic hydrocarbons with an open-shell biradical ground state. Chem. Asian J. 8, 28732885 (2013).

2. Treier, M. et al. Surface-assisted cyclodehydrogenation provides a synthetic route towards easily processable and chemically tailored nanographenes. Nat. Chem. 3, 6167 (2011).

3. Hachmann, J., Dorando, J. J., Avils, M. & Chan, G. K.-L. The radical character of the acenes: a density matrix renormalization group study. J. Chem. Phys. 127, 134309 (2007).

4. Bettinger, H. F. et al. Electronic structure of higher acenes and polyacene: The perspective developed by theoretical analyses. Pure Appl. Chem. 82, 905915 (2010).

5. Kuribara, K. et al. Organic transistors with high thermal stability for medical applications. Nat. Commun. 3, 723 (2012).

6. Lipomi, D. J., Tee, B. C.-K., Vosgueritchian, M. & Bao, Z. Stretchable organic solar cells. Adv. Mater. 23, 17711775 (2011).

7. Dabo, I. et al. Donor and acceptor levels of organic photovoltaic compounds from _rst principles. Phys. Chem. Chem. Phys. 15, 685695 (2013).

8. Bundgaard, E. & Krebs, F. C. Low band gap polymers for organic photovoltaics. Sol. Energ. Mat. Sol. C 91, 954985 (2007).

9. Anthony, J. The larger acenes: versatile organic semiconductors. Angew. Chem. Int. Ed. 47, 452483 (2008).

10. Lin, C., Wang, Y.-C., Hsu, S.-J., Lo, C.-F. & Diau, E. W.-G. Preparation and spectral, electrochemical, and photovoltaic properties of acene-modied zinc porphyrins. J. Phys. Chem. C 114, 687693 (2010).

11. Lee, J. et al. Singlet exciton ssion in a hexacene derivative. Adv. Mater. 25, 14451448 (2013).

12. Wilson, M. W. B., Rao, A., Ehrler, B. & Friend, R. H. Singlet exciton ssion in polycrystalline pentacene: from photophysics toward devices. Acc. Chem. Res. 46, 13301338 (2013).

13. Kivelson, S. & Chapman, O. L. Polyacene and a new class of quasi-one-dimensional conductors. Phys. Rev. B 28, 72367243 (1983).

14. Malloci, G., Mulas, G., Cappellini, G. & Joblin, C. Time-dependent density functional study of the electronic spectra of oligoacenes in the charge states

1, 0, 1 and 2. Chem. Phys. 340, 43 (2007).

15. Blase, X., Attaccalite, C. & Olevano, V. First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications. Phys. Rev. B 83, 115103 (2011).

16. Setten, M. J. v., Weigend, F. & Evers, F. The GW-method for quantum chemistry applications: theory and implementation. J. Chem. Theory Comput. 9, 232 (2013).

17. Jiang, D.-e. & Dai, S. Electronic ground state of higher acenes. J. Phys. Chem. A 112, 332335 (2008).

18. Angliker, H., Rommel, E. & Wirz, J. Electronic spectra of hexacene in solution (ground state. Triplet state. Dication and dianion). Chem. Phys. Lett. 87, 208212 (1982).

19. Qu, Z., Zhang, D., Liu, C. & Jiang, Y. Open-shell ground state of polyacenes: a valence bond study. J. Phys. Chem. A 113, 79097914 (2009).

20. Karakonstantakis, G., Liu, L., Thomale, R. & Kivelson, S. Correlationsand renormalization of the electronphonon coupling in the honeycomb Hubbard ladder and superconductivity in polyacene. Phys. Rev. B 88, 224512 (2013).

21. Capozzi, B. et al. Tunable charge transport in single-molecule junctions via electrolytic gating. Nano Lett. 14, 14001404 (2014).

22. Hohenberg, P. & Kohn, W. Inhomogeneous electron gas. Phys. Rev. 136, B864B871 (1964).

23. Kohn, W. & Sham, L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133A1138 (1965).

24. Blum, V. et al. Ab initio molecular simulations with numeric atom-centered orbitals. Comput. Phys. Commun. 180, 21752196 (2009).

25. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865 (1996).

26. van Lenthe, E., Baerends, E. J. & Snijders, J. G. Relativistic total energy using regular approximations. J. Chem. Phys. 101, 97839792 (1994).27. Kresse, G. & Furthmller, J. Efcient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169 (1996).

28. Shishkin, M. & Kresse, G. Implementation and performance of the frequency-dependent GW method within the PAW framework. Phys. Rev. B 74, 035101 (2006).

29. Blchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 1795317979 (1994).

30. White, S. R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. 69, 28632866 (1992).

31. White, S. R. Density-matrix algorithms for quantum renormalization groups. Phys. Rev. B 48, 1034510356 (1993).

Acknowledgements

We acknowledge the valuable discussions with S. Tretiak, O. Tal and T. Yelin. D.X. and M.A. acknowledge computational resources from HPC mso centre of Strasbourg and GENCI-CINES grant 2014-gem1100. Part of this work was performed on the computational resource bwUniCluster funded by the Ministry of Science, Research and Arts and the Universities of the State of Baden-WA14rttemberg, Germany, within the framework

programme bwHPC. R.K. and F.E. gratefully acknowledge the Steinbuch Centre for Computing (SCC) for providing computing time on the computer HC3 at Karlsruhe Institute of Technology (KIT) and the computing time granted by the High Performance Computing Center Stuttgart (HLRS) on the HERMIT supercomputer.

Author contributions

R.K. performed DFT calculations; P.S. performed the DMRG study; D.X. and M.A. did the GW calculations; R.K., D.X. and M.A. performed the PBE0 calculations; R.K. and F.E. analysed the data and wrote the manuscript.

6 NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6000 ARTICLE

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Korytr, R. et al. Signature of the Dirac cone in the properties of linear oligoacenes. Nat. Commun. 5:5000 doi: 10.1038/ncomms6000 (2014).

NATURE COMMUNICATIONS | 5:5000 | DOI: 10.1038/ncomms6000 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

& 2014 Macmillan Publishers Limited. All rights reserved.