Ab initio study of electron-phonon coupling in rubrene

The use of ab initio methods for accurate simulations of electronic, phononic, and electron-phonon properties of molecular materials such as organic crystals is a challenge that is often tackled stepwise based on molecular properties calculated in gas phase and perturbatively treated parameters relevant for solid phases. In contrast, in this work we report a full first-principles description of such properties for the prototypical rubrene crystals. More specifically, we determine a Holstein-Peierls--type Hamiltonian for rubrene, including local and nonlocal electron-phonon couplings. Thereby, a recipe for circumventing the issue of numerical inaccuracies with low-frequency phonons is presented. In addition, we study the phenyl group motion with a molecular dynamics approach.


I. INTRODUCTION
The scientific interest in organic semiconductors is continuously high and driven by many technological applications, where the molecular materials are employed in organic transistors, 1-5 organic light emitting diodes [6][7][8][9] , photovoltaic applications [10][11][12][13] and a variety of other electronic devices [14][15][16][17] . Thereby good charge transport properties of organic semiconductors is the key for their efficient application and a great effort is dedicated to improve carrier mobilities by means of chemical and structural modifications of these organic materials. In this process, theoretical input can provide guidelines towards high mobilities of charge carriers and additional functionality of the organic semiconductors. The understanding of several properties of these materials, however, remains incomplete as charge transport usually shows very different behaviour compared to conventional inorganic semiconductors. In addition, a complete quantitative characterization of the properties of the organic materials remains a challenge due to their complexity.
On the theoretical side, rubrene has first been studied as gas-phase molecules and later in solid state by means of semi-empirical and higher-level approaches 41,42 . However, to date a full ab initio characterization of the electronic properties and electron-phonon coupling (including Holstein and Peierls type of couplings) has not been achieved, which might be due to the sizeable structure of the rubrene unit cell consisting of 280 atoms and resulting in 840 vibrational modes. As an additional challenge, we note that standard density-functional theory (DFT) based methods have well-known difficulties with numerical accuracy when describing low-frequency vibrations.
Here we present a full first-principles characterization of rubrene crystals including all intra-molecular and intermolecular modes. We thereby demonstrate a practical way to remove the problem of inaccuracies with lowfrequency modes and imaginary frequencies that often occur in DFT-based methods. This enables their application also for systems with a large number of atoms.
The paper is organized as follows. In Sect. II we introduce theoretical and computational methods. Sect. III contains the resulting electronic structure, vibrational properties and electron-phonon couplings. Finally, we analyze the phenyl-group motions and possible flipping in a molecular dynamics (MD) study.

II. METHODOLOGY
A. Ab initio total energy approach We perform density functional theory 43 calculations of rubrene in crystal geometry. Our unit cell consists of four molecules of rubrene where pairs of molecules are related by a screw operation and other rubrene pairs by a non-primitive translation (see Fig. 1). Note that, while being computationally more demanding, a large unit cell has the advantage of weaker variations in the smaller Brillouin zone. We start from the experimental coordinates 20 , and perform a conjugate gradient optimization to obtain relaxed atomic coordinates and lattice constants of the crystal. For this relaxed unit cell, we calculate electronic properties, phonon modes and frequencies as well as electron-phonon coupling parameters.
All the simulations are performed with the Siesta code. 44,45 The calculations are done using the Local Density Approximation (LDA) using the exchangecorrelation (XC) potential of Ceperley and Alder 46 as parametrized by Perdew and Zunger. 47 It is known that intermolecular bonding in molecular crystals and simi- lar systems is partly due to van der Waals interaction which is only included in the homogenous limit in the semi-local XC functionals such as LDA or GGA-type functionals. [48][49][50] Still, the LDA approach used here has been shown to give reliable results for organic crystals 51 which is also the case for rubrene (see below).
Separable, norm-conserving pseudopotentials of the Troullier-Martins type 52 in the Kleinman-Bylander form, 53 are used to describe the effect of the core electrons. The basis sets in Siesta are strictly localized numerical atomic orbitals. We have used a split-valence double-ζ basis set including polarization functions, optimized for the bulk structure of rubrene. The parameters that define the basis are presented in Table I. Phonon calculations are performed using a finite differences approach, as described in Section II B. Molecular dynamics simulations are performed in the constant temperature ensemble using a Nose-Hoover thermostat and a time step of 0.5 fs.

Method I : direct diagonalization
For the simulation of the vibrational properties we calculate the force-constant matrix by displacing individual atoms along the Cartesian directions. The dynamical matrix is defined with the Hessian of the DFT total en-ergy according to where m i and m j are the masses of atoms i and j. u iα and u jβ are the displacements of these atoms along the Cartesian coordinates α and β, respectively. From Newton's equations of motion, the force F i on atom i upon displacement by u i is given as F i = − ∂E ∂ui . By approximation, we take finite differences with positive and negative displacements such that the dynamical matrix is calculated as The solution of the eigenvalue equation of the dynamical matrix yields eigenvectors e p (mode index p) with the Cartesian components e p iα (for atom i) and phonon frequencies ω p . The phonon normal modes ζ p are calculated as

Method II : frozen phonons
In order to obtain the phonon frequencies to a sufficient level of accuracy and minimize numerical inaccuracies (see below), we also employ an alternative method for their calculation, which improves the results of the calculations for the low-frequency modes. We take a two-step approach which involves the above direct diagonalization as a first step. As a second step, we displace all atoms of the crystal in the direction of the mode vectors ζ p of the given phonon mode p according to and re-calculate the total energies of the crystal with the frozen phonon. We perform these calculations for every normal mode for a set of different amplitudes λ and observe the change in energy E → E + ∆E(λ) which is given by the quadratic form Given the relaxed ground state for λ = 0, this allows to compute low-energy mode frequencies with better accuracy and without the problem of imaginary or negative ones as will be analyzed in Section III.B. Note that formally both approaches should give the same results if numerical inaccuracies were not present. Beyond the numerical calculation of electronic energy bands and vibrational properties, we derive an effective tight-binding model to parametrize electronic properties and electron-phonon coupling interactions in rubrene. Thereby we focus on the states derived from the HOMO of rubrene, i.e. the valence band structure in the crystal. Such material parameters are essential for simulating p-type charge transport.
We use a Holstein-Peierls model with the Hamiltonian which consists of an electronic part, a phononic part and a coupling part between electrons and phonons and explicitly reads with ε M N the transfer integrals between HOMO states M and N. a M (a † M ) and b Q (b † Q ) are the annihilation (creation) operators for electrons and phonons, respectively. Q is the coordinate of the phonon Q = (q, p) where p is the mode index and q is the wave vector. This model accounts both for intra-molecular (onsite/local) and intermolecular (non-local) electron-phonon interaction effects by linear coupling to the phonon operators (see below).
The electronic Hamiltonian (first term in Eq. (7)) depends on the full set of transfer integrals ε M N of molecules M and N. Due to the high symmetry of the unit cell, however, this set can be reduced such that the TBmodel of rubrene requires only few relevant transfer integrals. According to the assignment of indices (A to D) to the molecules in Fig. 2, the remaining symmetry-reduced electronic transfer integrals are ε AC , ε AD , ε AB , ε AA±b , and ε AA±2b . Here, b indicates a lattice vector in vertical direction and A + b denotes the orbital A in the neighbor unit cell. In order to find an analytically tractable form of the band structure, we simplify the Bloch-Hamitonian by assuming negligible coupling between molecules A and D (see results section), i.e. ε AD = 0). Taking into account the remaining terms, the band energies in k-space have the form Each of the four sign combinations (++, +−, −+, −−) in the second and third line of Eq. (8) gives rise to a band in the rubrene band structure. The set of transfer integrals will be determined by a least-squares fit to the ab initio band structure.

D. Electron-phonon coupling in rubrene
The electron-phonon coupling constants g Q M N ≡ g p M N (q) in the model are defined as the linear changes of the electronic matrix elements ε M N with the amplitudes λ of the phonon normal mode p with wave vector q and are calculated as In this definition, the coupling constants are dimensionless and explicitly depend on the phonon wave vector q. Note that a slightly different convention is used in a recent review. 55 Given the large supercell and the small Brillouin zone, this wave vector dependence, if necessary, can be captured by a simple model and does not need to be calculated for the large amount of modes. We are thus only interested in the ab initio fits for the (q = 0) case, which effectively reduces the computational effort.
To simplify the notation in Eq. (9), we introduce the dimensionless Holstein couplings constants according to which are averaged over the 4 molecules in the unit cell. Analogously, the non-local Peierls coupling constants are defined with i ∈ {AA ± b, AB, AC}. In consistency with the neglect of the transfer interal ε AD , we set the corresponding coupling constant g AD = 0. For similar reasons, g AA+2b which is much smaller (next nearest neighbor according to the lattice vector 2b) will not be considered for simplicity. In order to extract the full set of the above-defined electron-phonon couplings with ab initio methods, we perform DFT calculations in the frozen phonon geometry (Eq. (3)). Clearly, it would be too a formidable task to re-fit the band structure for all vibrations as was done for the ground state (λ = 0). We therefore take another route based on the Kohn-Sham Hamiltonian. We employ a basis transformation from the atomic orbital basis used in the Siesta calculations into the desired basis set of molecular orbitals. The Kohn-Sham HamiltonianĤ in the basis of atomic orbitals is described as The transition from the basis sets used in DFT and the desired orthogonal TB-model, can be performed in two steps. First, we project onto the molecular HOMO orbitals |ψ HOMO M known from gas-phase calculations where the index M indicates the molecule on which the HOMO orbital is located. |ϕ M µ are the basis functions associated to M and the coefficients a µM are easily obtained from Siesta calculations of single molecules. This allows expressing the Hamiltonian in the non-orthogonal HOMO basis while neglecting other types of molecular orbitals. This first step results in a Hamiltonian whose matrix elements are given as is the molecular overlap matrix. Coupling of the HOMO orbitals to other molecular states is neglected as we are only interested in the parameters of the HOMO-derived bands due to their energetic separation to other bands. The second step is the orthogonalization of the molecular orbitals for which we chose the Löwdin orthogonalization method 56 with the overlap matrix (12). We generate a new set of HOMO-like orbitals |Ψ M from the given initial set of normalized but non-orthogonal wave functions as This single step mixing can be repeated iteratively until the wave functions are orthogonal to the sufficient degree of accuracy. In the present case for the HOMO orbitals of rubrene, a single step of orthogonalization turns out to be sufficient as already the first iteration leads to negligibly small overlap matrices. This orthogonalization procedure, applied for undisplaced phonons, leads to the effective Hamiltonian H el in Eq. (6). By applying this approach for structures including the frozen-phonon distortion, we determine, by finite differences, the linear changes in the transfer integrals and obtain the Holstein and Peierls electron-phonon coupling constants in (7). Consequently, this transformation yields the coupling constants g Q M N in the orthogonal molecular HOMO basis.

A. Electronic band structure and transfer integrals
We first focus on the electronic properties of rubrene. The DFT band structure together with the TB-model  Fig. 3. The fit to the band structure results in the transfer integrals compiled in Tab. II, which reproduce the DFT band structure with good accuracy (rms deviation of 4.0 meV on a regular k-point grid). A comparison of the fits with and without the second neighbor term ε AA±2b in b direction shows slightly better results when taking this term into account (see inset of Fig. 3). This small change is consistent with an order-of-magnitude smaller value of ε AA±2b compared to ε AA±b . This additional parameter corrects the underestimation of the DFT-bands in the ΓY -direction.
Very similar values for the transfer integrals are obtained directly from the Kohn-Sham Hamiltonian (described in Section II D) by using the Löwdin orthogonalization of the HOMO bands. Both methods deviate by only about 3 − 4 meV and both sets of transfer integrals are consistent with the experimental band dispersion in the direction measured by ARPES experiments 57,58 and with theoretical calculations reported in literature 29,59 as is summarized in Tab II.

B. Vibrational properties
We turn to the calculation of the dynamical properties. By using Method I, we obtain all mode frequencies and eigenvectors. A selected list of low-frequency modes with ω (I) < 50 cm −1 summarized in Tab. III. It turns out that for some of the modes the standard diagonalization scheme leads to some unrealistic small or even imaginary frequencies. This general problem is due to residual numerical inaccuracies which may severely jeopardize the further calculation of electron-phonon coupling constants (Eqs. (10) and (11)), which requires division by the frequency. The reason for this behaviour originates from slightly noisy force constants that lead to the difficulties in the description of collective crystal vibrations in soft materials like organic crystals. These vibrations involve the motion of many atoms (at least for collective molecular modes such as libration modes or translations) which means that many force constants enter in the resulting vibrational frequency, as opposed to e.g. a CH-stretch mode, where the force constant related to the CH bond will mostly define the frequency. Considering that the forces ∆F iα for displaced atom i may have some small numerical error, this error can accumulate in a quantity like ω p that results essentially from all forces. If the mode frequency is low, the relative error can then be very large for such modes.
On the other hand, by inspecting all low-frequency modes in rubrene, we find that they show correct vibration patterns (e.g. for translation and libration modes) and are also orthogonal, which gives us confidence in the mode vectors ζ in contrast to the vibration frequency. Therefore we apply them as frozen modes in Method II for the re-calculation of the mode frequencies ω (II) and re-calculated the changes in the total energy ∆E around the equilibrium configuration for different normal mode amplitudes λ according to Eq. (5) for all vibrations. Fig. 4 shows an example to illustrate the difference of both methods for a particular mode (main frame). This specific vibration is an A g symmetric torsional mode of the phenyl rings. The energy obtained from the direct diagonalization of D iαjβ is ω (I) = 5.5 meV (44.4 cm −1 ), while the quadratic fit to the total energy in Method II yields the corrected mode energy ω (II) = 7.4 meV (59.6 cm −1 ). This value is increased by about 34 % compared to the original frequency and is a typical example of the general behaviour of low-frequency modes. Table  III compares a larger set of low-frequency modes in the two approaches and shows how strongly the frequencies can be corrected by Method II. A full overview over all modes (except the highest-frequency CH-stretch modes) is provided in the inset of Fig. 4. In this figure we plot  the relative frequency differences Apparently, with Method II we obtain a trend towards systematically higher energies especially for the lowfrequency modes which are corrected by up to 50 % (with few exceptions of even higher values). The median value for ∆ω rel is increased by 0.724 % indicating that the upper half of the frequency spectrum experiences only slight changes. This indicates that mainly collective modes are effected and should be corrected regarding their frequency, whereas ω (I) ≈ ω (II) for high-frequency modes, whose frequencies are described well by the initial diagonalization of the dynamical matrix. In the following we will only use the corrected frequencies ω (II) for the calculation of the electron-phonon coupling parameters.

C. Electron-phonon coupling
The calculation of the electron-phonon coupling in rubrene is performed with the frozen phonon approach. The effect of a specific mode ζ p on the band structure is visualized in Fig. 5. The HOMO bands are shifted proportional to the amplitude λ. This example illustrates the impact of a vibration of frequency 1593.3 cm −1 , for which we found a strong impact on the HOMO bands. This mode is an intramolecular C-C stretch mode and as such strongly changes the onsite energy. The changes in the band structure for different amplitudes λ would in principle enable the extraction of the electron-phonon coupling constants from the linear slope of the changes in the electronic energies (see e.g. inset of Fig. 5) 60,61 . In the considered case, the averaged electron-phonon coupling constant for the Γ point results in g 1593.3 0 = 0.21.
The practical calculation of the electron-phonon couplings of all phonon modes however is performed in the way described in Sect. II D, which avoids manual fits to the band structure for all modes. In order to present an overview over the total strength of the Holstein and Peierls coupling for all modes, we define the polaron bind- ing and its energy contributions E H and E P according to where E p H and E p P are the mode-resolved Holstein (H) and Peierls (P) contributions and E P consists of all Peierlscoupling constants for all nearest neighbors of a specific molecule.
The large number of electron-phonon coupling parameters is summarized in Fig. 6 by plotting E p H (red bars) and E p P (blue bars). Although the majority of the 840 modes have vanishing electron-phonon coupling, we identify a broad distribution of the remaining coupling modes over the plotted frequency spectrum. Modes with mainly non-local couplings dominate the low-frequency part (< 150 cm −1 ) with E p P values of up to 5 meV from strongly coupled B 3g -modes. Local couplings dominate the high-frequency part with a few strong coupling modes around 1600 cm −1 . This is consistent with gas-phase molecule simulations where the largest values of E p H are 15.45 meV, 7.24 meV, and 3 meV for the vibrations at 1588 cm −1 , 1007 cm −1 and 602 cm −1 , respectively (green bars in Fig. 6). We also observe at low energies that the strong-coupling modes from gas-phase are splitted in energy in the crystal phase and the distribution of the respective electron-phonon coupling is strongly broadened. A selected list is contained in Tab. IV and compared to literature values.
From our results only inversion symmetric modes with either the A g or B 3g symmetry contribute to g p H and affect the onsite energies of the molecules. Our data suggest that anti-symmetric modes in general play a subordinate role regarding electron-phonon interaction in rubrene in accordance to general symmetry arguments.
The larger effective Peierls mode frequency ω P is a results from the observed difference of our low-frequency electron-phonon coupling constants and the reference values. From the given mode patterns of relevant modes, we observe that the inter-molecular electron-phonon coupling is associated with the motion of the phenyl rings.
In particular, modes with torsion of the phenyl rings contribute here (e.g. mode from fig. 4 with ω p =59.6 cm −1 with dimensionless coupling constants of g p 0 = 0.19 and g p AA±b = 0.12). On the other hand, phenyl ring wagging (or flipping) modes that move perpendicular to the tetracene plane, couple only weakly to the HOMO. The impact of those flipping motions, which has been discussed in literature 62 is investigated in the following section.

D. Flipping motion of phenyl groups
In this section we analyze if the dynamics of the phenyl groups may be stronger than suggested by the weak electron-phonon coupling associated to these vibrations. Motivated by the work of Kloc et al. 62 , we performed MD simulations as described in Sect. II A. The authors of Ref. 62 suggested a model where two phenyl groups from the same side of rubrene change their positions relative to the tetracene backbone from above to below the tetracene plane and vice versa. Indeed, phenyl groups are very flexible and they vibrate around their equilibrium positions which are either below and above the backbone plane on each side due to mutual repulsion. We investigate this situation by calculations with a supercell, which is twice the original cell (i.e., consisting of 8 molecules of rubrene). For one of the eight molecules in the supercell we changed positions of the phenyl rings to the other side of the tetracene backbone. We let the structure relax its energy in the DFT simulations and observed that this starting configuration is not a stable minimum. Indeed the phenyl group moved back to their original position of the undistorted crystal.
In order to analyze the impact of finite temperature on the question of the phenyl group dynamics, we studied their motions with MD simulations. We expected that with increasing temperature the phenyl groups would vibrate so much that they would be able flip to the other side of the tetracene plane. To analyze such events in these simulation, we define the bonding angles α 1 and α 2 according to Fig. 7 and follow the change of the angles over time for various temperatures up to 500 K. Figure 8 shows the dynamics of the phenyl groups for the largest temperature and all eight molecules. While the phenyl groups are vibrating strongly, there are no changes in the sign of angles which would indicate that the phenyl groups flipped to the other side of the tetracene plane. The fact that we do not observe such flipping even at such large temperature, suggests that they are likely unable to cross the backbone plane at am-   Fig. 8 but with weaker oscillations in the dynamics of α 1 and α 2 were obtained (not shown as plot), thus confirming the qualitative results. The analysis of the average angles by interpreting the MD trajectory as statistical ensemble and averaging over α 1 and α 2 , is summarized in Tab. V for all temperatures. While the average angle is rather independent of T , its standard deviation increases with T and reaches 22% at 300 K. Finally, we note that we do not discard that these flipping processes happen at larger time scales, but we do not observe them in the time scale of one picosecond, and we never found such flipping configurations to be a (meta)stable minimum.

IV. CONCLUSIONS
A protocol for reliable but efficient calculations of material parameters from DFT has been introduced and applied for rubrene. In particular, we have used DFT-based methods to compute all relevant electronic, phononic and electron-phonon interaction parameters of rubrene, which serves as a prototype of a complex organic material build from a molecular core and functional side groups. More specifically, we have determined tightbinding parameters for a Holstein-Peierls type Hamiltonian for electron-phonon coupling which is directly applicable for charge transport modeling, or more generally for studies of macroscopic material specific properties, such as density of states, conductivity, carrier mobility, etc. The comparison of the parameters with other work, shows good agreement. The accurate treatment of electrons and phonons in organic small-molecule systems together with sophisticated macro-scale simulations is the basis for advancing the level of transport modeling in these systems.