The role of polymorphism in organic thin films: oligoacenes investigated from first principles

We investigate the cohesive properties of oligoacenes within the framework of density functional theory including van der Waals interactions. In comparison, we evaluate the local density approximation as well as gradient corrections, but also a widely used semiempirical procedure accounting for the long-range dispersive forces, in terms of their performance for the energetics of such weakly bound systems. Besides the cohesive energies we discuss in detail the surface energies which, in turn, allow for obtaining the crystal shapes based on Wulff's construction for the oligomer series from naphthalene to pentacene. In particular, we focus on comparing two different pentacene polymorphs, i.e. the so-called bulk structure and the thin film phase, the latter being predominately found in thin film growth. We not only study the impact of molecular conformation on the details of these polyhedra, but also the influence of the number of considered index planes and the role of the underlying exchange correlation functional. Based on the relaxed crystal structures for the two polymorphic phases, we compute the electronic band structures as well as the optical spectra. To account for excitonic effects in the latter, we solve the Bethe–Salpeter equation for the electron–hole pairs, thereby considering the coupling between resonant and anti-resonant terms.


Introduction
Exploiting the toolbox of organic chemistry, a huge variety of organic molecular crystals can be grown to exhibit certain physical properties. The molecular character of the building blocks enables functionalization of the material and hence the tailoring of its properties. One particular class of such structures are organic semiconductors based on π-conjugated molecules. Typical representatives like the oligophenylenes, oligoacenes or oligothiophenes exhibit band gaps in the visible range of the light spectrum making them suitable as active layers in opto-electronic devices. Owing to the fact that the electronic band gap scales with the inverse molecular length, some of the most promising molecules consist of 4-6 aromatic rings, i.e. are roughly 1.5-2.5 nm, long. This molecular size, in turn, makes the structural, electronic and optical properties highly anisotropic.
Another characteristic feature of such materials is their bonding. While strong covalent bonds are responsible for intramolecular interactions, the molecules are weakly bound within their crystalline phases, and hence known as typical van der Waals (vdW) systems. This fact has a tremendous impact on thin film growth. On the one hand, shallow energy potentials lead to polymorphism, since different molecular conformations exhibit only small energy differences. Moreover, the pronounced anisotropy and complex nature of the molecules as the building blocks make organic film growth fundamentally different from the inorganic counterpart. As a matter of fact, molecule-molecule as well as molecule-substrate interactions strongly depend on the relative orientations of the interacting objects. The deposition of a molecular layer on top of a substrate is determined by a delicate balance between the weak interactions among the 3 molecules as well as between molecule and substrate. It is this interplay which governs thin film morphology. For example, sexiphenyl can orient differently on a single substrate, depending on the growth conditions [1]. A more recent report by Hlawacek et al [2] has revealed a multilayer growth of sexiphenyl on mica in the form of terraced mounds, where a combined experimental and theoretical analysis has revealed a sizable Ehrlich-Schwoebel barrier. Depending on the nature and temperature of the substrate, the orientation of deposited pentacene molecules can vary from islands of standing molecules on inert and flat substrates, to flat lying molecules on metal surfaces [3].
Thinking about possible applications, the above-described situation does not easily allow for well-reproducible materials with distinct properties. In particular, the optical spectra strongly depend on the orientation of the material with respect to the polarization direction of light [4]. Due to the anisotropy mentioned above, the absorption onset sensitively depends on the orientation of the molecules. In the case of oligoacenes, even the nature of the lowest optical transitions is altered from a strongly bound Frenkel-type exciton to free electron-hole (e-h) pairs when going from polarization perpendicular to parallel to the long molecular axis [5]. This fact has an important implication, e.g. for the application in solar cells, where both, an absorption range in the frequency spectrum of sunlight is equally important as the desire for highly mobile free charge carriers.
For all these reasons, a profound understanding of organic thin film growth is the ultimate goal of many current investigations. From a theoretical point of view, a reliable description of the structural properties of organic semiconductors is a prerequisite towards modeling organic film growth. By tackling the energetics of molecular crystals as well as their interaction with metal substrates conclusions can be drawn in terms of thin film growth close to equilibrium where the surface energies of the respective constituents as well as the interface energy between the film and the substrate play important roles [6,7]. For both scenarios, it could be shown recently that the inclusion of vdW interaction is not only crucially important [8]- [11] but makes density functional theory (DFT) a powerful tool in this context.
In this paper, we focus on the effects of polymorphism on the structural, electronic and optical properties of pentacene, which is a prominent material due to its high-charge carrier mobility [12]. To this extent, we study two polymorphic phases, i.e. the Campbell bulk structure and the thin film phase in terms of their structural conformation and the resulting physical properties. We discuss the respective electronic band structures and finally show the impact on the optical spectra. The detailed investigation on the energetics in terms of cohesive and surface energies as well as crystal shapes also contains results for several members of the oligoacene series, which are naphthalene, anthracene and tetracene.

Theoretical background
Computational materials science largely relies on DFT for a quantitative description of physical and chemical properties of condensed matter. On the one hand, all kinds of band structure methods are based on DFT despite the fact that the interpretation of the electronic DFT levels is not rigorously justified. In contrast, the electron density obtained by the self-consistent solution of the Kohn-Sham equations is, in principle, identical to that of the corresponding interacting many-body problem (provided that the approximation made for describing exchange and correlation (xc) effects account for the actual situation of the material under investigation). As a consequence, the total energy, being a functional of the density, should be most reliable. 4 Therefore, in particular, structural properties are usually very well reproduced forming a solid basis for the success and predictive power of DFT.

Energetics including vdW interactions
For organic materials, the situation is more difficult. While the covalent intramolecular bonds can be well handled by conventional xc functionals such as the local density approximation (LDA) or the generalized gradient approximation (GGA), the intermolecular interactions are much weaker and of non-local character, i.e. governed by vdW dispersion forces. This nonlocal correlation effect has not only been a puzzle for electronic structure theory for a long time, but is a very active field of ongoing research. To review all the current methodologies, ranging from simple models to quantum chemistry approaches, would go far beyond the scope of the paper.
A simple solution to this problem are semiempirical corrections to be added on top of DFT calculations. Such methods are easy to handle and computationally efficient but often suffer from the lack of transferability. Still they are useful for obtaining a quick estimate of vdW energies, in particular for very large systems. A widely used parameterization of such a scheme has been proposed by Grimme [13,14]. While we generally stick to the DFT approach described below, we will show some results obtained with the latter method for comparison.
The most promising and successful attempts to overcome the problem of treating sparse matter within DFT, are based on the adiabatic connection dissipation fluctuation theorem (ACFDT). An approximate solution for it is provided by the random phase approximation (RPA), which is used in connection with Hartree-Fock exchange. It has been applied to a variety of systems, including molecules [15] and solids [16,17]. This scheme is, however, computationally very demanding and out of reach for systems like the molecular crystals under investigation. Here, the recently developed vdW density functional (vdW-DF) [18] represents a true alternative, being derived from the ACDFT, and hence similarly accurate (from what we can tell from the experience so far), but computationally much less involved. It describes dispersion interactions in a seamless fashion, correctly yielding the known asymptotic power law.
Since the non-local dispersion forces are usually considered to arise from non-overlapping densities, the total correlation energy is divided into two contributions [19], where E nl c [n] accounts for the long-ranged dispersive interactions. The first term, E 0 c [n], is, in principle, semilocal, but is treated by the LDA. This appears to be a reasonable approximation since it is exact in the limit of slowly varying densities. Hence, the recipe for treating correlation effects within the vdW-DF is ensuring that no double counting of non-local effects will occur. The non-local term can be considerably simplified by noting that the long-range interactions are less sensitive to the details of the system's dielectric response than the short-range terms, and hence the polarizibility is treated in a plasmon-pole model. Starting from the ACFDT formula for the correlation energy one arrives, after some tedious algebra [18], at an expression of the form 5 The details of the interaction between the electron density at the points r and r are hidden in the kernel φ, which is, by construction, a generalized function of space coordinates, the electron densities n and their gradients, i.e. φ = φ(|r − r |, n(r), n(r ), ∇n(r), ∇n(r )). Introducing the quantities where q 0 also depends on the density and its gradients, the kernel φ thus depends on r and r only through d and d , so that it can be tabulated in advance in terms of these two variables, or better, in terms of the sum and difference of these two, i.e. D and δ defined by When both d and d are large, the asymptotic form of φ(d, d ) is with C = 12(4π/9) 3 me 4 . Thus, the interaction energy has the correct r −6 dependence for large distances. Another important feature of the kernel function is that E nl c is strictly zero for systems with a uniform electron density, as was initially imposed by equation (1). Moreover, the form of equation (3) keeps the spirit of DFT in the sense that the total energy requires only the knowledge of the electron density and does not depend on any adjustable input parameter.
The theory of the vdW-DF discussed above is an approximation for the correlation energy only. In practice, one needs to approximate the exchange energy as well. The standard scheme is to use a GGA for exchange, where it is important to choose a flavor, the exchange part of which does not produce binding that is not present when the exchange is treated exactly. In fact, such a spurious binding from exchange has been found for rare gas dimers, a feature absent for the exact Hartree-Fock exchange. Since revPBE [20] does not exhibit this shortcoming, it was suggested to be used in numerical calculations.
Summarizing, the total xc energy in vdW-DF theory is defined as The vdW-DF is carried out in a post-scf manner. First, the electron density is obtained selfconsistently by using some GGA, and, in the second step, the non-local correlation energy E nl c is evaluated using this density as an input. In the final step, the vdW-DF total energy is calculated as Although in most cases PBE is taken as the functional to carry out the self-consistency cycle, it should be emphasized that the actual choice of the GGA flavor does not noticeably influence the non-local correlation energy [9], and even LDA leads to practically the same result. Finally, we note that self-consistency has little effect on the atomic interaction energy as shown by Thonhauser et al [21] for the systems investigated so far. This should be particularly true also for organic crystals where we do not expect any significant charge redistribution arising through the non-local interactions.

Optical properties
To properly account for many-body effects in the optical spectra of solids, one has to invoke many-body perturbation theory, i.e. the combination of the GW approach with the solution of the Bethe-Salpeter equation (BSE) [22]. This procedure has turned out to be a benchmark for ab initio calculations of optical properties of semiconductors and insulators. While the former gives rise to the quasi-particle band structure, the BSE accounts for the e-h interaction during the excitation process. The BSE represents the equation of motion for the e-h two-particle Green's function, expressing the linear response to an optical perturbation. Solutions of the BSE in an ab initio framework have shown that e-h interactions are indeed important in order to correctly account for bound excitons as well as a redistribution of oscillator strengths in the optical spectra of semiconductors and insulators. A variety of results are so far available for inorganic [23]- [27] and organic semiconductors [5], [28]- [34].
Starting from the BSE in its integral form it can be transformed into a matrix eigenvalue equation by expanding all quantities in terms of single-particle electron and hole states ψ ck and ψ vk , respectively [35]. From the resulting matrix eigenvalue problem the excitations energies (eigenvalues) as well as the e-h coupling coefficients (eigenvectors) can be obtained. The matrix structure of this effective e-h Hamiltonian H is given by the following expression [36]: where the diagonal blocks R and the coupling blocks C, respectively, are defined according to Here, E c 1 k 1 and E v 1 k 1 , respectively, denote the electron and hole quasi-particle energies, while v 1 c 1 k 1 ,v 2 c 2 k 2 represents the e-h interaction kernel given by the sum of the bare exchange and the screened direct interaction [22]. Note that the diagonal blocks contain the difference in the quasi-particle energy and are therefore greater than the single-particle band gap E g . Hence, the coupling matrices C can be neglected when e-h interactions are small compared to E g , which results in the so-called Tamm-Dancoff approximation (TDA). Employing the TDA reduces the size of the eigenvalue problem of the full BSE matrix H by a factor of 2 and leaves an eigenvalue problem for the Hermitian matrix R.

Implementation of vdW-DF
The non-local contribution to the correlation energy given by equation (3) implies integration in a six-dimensional space. Moreover, the long-range character of the interaction requires to consider (r, r ) pairs far beyond one unit cell. Standard integration schemes on regular three-dimensional grids would hamper the calculation of large periodic systems due to their unfortunate scaling behavior. We have thus chosen an alternative route via the Monte Carlo integration technique using the CUBA library for multidimensional integration [37]. In particular, we make use of a routine called DIVONNE, which lead to the best results in terms of accuracy, stability and timing. Details about this efficient implementation will be published Table 1. Lattice parameters for the oligoacene series. The values are given with respect to the internationally adopted reference frames [51] as used more recently for the thin film phase [48]. For comparison, the same data are provided in the original conventions, labeled as 'old'.  Table 2. Calculated cohesive energies (in eV) and surface energies (in mJ m −2 ) of the oligoacene series using vdW-DF. elsewhere [38]. In general, with the only input being the charge density, our code can be used on top of all major DFT packages. In the work shown here all self-consistent total-energy calculations were carried out with PWscf [39].

Semiempirical vdW corrections
A computationally less demanding and straightforward way to deal with vdW systems in the framework of DFT and standard approximations for the xc potentials is to add an empirical potential of the well-known asymptotic form C 6 · R −6 to the DFT energy, where R denote interatomic distances and C 6 dispersion coefficients. A systematic work along these lines was performed by Grimme coming up with two generations of such potentials, named Grimme'04 [13] and Grimme'06 [14] below. These are widely used in literature and hence we apply them for comparison with the vdW-DFT. In this approach, the vdW-corrected total energy is given by where E DFT is the energy as obtained from a self-consistent DFT run, and E disp is an empirical dispersion correction given by Here, N at is the number of atoms in the system, C i j 6 denotes the dispersion coefficient for an atom pair with indices i and j, s 6 is a global scaling factor depending on the underlying xc functional, and R i j is the corresponding interatomic distance. The last fraction represents a damping function, introduced to avoid near-singularities for small R, with R r being the sum of atomic vdW radii and d a damping factor.
The calculations carried out in this work have been performed by PBE. The vdW corrections E disp have been taken into account by adopting a cluster of 10 × 10 × 10 unit cells. This size has been found to be sufficient to give converged results. The actual parameters are taken from the original [13,14].

Cohesive properties
The total-energy calculations and structure relaxations have been performed within the DFT plane-wave package PWSCF [39]. We utilize ultrasoft pseudo-potentials [40] with a planewave energy cutoff of 40 Ryd. Bulk energies are computed by adopting the experimentally known space groups and lattice constants and relaxing the atomic positions. Isolated molecules are treated by supercells, while surface energies are calculated by the repeated slab approach. Concerning the slab thickness, only one molecular layer turned out to be sufficient to converge surface energies to within 10 meV for the low-index planes. However, for the higher-index surfaces (102) and (10-2) presented in table 3, we use two molecular layers due to the increasing complexity of the surface morphologies. Vacuum layers of 10 Å have been found to be adequate to prevent interaction between the translational images for the isolated molecules as well as slab geometries.

Optical properties
In the usual implementation of the BSE approach a coupling between the resonant and antiresonant blocks of the Hamiltonian are neglected since these terms have been found to have negligible effects on the transition energies in conventional, i.e. inorganic semiconductors such as bulk Si or GaAs. The inclusion of these terms, which proves to be relevant for the systems studied in this work poses a technical problem which is connected with the non-Hermiticity of the full Hamiltonian of equation (11). Therefore, in order to go beyond the TDA, we express the frequency-dependent macroscopic polarizability tensor following a route proposed by Schmidt et al for the TDA-BSE case [41]. Here, the macroscopic polarizability is obtained without requiring the eigenvalues of H explicitly, i.e. without the diagonalization of the non-Hermitian matrix H . Rather, it can be shown that the polarizability can be equivalently obtained from the time evolution of the polarizability and a subsequent Fourier transform [41]. In addition to this time evolution method, we have also employed a standard diagonalization scheme for the TDA-BSE Hamiltonian in order to also access the eigenvectors necessary for the evaluation of e-h wave functions.
All optical calculations presented here have been obtained by utilizing the full-potential linearized augmented plane wave (FP-LAPW) plus local orbitals method as realized in the WIEN2K code [42]. The implementation of the TDA-BSE scheme within the FP-LAPW method has been described elsewhere [27]. Convergence with respect to the k grid and number of valence and conduction bands has been checked. For all spectra we have chosen a broadening of 0.2 eV, which typically requires time steps of t ≈ 0.01 fs and around 3000 time steps within the time evolution scheme.

Oligoacene crystals
The oligoacens are typical representatives of organic semiconductors, which form molecular crystals consisting of π conjugated rod-like molecules of various sizes. These building blocks are presented in figure 1 for the series starting with naphthalene (labeled 2A), a molecule consisting of two benzene rings. The larger representatives depicted are anthracene (3A), tetracene (4A) and pentacene (5A). The last deserves special attention, since the growth of pentacene films shows a variety of crystalline phases, where three different polymorphic structures have been reported [43,44]. Only for two structures, namely the Campbell or bulk phase [45], and the so-called single crystal phase [46,47], have complete structural solutions been known. Although the thin film (TF) phase is crucial for the charge transport within thin film transistors, the geometry of this structure has been resolved only recently [48]- [50].
The smaller molecules crystallize in monoclinic space groups, only tetracene and pentacene exhibit triclinic symmetry, however, with angles close to 90 • . They all shape a herringbone pattern with herringbone angles of about 50 • . Perpendicular to the ab plane they form layers of almost upright standing molecules, as shown for pentacene in the right panel of the figure. Besides the herringbone angle, θ , two more angles, i.e. χ, and δ, are used to describe the orientation of the molecules inside the crystal. In the triclinic phase, the long molecular axes of the molecules are nearly parallel to each other, hence enclose an angle, δ, which is almost zero. The corresponding crystalline parameters are given in table 1. Most of the literature data on the whole acene series, including the bulk phase of 5A, have been given with respect to a setup, where the spacing in the b-direction is the smallest. In contrast, newer crystalline data, including the ones for 5A in the TF phase, are sticking to the so-called Niggli cell [51]. For better comparison of the two polymorphs we have transformed all the structural parameters to the new definitions. In these conventions, which will be used in the following throughout the whole paper, it becomes clear that the lattice spacings in the ab plane are very similar for all members of the series, while c naturally scales with the length of the molecule. Additionally, one can observe the general trend within the series to orient the long molecular axes along one direction. In table 1, this corresponds to a decrease of δ with the oligomer length. This behavior leads to a reduction of b while a stays almost constant, and probably is the main driving force for the change of the group symmetry from the monoclinic to triclinic one.
The pentacene thin film phase is again characterized by the familiar herringbone packing with a herringbone angle close to that of the bulk crystal phase. The most pronounced difference to the bulk phase is in the value of the tilt angle between the long molecular axis and the normal to the ab plane as can be seen in figure 2. While this angle is only about 3 • for the TF structure, it is roughly 22 • for the bulk polymorph. The main reason for this behavior can be found in the smaller interlayer separation of the bulk phase (1.45 nm) as compared with the thin film phase (1.543 nm). On the other hand, the herringbone angle θ ≈ 54 • of the latter is close to the corresponding quantity for the oligoacene series (≈ 52 • ).

Cohesive energies of oligoacenes
Until recently, one had to rely on the experimentally determined lattice parameters when calculating structural properties of molecular crystal. Based on them, the internal degrees of freedom could be relaxed, i.e. the geometry of the molecules inside the crystal as well as the orientation of the molecules with respect to the crystal axis [53,54]. In view of the many open questions related to film growth, to reliably tackle the energetics of such weakly bound systems is, however, a must. Hence, the advent of the vdW-DF as described above opened new perspectives. Here, we start to discuss the cohesive energies of the two 5A polymorphs, evaluating the quality of different xc functionals. To this extent, we also include the shorter members of the acene series, which have been described in [8]. The cohesive energy is given by where E mol and E bulk are the total energies of the isolated molecule and the bulk crystal, respectively, and n denotes the number of molecules per unit cell. With this definition, the cohesive energy is positive for any stable crystal. In figure 3, the results are presented together with available experimental values estimated from the enthalpy of sublimation [55]. Focusing first on the standard xc functionals, one can see that LDA generally underestimates cohesive energies by as much as 25%. Relying on purely local correlations, this discrepancy is as little surprising as the fact that LDA reproduces the linear trend with increasing molecular length. Through its well-known overbinding effect it mimics a compensation for the missing dispersion interactions. It is noteworthy that GGAs completely fail in describing the cohesive properties of these compounds. The observed values are almost zero for PBE, and even negative for revPBE. Thus, although LDA and GGA reveal reliable results for intermolecular conformations for fixed lattice parameters [53,54], these (semi)local functionals break down when it comes to intermolecular bonding. In contrast, the usage of the vdW-DF functional gives cohesive energies in excellent agreement with experiment, as already discussed in [8]. The right panel of figure 3 zooms into the observed range of the binding energies to better distinguish between the various results. Among the two parameterizations of semiempirical vdW contributions, Grimme'04 is even worse than LDA, whereas Grimme'06 comes closer to the ab initio results and experiment, but still underestimates the bond strength. The calculated cohesive energy of the TF polymorph is The right panel zooms into the energy range of interest, containing besides the results for 5A in the Campbell phase also those for the TF structure (marked by the yellow frame). For comparison, the data obtained from semiempirical vdW corrections by Grimme [13,14] are given.
only slightly higher than that of the Campbell structure. It suggests that these two phases can coexist in thin films and is in line with experimental observations: the TF structure is observed only for small thicknesses, while with increasing number of deposited layers the bulk phase is evolving.

Surface energies of oligoacenes
While the cohesive energies could be critically evaluated against measured data, surface energies for the molecular crystals under investigation are hardly accessible by experiment. Hence, theory plays an important role to predict these values for which, in turn, a reliable computational toolbox is required. Not surprisingly, the surface energies, again demonstrate the need for the inclusion of vdW interactions. In equation (16), E slab is the total energy of the slab, A is the surface area, and the factor two accounts for the fact that the latter contains two surfaces. In figure 4, the formation energy E is plotted as a function of vacuum distance, which is inserted between two molecular layers. The upper panels show the curves for the (100) surface, while the middle and the lower panels display the corresponding information for the (010) and (001) index planes, respectively. Before discussing the surface energies of the two phases, we evaluate the overall performance of different functionals. As described already elsewhere [8], PBE not only gives almost zero binding energy at the experimental lattice spacing, but exhibits some binding at much too large distances. Since zero vacuum size indicates the experimental lattice parameter, the aforementioned overbinding of LDA in terms of lattice spacing is visible, which at the same time underestimates the surface energies compared to vdW-DF by approximately 25-30%. In contrast, they are almost zero for PBE, which would predict unstable structures. Among the two semi-empirical potentials, Grimme'04 excellently reproduces the experimental unit cell parameters, but gives rather low values for the surface energies, while Grimme'06 improves the surface energies but worsens the equilibrium distances. In all cases, the latter behaves very similar to LDA. Since the non-local dispersion forces are solely responsible for the binding, the usage of a semi-empirical approach represents a practical tool to quickly estimate the cohesive properties of a molecular crystal, especially if it consists of a large number of atoms. In that sense, a comparison to the ab initio results could also contribute to finding optimal materialdependent parameters for the dispersion energy E disp .
The energy at infinite distance (10 Å) is set to zero, hence the respective surface energy can be obtained from the depth of the energy surface. The (100), (010), (001) and (110) index planes have been shown [56,57] to give the main contribution to the equilibrium crystal shape (ECS) for pentacene. Hence, they are displayed for the whole nA series in table 2. For monoclinic crystals, the planes (110) and (1-10) are equivalent, which is not the case for the triclinic structures 4A and 5A. Therefore, the corresponding surface energies are, in principle, different.
The calculated values for both these surfaces exhibit, however, negligibly small differences between the two phases (less than 1%) due to the small deviation of their crystal structures from the monoclinic ones. For this reason the (1-10) results are not displayed.
A common feature for all compounds is that the (001) plane exhibits the lowest surface energy. The nature of this anisotropy can be understood in terms of the crystalline packing in the (001) direction. Hence, on substrates with comparably small substrate/molecule interactions, thin films are expected to wet the substrate and to preferentially orient along the (001) direction, since the total surface free energy is minimized this way. When comparing the 5A polymorphs, there is a difference in the (001) surface energy. In fact, E 001 is larger for the TF phase than for the Campbell structure. Here one has to consider that there are two determining factors, which are the formation energy, E, as well as the surface area, A (see equation (16)). As the bottom panels of figure 4 show, the interaction energy between the (001) plans is slightly larger in the Campbell structure, which is a consequence of the smaller interlayer separation of 14.43 Å compared with 15.44 Å in the TF phase. However, the (001) surface area is smaller for the latter polymorph resulting in a higher surface energy. For the other three surfaces, the difference in the surface energies are less pronounced.
Although the four surface planes discussed above are the most important ones, the crystal shapes are not solely determined by them, as we will see below. To demonstrate this, we have also calculated the surface energies for the two pentacene phases for a larger set of indexes which are listed in table 3.
These results demonstrate that the differences between the surface energies computed by LDA and vdW-DF are no longer given by a rough proportionality as might be induced from the behavior of the low-index surfaces listed in table 2. Rather we observe that the differences are enhanced for higher index planes. The deviations may be as large as a factor of 1.8 as found for the (01-1) surface of the Campbell structure or the (102) plane for the TF phase. We assume this trend to be connected with the increasing corrugation of the latter surfaces for which interactions from more distant molecular planes start to play an important role. These interactions are better captured within the vdW-DF approach than within the LDA.

ECSs
Using the surface energies as summarized in table 2, one can obtain the ECSs based on Wulff's construction [6]. The results are visualized in figure 5. The finding that the (001) surface energy is minimal for all the oligomer series is illustrated by the fact that the (001) crystal faces have the largest area. Compared to experiment, good agreement is found with the shape of anthracene grown on graphite [58]. In particular, the (001) orientation of the crystal parallel to the substrate, the appearance of approximately hexagonally shaped facets of the (100) and (110) planes, and the small (010) facet is consistent with our calculated ECS. Going from naphthalene to anthracene, the shape is changing from octagonal to hexagonal. This finding is connected with the decreasing size of the (100) plane within the whole series, which, in turn, is caused by the increase of the (100) surface energy as discussed above.
While the overall ECSs of the two polymorphs given by LDA are very similar (top panels), more pronounced differences become visible in vdW-DF case (bottom panels). In the latter,  the (-111) and (1-11) facets dominate over (1)(2)(3)(4)(5)(6)(7)(8)(9)(10) in the thin film polymorph. Interestingly, the crystal shapes of the Campbell structure are almost the same for both functionals. This is a consequence of the fact that the LDA surface energies are underestimated by an almost constant factor. This proportionality, however, is not given for the TF phase, where the shape obtained by LDA is rather different from that produced by vdW-DF. This not only leads to important changes in the ECS of the TF phase, where (111), (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11) and (-111) start to dominate over (110) and (1-10) facets, but also highlights the importance of accounting for vdW interactions in the description of organic film morphologies.

Electronic structure
The electronic structure of pentacene has attracted considerable interest not the least due to the high mobilities reported in 5A thin films. Here, we report on the band structures of both pentacene polymorphs, the Campbell bulk phase and the thin film phase as depicted in figure 7 [48]. The band structure calculations were carried out for a path in the Brillouin zone connecting the high-symmetry points X, , Y, C and Z with the internal coordinates of these points being (0.5,0,0), (0,0,0), (0,0.5,0), (0.5,0.5,0) and (0,0,0.5) in units (2π/a,2π/b,2π/c) as indicated in the top of the figure. Note that the Z -direction is normal to the ab plane, i.e. perpendicular to the pentacene layers. Since all triclinic angles are close to 90 • , the X -and Y -directions are almost parallel to the crystal ab plane, hence, reflect the in-plane dispersion.
The subbands corresponding to the uppermost valence band (VB) arising from the highest occupied molecular orbital (HOMO) and lowest conducting band (CB) pair which are due to the lowest unoccupied molecular orbital (LUMO) are displayed in green in this figure. The band structure of the thin film phase exhibits a more dispersive character. In particular, we observe a valence bandwidth twice as large as compared to the corresponding bulk polymorph. The largest splitting of the VB and CB of the thin film phase is observed at the point, while in bulk pentacene this is at the C point. The top of the VB of the thin film polymorph is situated at (0.7, 0.0, 0.0). It is noteworthy that the Kohn-Sham band gap of the thin film phase (0.70 eV) is only slightly smaller than that of the bulk phase (0.74 eV), for which similar values have been reported previously [33,59]. In both cases, the gap is a direct one.
The denser packing of the molecules in the ab plane leads to stronger intermolecular interaction. This is illustrated in of figure 8 which shows the square of the HOMO wave function at the point in this plane. One notices a stronger overlap of the wave functions in the TF structure compared to the bulk phase. As a result, the bandwidths of both the conduction and the valence bands are significantly enhanced in the TF phase compared to those of the bulk polymorph. This finding has important consequences for the electron and hole mobility in pentacene organic field effect transistors. It shows that the growth of the thin film phase close to the gate dielectric is really crucial in order to maximize the performance of such devices.
We conclude this section by noting that the small dispersion observed in the bulk phase is in line with recent angular resolved photoemission experiments on highly ordered films of pentacene, which have been clearly identified as belonging to the single crystal phase [60]. In this work, it was shown that the measured band dispersion of the HOMO agrees both in magnitude and direction with calculations of the band structure for the single crystal phase. The higher band dispersion observed by other groups [61,62] might indeed be due to the presence of the thin film phase though other explanations related to film morphology have also been offered [63].

Optical properties
The calculated optical spectra of pentacene are depicted for the two polymorph phases in figure 2. Thereby the imaginary part of the frequency-dependent dielectric tensor is obtained by solving the BSE for the e-h pair according to equations (11) and (12). By successively improving the e-h Hamiltonian H , we are able to assess the influence of commonly employed approximations leading to four levels of sophistication. Firstly, by setting the e-h kernel equal to zero the electron and the hole are regarded as non-interacting leading to the independent partial approximation or random phase approximation (RPA). Taking into account the e-h exchange interaction in is equivalent to the RPA including local field effects (RPA + LF). As a third step, we add the direct, statically screened, e-h interaction to the kernel (TDA-BSE). It is this latter term, which can lead to the formation of bound excitons. Finally, we also take into account the coupling between resonant and anti-resonant terms in the Hamiltonain, that is we go beyond the TDA. Optical spectra for these four levels of approximation are compiled in figure 9 for the pentacene bulk crystal phase (top row) and the pentacene thin film phase (bottom row). Due to the triclinic crystal structure there are six independent components of the dielectric tensor of which we show the three diagonal ones labeled as x x (left column), yy (middle column) and zz (right column), respectively. It is apparent that the anisotropy of the crystal structure is also reflected in the optical properties. The energetically lowest transition of pentacene, appearing around 2 eV in the RPA is polarized parallel to the short molecular axis, followed by a stronger transition around 4 eV in RPA polarized parallel to the long molecular axis. Due to the distinct orientations of the pentacene molecules with respect to the molecular layers, i.e. almost perpendicular in the thin film phase and inclined by 22 • in the Campbell phase, only traces of the 4 eV transition can be found in the x x and yy components of the thin film phase, whereas it appears as rather strong contribution to x x and yy in the bulk crystal phase. This argument not only holds for the spectra computed within the RPA but also for the higher levels of approximation. As a consequence one expects the polymorphism in pentacene thin films consisting of upright standing molecules to be detectable by optical spectroscopy. While the 4 eV peak should clearly show up in normal incidence reflectance spectroscopy in the bulk crystal phase, only a weak signature of this transition should be found in the TF structure.
When turning to the influence of the various approximations in the e-h kernel mentioned above one observes a rather strong effect of the e-h exchange interaction. Due to its repulsive character it shifts the RPA spectrum (filled gray area) to higher energies and also modifies the intensities of transitions rather strongly with the tendency to transfer oscillator strengths to higher frequencies. The direct e-h interaction has the opposite effect: it is attractive and leads to red-shifted spectra and transferral of oscillator strengths to lower energies. For the zz Figure 10. The e-h wave function as isosurface plots (yellow indicating positive, blue negative values) of the energetically lowest exciton for pentacene in the bulk (top) and TF phase (bottom). The left column shows side views on the herring bone layers, while the right column displays a top view highlighting the spatial extent with one layer. The hole coordinate has been fixed at the position indicated by the green cross, the unit cell vectors are shown as red arrows.
polarized spectrum the effects of these competing terms in the e-h kernel almost compensate each other resulting in an optical spectra (TDA-BSE, red line) rather similar to the original RPA spectrum with peaks slightly above the RPA ones. On the other hand, in the case of the lowest energy transition appearing in the x x and yy components in the bulk dielectric function, the inclusion of the attractive e-h interaction produces a bound exciton state within the band gap seen as small peak in the x x and yy components around 1.6 eV. However, unlike the situation in phenylenes [64] and thiophenes [31] the oscillator strength of this bound exciton peak decreases compared to the RPA result where the effect is more pronounced in the x x component and comparably small in yy.
We analyze these bound exciton states in more detail by evaluating the corresponding e-h wave function which can be constructed from the eigenvectors of the BSE Hamiltonian and the Kohn-Sham valence and conduction states [27]. They are displayed in figure 10 for both, the bulk crystal phase (top row) and the thin film phase (bottom row). Two different projections are shown, a side view emphasizing the layered crystal structure (left) and a top view bringing forward the herringbone arrangement of the molecules within the layers (right column). Note that an e-h wave function generally contains two spatial coordinates, one for the electron and one for the hole. In order to visualize such a quantity it is common practice to fix one of the spatial coordinates to a specific point in space. We choose to fix the hole coordinate in the middle of a pentacene molecule at the center of the π cloud slightly above the molecular plane indicated by the green crosses in figure 10. The electron distribution with respect to this fixed hole position is shown as isosurface plot, where the yellow (blue) color indicates positive (negative) sign in the wave function. From the plots it is apparent that the lowest exciton in pentacene is well localized within one molecular layer. A top view on the layers reveals that the exciton is considerably extended within one layer resulting in a two-dimensional character of the exciton wave function. We find the exciton to be rather isotropic in the ab plane for the bulk crystal phase while there is a pronounced anisotropy for the thin film phase with a larger extent of the exciton in the crystal a direction.
Finally, we turn to the effect of the resonant/anti-resonant coupling terms in the BSE Hamiltonian depicted by the orange curves in figure 9. The overall effect of the off-diagonal coupling matrices C in equation (11) is a decrease of the transition energies, which can be understood from the matrix structure of the Hamiltonian. However, the magnitude of this effect has different strength for the three polarization directions. While the main transition for zz polarization is red-shifted by roughly 0.3 eV with respect to the TDA-BSE, there is hardly any change in the energetic position of the lowest exciton peak observed in the x x and yy components of the dielectric tensor.

Conclusions
In summary, we have investigated the energetics of the oligoacene series in the framework of DFT including vdW interactions. This approach enables us to tackle the prediction of surface energies for these materials, which so far have been inaccessible by experimental techniques. Although we confirm that, as in previous studies, LDA systematically underestimates these quantities by roughly 25-30% for the simplest surfaces as (100), (010), (001) and (110), we find that the discrepancies become more dramatic for other surface planes. This fact has also implications on the crystal shapes, which we obtain by Wulff's construction. Taking into account only the four low index planes mentioned above, the polyhedra change from octagonal in case of naphthalene to hexagonal for the longer molecules. However, we also show that these shapes are strongly influenced by the higher index planes, revealing distinctively more facets, where their size depends on the underlying xc functional. An even more pronounced change in crystal shape is found when adopting the thin film phase of pentacene. The reasons are traced back to the different molecular conformation compared with the Campbell structure. Besides the cohesive properties of these two polymorphs, we also investigate the electronic and correspondingly the optical properties. By comparison of the band structures we can understand the higher mobilities found in the thin film phase, but also specific features in the dielectric tensors. The latter are calculated by solving the BSE to account for excitonic effects, where the coupling between resonant and anti-resonant terms turn out to be important to correctly describe the optical spectra. The e-h interactions play a dominant role in both phases when the light is polarized perpendicular to the long molecular axis. For these in-plane components, we predict a pronounced difference between the two structures, which is a direct consequence of their distinction in terms of the angle between the molecular axis and the surface normal, which is 22 • in the Campbell phase and 3 • only in the TF structure. In fact, due to the triclinic symmetry, the zz component shows up in the other diagonal terms. Hence, the presence or absence of the Campbell phase in thin films should by easily detectable in optical spectroscopy.