Discovering structure–property relationships for the phonon band structures of hydrocarbon-based organic semiconductor crystals: the instructive case of acenes

Abstract

By studying the low-frequency phonon bands of a series of crystalline acenes, this article lays the foundation for the development of structure–property relationships for phonons in organic semiconductors. Combining state-of-the art quantum–mechanical simulations with simple classical models, we explain how and why phonon frequencies and group velocities do or do not change when varying the molecular and crystal structures of the materials.

---

Introduction

Organic semiconductors (OSCs) are a scientifically and technologically appealing class of materials, which allow the realisation of cheap electronic systems with specifically tuned properties, and which exploit the incredibly versatile toolbox of organic chemistry.^1,2 Devices built from such materials include organic light-emitting diodes,³ solar cells,^4,5 or thin-film transistors.⁶ Many of the physical properties of OSCs are crucially influenced by phonons, which are the quasiparticles of lattice vibrations. They impact thermal properties, such as the heat capacity or free energy, including the phonon contribution to a material's entropy. These quantities are crucial for assessing the relative thermodynamic stability of polymorphs,^7–10 which are often energetically very close.^11,12 Moreover, essentially all practically relevant transport properties of an OSC crystal are determined by phonons. This applies to the crystal's thermal conductivity,^13–15 its ability to convert heat to electric energy (i.e., its thermoelectric performance),^16–18 and, via the particularly strong electron–phonon coupling, its charge-carrier mobility.^19–31 For the latter, certain killer phonon modes have already been identified in the literature.²⁹ Due to their particularly strong coupling with parameters relevant to charge transport (like the intermolecular transfer integrals) they have been identified as the key reason for charge-carrier mobilities in most OSC materials being much lower than expected.²⁹

Consequently, a precise understanding of phonons in OSCs is key to understanding the above-mentioned properties. All harmonic characteristics of phonons are encoded in the respective phonon band structures. Of particular relevance in this context is the development of a clear picture of how phonon properties and band structures are affected by aspects like the molecular packing and the structure, length, and shape of the molecules that constitute the OSC. Unfortunately, so far, comprehensive structure–property relationships for phonons in OSCs are largely missing.

In order to provide a starting point for their development, in the following, we present a detailed discussion of phonons of a series of acenes, which are arguably one of the most central families of organic semiconducting materials. Unfortunately, an experimental determination of phonon band structures is extremely challenging for these materials, as for the commonly employed neutron scattering experiments, one requires comparatively large single crystals, which also need to be fully deuterated to reduce the fraction of incoherent scattering events. Consequently, there are only very few, scattered, examples of crystals consisting of conjugated molecules for which phonon band structures have been measured. These comprise, for example, deuterated naphthalene^32,33 and anthracene,³⁴ for which the dependence of the bands on temperature,^32,35–37 pressure,³⁸ and on the anharmonic nature of the phonons³⁹ have been determined. Thus, for our discussion we will rely on suitably benchmarked dispersion-corrected density-functional-theory (DFT) calculations. These have the additional benefit of providing direct access to the atomistic motions that describe specific phonon modes. The methodology employed here can be used with great confidence, bearing in mind that it shows excellent agreement with the experimentally determined phonon bands of the above-mentioned deuterated naphthalene.^40–42 Most challenging for the simulations is the accurate description of the intermolecular phonon region (up to ∼6 THz ≈ 200 cm⁻¹), as it is strongly impacted by intermolecular van der Waals (vdW) interactions that are not properly accounted for in semi-local DFT. Therefore, we benchmarked a variety of vdW corrections against available experimental Raman data⁴³ and phonon band structures⁴⁰ for OSCs in that region. For the corrections applied here, RMS deviations between measured and simulated phonon frequencies on the order of ∼0.13 THz (∼4.3 cm⁻¹) have been obtained.^40,43 Important ingredients for phonon calculations are also the crystal structures of the studied materials, where we started from (existing) OSC crystal structures documented in the Cambridge Structural Database (see below).⁴⁴

In the following we provide a systematic discussion of the phonon properties of a series of acenes from benzene to pentacene. This involves characterisation of the various phonon modes, especially with respect to their evolution with molecular length. Also here, the primary focus will be on the low-frequency region (≤7 THz ≈ 230 cm⁻¹), which comprises all intermolecular phonons and the lowest intramolecular bands, and which is particularly relevant for transport processes, as the corresponding modes are thermally occupied.⁴⁰ The observed modes comprise molecular translations (the killer phonon modes) and rotations as well as molecular bending and torsional motions. In addition, we analyse the group velocities of the phonons, concentrating primarily on the acoustic bands (i.e., the sound velocities in the materials, which can reach several kilometres per second), and finally we will discuss the impact of the phonon properties on the temperature dependence of the heat capacity as an illustrative example of a thermodynamic property. In the spirit of structure–property relationships, our results show that the evolution with molecular length differs considerably between different phonon modes (i.e., different types of vibrations). In the majority of cases the evolution can, however, still be explained by simple classical models building on the changes in molecular masses, moments of inertia, unit-cell sizes, and estimated intermolecular force constants. These “molecule-based” trends are then modified by the differences in the packing motifs of the different acenes.

Methods

All phonon band structures were calculated using the Phonopy package⁴⁵ (version 1.13.2) based on atomic forces calculated within density functional theory^46,47 (DFT) employing the Vienna ab initio simulation package^48–51 (VASP; version 5.4.1). DFT calculations were carried out using the PBE functional,⁵² the standard pseudopotentials⁵³ for the projector-augmented wave method,⁵⁴ and a plane wave energy cutoff of 900 eV. Given its outstanding performance compared with other a posteriori van der Waals corrections for calculating the properties of low-frequency phonons in organic semiconductors,^40–43 we employed the D3 van der Waals correction⁵⁵ with Becke–Johnson damping⁵⁶ (D3–BJ).

Prior to phonon calculations, the geometries of the various systems originally obtained from the Cambridge Structural Database⁴⁴ were optimised in terms of their atomic positions and lattice parameters based on the Rose–Vinet equation of state⁵⁷ to avoid Pulay stresses,⁵⁸ as described in more detail in ref. 40. Phonon band structures were calculated for suitably converged supercells and using a finite displacement approach (displacement distance of 0.01 Å). Further details regarding the employed simulation parameters, the meshes for the discrete q-point sampling and the calculation of densities of states per group velocity and frequency can be found in Section S4 (ESI†). The produced output and the required input to reproduce our results can be accessed via the NOMAD database.⁵⁹

The phonon displacement patterns and the animations in the ESI† were generated by visualising the geometries displaced along the normal mode coordinates using the Ovito software (version 3.3.1).⁶⁰ The structural view of pentacene in Fig. 3(c) was plotted using Vesta (version 3.5.5).⁶¹

Results and discussion

Structures of the studied systems

All acenes discussed below crystallize in a herringbone pattern with two molecules per unit cell (Z = 2). The molecules are arranged in planes parallel to the a₁ and a₂ axes. The choice of the order is such that a₁ (a₂) denotes the longer (shorter) distance of symmetry equivalent molecules, while the a₃ axis is roughly parallel to the long molecular axis (for more details see Section S3, ESI†). Benzene⁶² (1A), naphthalene⁶³ (2A) and anthracene⁶⁴ (3A) crystallize in monoclinic lattices, while crystals of tetracene⁶⁵ (4A) and the two considered pentacene polymorphs^66,67 (5A and 5A-II) belong to the triclinic crystal class. In passing we also note that an orthorhombic benzene polymorph⁶⁸ (1A-o) exists, which is thermodynamically more stable at low temperatures (and zero pressure).^69,70 Due to its fundamentally different structure, its properties are described only in Section S9 (ESI†). Moreover, although at room temperature polymorph 5A-II is more stable than 5A, we will primarily focus on 5A as again its structure is more consistent with the structures of the shorter acenes, as explained in Section S3 (ESI†). The latter also contains an in-depth comparison of the lattice parameters and the arrangements of the two molecules per unit cell in the different acenes. This comparison shows that there are two groups of systems: in 1A, 2A, and 3A, the two molecules per unit cell tilt symmetrically into directions above and below the a₁ axis (space group number 14); conversely, in the longer acenes, both molecules tilt essentially in the same direction such that their glide planes and screw axes are lost (space group number 2) and the surroundings of the two molecules per unit cell are no longer symmetry-equivalent. 1A also stands out somewhat within the first group of systems due to the significantly increased angle between the short molecular axes of the two molecules per unit cell, and due to particularly large normal distances between the planes of equivalent molecules in adjacent layers (for more details see Section S3 (ESI†) and the animations of the phonon modes). All these structural differences are insofar relevant, as they affect the intermolecular interactions, whose impact is superimposed on trends in the phonon characteristics that follow from purely molecular properties, such as molecular lengths, masses, and moments of inertia.

Low-frequency phonon band structures

In view of the above discussion, we will focus on the low-frequency region up to ∼7 THz comprising those phonons that are most relevant, e.g., for heat and charge transport. For a more quantitative discussion see ref. 40. Exemplary phonon band structures for monoclinic benzene,⁶² anthracene,⁶⁴ as well as for the triclinic polymorph I of pentacene⁶⁶ are shown in Fig. 1(a–c). The phonon band structures of naphthalene,⁶³ tetracene,⁶⁵ the thermodynamically more stable^69,70 orthorhombic benzene (with four molecules in the unit cell),⁶⁸ and of polymorph II of pentacene⁶⁷ are given in Section S6 (ESI†). Up to ∼4 THz, the majority of the bands represent intermolecular phonons, which are characterized by motions of the molecules as more or less rigid objects. Of the twelve intermolecular bands (six per molecule in the unit cell), there are three acoustic and nine optical branches, where the latter comprise three translational and six rotational modes. These intermolecular bands are typically characterized by comparably large so-called participation ratios (PRs; see the colouring of the bands in Fig. 1(a–c)). The PRs range between 1 (corresponding to a mode in which each of the N atoms in the unit cell moves with the same amplitude) and N⁻¹ (corresponding to a mode in which only one atom in the unit cell moves) and they describe the degree to which all atoms of the unit cell participate in the respective motion. Thus, they are a measure of the degree of delocalization of a mode (a rigorous mathematical definition of the PRs is given in Section S4.5, ESI†).^71–74 As these intermolecular modes involve the motion of rather large masses (i.e., entire molecules), which are bonded by relatively weak van der Waals forces, it is not surprising that they dominate the phonon properties at low frequencies. Nevertheless, for longer acenes, one still observes an increasing number of bands in the displayed spectral range despite the constant number of intermolecular phonon bands. This is a result of the increasing number of intramolecular modes that are shifted into the low-frequency region with increasing molecular length due to the associated higher molecular flexibility. This finding is in agreement with experimental observations for anthracene in the literature.³⁴

Fig. 1 PBE/D3–BJ-calculated phonon band structures coloured according to the mode participation ratios of (a) the monoclinic polymorph of benzene, (b) anthracene, and (c) polymorph I of pentacene. For a description of the corresponding high-symmetry paths in the first Brillouin zones and the mathematical definition of the participation ratio see Sections S6 and S4.5 (ESI†). Panels (d) and (e) show the evolution of intramolecular (highlighted in (d)) and intermolecular (highlighted in (e)) frequencies at the Γ-point as a function of the molecular length. The frequencies are displayed as short horizontal bars, with the connecting lines serving as a guide to the eye. The horizontal dashed line in (d) marks the frequency interval shown in (e). The types of translational rigid intermolecular modes (TRIMMs) and rotational rigid intermolecular modes (RRIMs) are labelled with the lattice vectors a_i and the molecular directions (long axis, short axis and the axis normal to the π-plane) along/around which the vibrations primarily occur. Animations of the displayed modes can be found in the ESI.† The shaded areas act as a guide to the eye and emphasise the observed trends for the selected intermolecular modes.

To be able to analyse the situation more systematically, Fig. 1(d) and (e) show the evolutions of the Γ-point frequencies of all low-frequency modes. Panel (d), which is plotted over a wider frequency range, highlights the different intramolecular modes, while panel (e) focuses on the intermolecular ones.

Intramolecular modes

As far as the intramolecular modes are concerned, in the low-frequency region, one is not dealing with stretching or bending motions of specific bonds, but rather with distortions of the entire molecular backbone. This is again a consequence of the larger associated masses and smaller associated force constants. As quantified below and shown in Fig. 1(d), the energies of these modes drop significantly for longer acenes. Consequently, for benzene there is still a pronounced band gap between intra- and intermolecular modes, which amounts to ∼8.1 THz. However, already for naphthalene the intramolecular bands approach the intermolecular region and for longer acenes the intra- and intermolecular bands overlap. To be more specific, the intramolecular modes of up to ∼14 THz comprise out-of-plane and in-plane backbone bending as well as backbone torsion vibrations. As there are two molecules per unit cell, each intramolecular vibrational mode results in two bands, with one corresponding to an in-phase (IP) and the other to an antiphase (AP) motion of the molecules. The ensuing displacements for the case of naphthalene are illustrated in Fig. 2(a) (with animations in the ESI†). In passing we note that for the in-plane bending and backbone torsion vibrations, the in-phase and antiphase vibrations are very close in energy (i.e., their difference is sometimes rather poorly resolved in the plot in Fig. 1(d)).

Fig. 2 Displacement patters of the lowest intramolecular (a) bending, (b) torsional, and (c) in-plane bending modes for the two molecules in the unit cell of naphthalene. The equilibrium positions of the molecules (C grey, H white) are shown as grey and white sticks and semi-transparent balls, while the displaced positions are shown as coloured sticks (red and blue for the positive and negative amplitudes, respectively). The arrows indicate the most important motions. The effective bending length, l_eff, is indicated in (a) as the distance between the two nodes of the bending amplitude. The nodal plane of the torsional mode in (b) is indicated by the solid black line. (d) Decrease in the frequency of the (in-plane) bending and torsional modes as a function of the size of the molecule for the five studied acenes showing the in-phase (IP) and antiphase (AP) frequencies. The models (corresponding to the dashed lines) are explained in more detail in Section S7 (ESI†). (e–g) Displacement patterns of the translational rigid intermolecular modes (TRIMMs) of naphthalene (e) along the lattice vector a₃, (f) along a₂, and (g) along a₁. The coloured semi-transparent regions emphasise the parts in space where the most pronounced changes of the (geometric) molecular overlap are observed. (h) Calculated effective intermolecular mode force constant values for the three TRIMMS in the studied molecular crystals as a function of the molecular mass. The dashed lines correspond to fitted linear functions through the origin for the TRIMMS along a₃ (considering all systems for that fit) and a₂ (omitting benzene for the fit for reasons explained in the main text). For the TRIMM along a₁, a constant was fitted.

To explain the shift of the bending modes with acene lengths, we resort to the model of a classical bending beam, whose eigenfrequency is inversely proportional to the square of its length (which, as argued in Section S7.1 (ESI†), is associated with the distance between the nodes of its fundamental vibration, l_eff). Taking naphthalene as a reference, this yields the following relation for the frequency ratios of the corresponding modes in acenes:

The evolution of ω_n obtained from eqn (1) perfectly matches the trend of the band positions in the actual band structures, as shown in Fig. 2(d). An analogous situation also applies to the in-plane bending modes, although with a slightly altered effective length and a higher stiffness (as motivated in Section S7.1, ESI†). Interestingly, the model from eqn (1) primarily follows the evolution of the in-phase modes. This is most pronounced for the backbone bending case, where an increase in the difference between the in-phase and antiphase energies with molecular length is clearly resolved. The differences can be understood as a consequence of the significantly increased intermolecular interactions for antiphase vibrations, which increases their apparent stiffness and, thus, the associated frequencies. This effect becomes particularly relevant for the larger displacement amplitudes of the central regions of the molecules.

A similar classical model can also be used to explain the decrease of torsional frequencies with n, as shown also in Fig. 2(d). The crucial factors here are (i) the torsional moment of inertia and (ii) the torsional stiffness. Both quantities not only depend on the molecular lengths, but also on whether the nodal plane of the torsional vibration at the centre of the molecule cuts through a phenylene ring or lies at the interface between two rings (as detailed in Section S7.2, ESI†). The latter gives rise to the odd–even effects observed in the evolution of the model data and the actual Γ-point frequencies, which are visible in Fig. 2(d).

Note that all the intramolecular modes discussed above show significant shifts to higher frequencies in the crystalline systems compared with the equivalent eigenmodes in the isolated molecules. This is the consequence of the increase of the associated force constants due to the intermolecular interactions in the crystals (a more detailed discussion of these shifts can be found in Section S11, ESI†).

Intermolecular modes

Concerning the modes characterized primarily by motions of rigid molecules, several different cases need to be distinguished: the in-phase translational oscillations of the two molecules in the unit cell form the acoustic bands, which in all acenes display an almost perfectly linear dispersion over a comparably large frequency range (see Fig. 1(a–c)). Moreover, their band widths are relatively similar for different directions in reciprocal space. A more in-depth discussion of band dispersions and differences between transverse and longitudinal modes will be provided below, when analysing the corresponding group velocities. At this stage we would, however, like to comment on the widths of the acoustic bands. Strictly speaking, they decrease with the lengths of the acenes, but with the aid of the participation ratios (i.e., by following the nearly linear band fragments with particularly high PRs) one can conclude that this is primarily a consequence of avoided band crossings, which occur as a consequence of phonon hybridisations involving other intermolecular bands belonging to the same irreducible symmetry representation. This effect is particularly pronounced for the triclinic crystals of tetracene and pentacene. Avoided crossings cause the opening of small band gaps, whose widths depend on the hybridisation strength. They are unfavourable for heat transport because of the (partial) loss of more strongly dispersing bands.^75–77

Antiphase translational modes with the two molecules moving in opposite directions form optical bands, which are also characterized by particularly high degrees of delocalisation (PR ≈ 1). Three such bands exist in the systems considered here, and the evolution of their Γ-point frequencies is shown by the red symbols in Fig. 1(d). These translational rigid intermolecular modes (which we will refer to as TRIMMs) have a pronounced impact on the materials' charge-transport properties, as they comprise motions that result in massive variations of intermolecular electronic couplings.¹⁹ Thus, they are a major source of dynamic disorder and have been associated with killer phonon modes.²⁹ Consequently, pushing them to higher frequencies to reduce the thermal occupation of these modes is expected to improve the carrier mobilities.^22,24–29 Extending the conjugated backbones apparently does not have this effect, as illustrated in Fig. 1(e).

Interestingly, the evolution of each of the three TRIMMs with chain length appears to be fundamentally different with either decreasing or constant energy values. To understand this, one has to analyse the corresponding real space molecular displacements, which are shown in Fig. 2(e–g) (for animations see ESI†). It turns out that in all systems, the lowest frequency TRIMM corresponds to a translational vibration of the molecules along the a₃ axis, which is largely parallel to the long molecular axes (a₃/long; see Fig. 2(e)). For naphthalene and longer acenes, this translation is combined with a minor bending of the molecular backbones. The next higher TRIMM is dominated by an antiphase translation of the two molecules in the a₂ direction, which is the lattice parameter at the smallest angle to the short molecular axes (a₂/short; see Fig. 2(f)). Typically, these vibrations involve hardly any distortion of the molecules; only in tetracene a particularly strong hybridisation with the in-phase bending mode occurs due to the near degeneracy of that mode and the TRIMM. The highest-energy mode corresponds to an antiphase translation in the a₁ direction, which primarily causes a periodic variation of the distances normal to the molecular π-planes (a₁/normal; see Fig. 2(g)). Here, for anthracene and tetracene, significant hybridisations occur with the very-close-lying in-phase and antiphase molecular bending vibrations.

To understand the evolution of the energies of the modes with acene length, it is again useful to consider a simple mechanical model: in a spring–mass system, the (angular) frequency, ω, depends on the mass, M, and the effective spring constant, k_eff, as

ω² = k_effM⁻¹ (2)

As the masses of the different acenes are known, eqn (2) can be used to determine the effective force constants for the different TRIMMs. These are plotted in Fig. 2(h) as a function of the associated molecular mass values. Notably, the effective intermolecular force constant for the highest TRIMM (a₁/normal) scales linearly with the length and molecular mass (see the linear fit through the origin in Fig. 2(h)). This is plausible, as for a vibration in which primarily the distance between neighbouring acenes is varied, the forces (determined by van der Waals attraction and Pauli repulsion) should scale with the intermolecular overlap. Thus, longer molecules show stronger interactions and larger intermolecular force constants for this type of motion. According to eqn (2), a simultaneous linear increase in molecular mass and intermolecular force constant results in a constant vibrational frequency, as observed in Fig. 1(e) for the a₁/normal TRIMM.

The a₂/short TRIMM corresponds to a vibration that primarily induces a change in molecular overlap along the entire length of the molecule, as shown in Fig. 2(f). This suggests that also for the a₂/short TRIMM the effective force constant should increase essentially linearly with the acene length and, thus, with the molecular mass, and the associated frequencies should again stay largely constant. This is indeed observed in Fig. 2(f) and (e) for all acenes apart from benzene for which the force constant is higher than a linear fit through the origin would suggest. This then causes the distinctly higher frequency of the respective mode for benzene shown in Fig. 1(e). We attribute this to peculiarities of the benzene structure, especially the significant increase of the angle between the short molecular axes of the two molecules per unit cell, which results in a steeper increase of intermolecular interactions upon alternatingly displacing the rows of molecules along a₂ (see Section S3, including the provided animations, ESI†).

Finally, for the lowest-energy TRIMM the associated displacement along the long molecular axis should induce at most a weak length dependence of changes in the geometric overlap between neighbouring molecules, as here primarily the overlap of the outermost rings is modified (see Fig. 2(e)). This suggests a length-independent effective force constant, which is indeed observed with the exception of tetracene, for which the force constant is somewhat smaller than expected. Consequently, the vibrational frequency is expected to drop with the one over the square root of the molecular mass, fully consistent with Fig. 1(e), with tetracene again being an outlier here with a too small frequency. The reason for that is not fully understood, but it should be mentioned that tetracene is the first system in the series for which a fundamental change in the tilting direction of the two molecules is observed (see the above description of the geometry and Section S3, ESI†). In addition, tetracene lies at the transition point between systems in which the a₃/long TRIMM couples with an antiphase (for the shorter acenes) vs. an in-phase bending motion (for pentacene).

The dispersions of the TRIMM bands are rather similar for benzene, naphthalene, anthracene and pentacene-I with band widths of around 0.2–1.2 THz (the largest band widths are found for the a₁/normal TRIMMs). A fully quantitative determination of the band widths is, however, complicated by the occurrence of avoided crossings. The rather significant band dispersion for TRIMMs is not unexpected considering that the relative phase of translational motions in neighbouring unit cells (expressed by the wavevector q) strongly impacts the relevant intermolecular distances and, thus, also the intermolecular interactions and the frequencies of the vibrations.

The last class of low-frequency vibrations is the in-phase and antiphase rotations of the molecules, which in the following will be referred to as rotational rigid intermolecular modes (RRIMMs). As shown in Fig. 1(e), the two highest RRIMM bands correspond to rotations around an axis largely parallel to the a₃ direction and, thus, largely around the long molecular axes especially for naphthalene to tetracene (see animations in the ESI†). For benzene, the rotational axis is strongly inclined and in the two pentacene polymorphs there is a quite significant hybridisation with very close-lying intramolecular modes (see Fig. 1(d, e) and Section S9.2, ESI†). For a₃/long-axis rotational modes, one would expect the moment of inertia (as the quantity analogous to the molecular mass in eqn (2) for a rotational mode) to increase linearly with molecular length and a similar evolution would be expected for the restoring torque. Indeed, when plotting the effective rotational force constant (again in analogy to eqn (2)) as a function of the moment of inertia around the long molecular axis (see Section S7.3, ESI†), one observes a linear dependence similar to the a₃/long TRIMM. This is consistent with the rather system-independent frequencies of the in-phase and antiphase a₃/long RRIMMs. Benzene is again an outlier, especially for the in-phase mode, which we primarily attribute to the strong inclination of the rotational axis relative to the long molecular axis mentioned above. The variations in frequency are, however, more pronounced than for the a₃/long TRIMM. This is insofar not surprising, as molecules tilt differently in the different systems: for 1A, 2A, and 3A the tilting directions of the long molecular axes of the two molecules are different with rather significantly varying tilt angles; conversely, both molecules tilt in the same direction for 4A and 5A (see above and Section S3, ESI†). This then causes some variation of the orientation of the rotational axis between the different systems.

The situation becomes even more involved for rotations, which are nominally around the short molecular axes or around axes normal to the molecular plane (i.e., around axes nominally parallel to a₂ and a₁). There, one faces two major problems: first, it is difficult to formulate expectations based on classical models, as the moments of inertia and the restoring forces for such vibrations are expected to have a strongly non-linear dependence on the molecular size. The latter is primarily due to the rather different displacement amplitudes for the outer rings in the different systems, which then cause large differences in the intermolecular interactions. Second, an in-depth inspection of the associated atomic displacements shows that the orientations of the rotational axes for these RRIMMs vary strongly between the systems. Therefore, we provide a detailed discussion of the a₁/normal and a₂/short RRIMMs only in Section S7.3 (ESI†).

As far as the observed band dispersions are concerned, the bands associated with the RRIMMs are not necessarily flat, despite the considerably reduced participation ratios. This can again be attributed to the significant impact of the vibrational displacements on the intermolecular interactions triggering a strong dependence of the frequencies on the phase shift (and, thus, on q). Still, in view of the significantly changing structures of the bands for the different acenes it appears hardly possible to extract general trends. Therefore, we will focus the subsequent discussion of group velocities as a measure of band dispersions on the particularly relevant acoustic phonons and only at the end of the next section will return to higher-lying, optical bands.

Group velocities

Besides their fundamental relevance, group velocities are also of distinct practical importance, as they significantly impact all energy-transport related phonon properties (such as the thermal conductivity). Their full analysis is, however, quite involved, as the group velocity, v_g, is a three-component vector for each phonon. Thus, in the following we will resort to discussing (vector) norms of the group velocity, v_g, and, as a first step, will focus on acoustic phonons in the long-wavelength limit (i.e., at small wave vectors, q). These are representative of the situation for a large fraction of reciprocal space, as the acoustic bands display a linear dispersion (and, thus, a constant group velocity) over a rather wide frequency range (see Fig. 1(a–c)).

Of particular interest are the longitudinal acoustic (LA) bands, which typically display the largest group velocities in a solid and which also correspond to the (longitudinal) sound velocities. At this point it is worth mentioning that purely longitudinal or transverse phonon polarisations with entire lattice planes moving in phase are observed only in high-symmetry directions of highly symmetric crystals, while, in general, the polarisations have a mixed longitudinal–transverse character.^78,79 Nevertheless, in the following, we will refer to longitudinal acoustic and transverse acoustic (TA) modes, whenever they have a primarily longitudinal or transverse character. Prior to discussing the absolute magnitudes of the LA group velocities and their evolution with molecular length, it is interesting to analyse their directional dependence. Fig. 3(a) shows the directional dependence of the deviation of v_g,LA from its mean value of 4.30 km s⁻¹ for 5A (equivalent plots for the other acenes are shown in Section S8, ESI†). This schematic testifies to a pronounced anisotropy of v_g,LA, which is not surprising, considering the highly anisotropic crystal structures of the acenes. As quantified in Fig. 3(b), v_g,LA in 5A varies between 2.77 and 5.33 km s⁻¹, i.e., by nearly a factor of 2. It is largest in a direction essentially parallel to the long molecular axis (and the a₃ axis). This can again be explained by a classical model. The group velocity of acoustic phonons in high-symmetry directions in the long wavelength limit is approximately given by^78–80

with k representing the force constant describing the interactions between the oscillating objects, and M and L representing their mass and distance, respectively. The latter is particularly large along the a₃ axis with the effect being most pronounced for 5A. The distance dependence can be rationalised by the fact that L represents the period of the weak intermolecular van der Waals bonds that are compressed and expanded during wave propagation and, thus, delay the speed of elastic energy migration. Moreover, an in-depth analysis of the displacements shows that at higher q-values (i.e., for large phase shifts between neighbouring unit cells) the phonons also comprise some compression and expansion of the particularly stiff covalent C–C bonds within the molecules. Group velocities in directions parallel to the short molecular axis and perpendicular to the π-planes are intermediate (see the red and green arrows in Fig. 3(a)), while particularly low group velocities are found in directions parallel to a₁ and a₂. This also applies to directions parallel to their projections into the plane perpendicular to a₃ (termed a₁′ and a₂′), which is illustrated in Fig. 3(c). The intermolecular distances between the molecules and their periodic replica in these directions (as a measure of L) are comparably small and, considering the arrangement of the molecules in the plane perpendicular to a₃ (see Fig. 3(c)), it is, indeed, plausible that a longitudinal compression of the unit cell along a₁′ and a₂′ would be energetically not very costly (i.e., the associated force constant would be comparatively small). This is at least suggested by the animations of the displacements for wave propagation in the a₁′- and a₂′-directions (see the animation in the ESI†). In the a₂-direction, they primarily involve π-planes gliding past each other, and in the a₁-direction they are determined by motions of comparably distant and only partially overlapping π-planes towards and away from each other.

Fig. 3 (a and d) Spatial distributions of the variations of the sound velocities, v_g, projected onto unit spheres for the (a) (primarily) longitudinal acoustic (LA) and (d) lower (primarily) transverse acoustic (TA₁) modes in pentacene relative to their arithmetic mean values, v̄_g (4.30 km s⁻¹ and 2.47 km s⁻¹). To facilitate the orientation in space, the directions of the lattice vectors (a₁, a₂, and a₃) are indicated by the black arrows. In addition, the long molecular axes (blue arrows), the short molecular axes (red arrows) and the normal vectors to the π-planes (green arrows) of the two molecules per unit cell are indicated. The three arrows with a slightly darker shade belong to one molecule, while the lighter arrows belong to the other molecule. (b) Average sound velocities of the LA and TA_1/2 phonon branches (bars) for the studied acenes from benzene (1A) to pentacene (5A). The “error margins” indicate the minimum and maximum observed sound velocity values of the respective type of phonon band on the sampled sphere. The distances between the upper edges of the bars and the horizontal coloured dashed lines correspond to the standard deviations of the distributions of sound velocities of the respective type. (c) Polar plot of the LA sound velocity (units: km s⁻¹) in pentacene in a plane perpendicular to the lattice vector a₃ overlaid with the crystal structure seen along this direction. The projections of the remaining two lattice vectors in this plane are labelled with a₁′ and a₂′. (e and f) Two-dimensional densities of states per group velocity and frequency for pentacene. The white dashed box in (f) indicates the region of the expanded area shown in (e).

The above-described anisotropies prevail for the shorter acenes, albeit with a decreasing magnitude of the directional variation of v_g,LA (see Fig. 3(b) and the polar plots in Section S8, ESI†). The only exception is benzene, for which the maximum of v_g,LA is no longer parallel to the molecular axis, which is defined as the “long axis” in accordance with the longer acenes. This is a consequence of the molecule not having its largest extent in that direction (see the discussion in Section S2, ESI†).

Based on the above arguments, one can also understand the variation of v_g with molecular length (see the “error margins” of the green bars in Fig. 3(b)). The minimum group velocities (the lower margin) remain essentially constant, as the situation does not change very significantly as a function of the molecular length for wave propagation close to the a₁,a₂-plane. Here, the force constants and the molecular masses both increase linearly with molecular length, analogous to the above explanation for the system-independent Γ-point energy of the a₁/normal and a₂/short TRIMMs. These modes involve conceptually similar displacements as the LA modes in the a₁- and a₂-directions (albeit of the two molecules within the unit cell relative to each other rather than of molecules in neighbouring cells). Conversely, the increase of the maximum (and concomitantly average) values of v_g can be explained by an increasing intermolecular distance in a₃ direction, which enters linearly into eqn (3), and which is only partly compensated by the increasing molecular mass entering only to the power of −1/2.

As expected, the group velocities associated with the transverse acoustic (TA) modes are significantly smaller than for the longitudinal modes. Interestingly, one typically observes similar minimum, maximum and average values for 1A, 2A, and 3A, while for longer acenes, the minimum values of v_g,TA1 and v_g,TA2 drop and the maximum values increase, such that the anisotropy becomes particularly large. For example, for the TA₁ mode in 5A the ratio between the minimum and maximum values of v_g,TA1 becomes as large as ∼7.1. Moreover, the directions of minimum and maximum group velocities appear in fundamentally different directions than for the longitudinal modes. This is exemplarily shown for the TA₁ mode of 5A in Fig. 3(d), for which the group velocities are particularly small for wave propagation along the short and long molecular axes as well as perpendicular to the π-planes (the corresponding plots for the TA₂ mode and for the other acenes can be found in Section S8, ESI†).

The above considerations exclusively concern the acoustic phonons in the long wavelength limit. However, already the bands shown in Fig. 1(a–c) indicate that in acenes there are also several significantly dispersing optical phonon bands. This applies particularly to the low-frequency region, as shown in Fig. 3(e) for 5A. Interestingly, Fig. 3(f) shows that there are also rather significantly dispersing bands at even higher frequencies (up to ∼35 THz). Examples for high-lying optical bands with maximum v_g values above 1 km s⁻¹ in pentacene comprise intramolecular bending modes around the molecules’ long axes (at ∼13.5 THz and ∼13.7 THz at Γ), out-of-plane bending of the hydrogen atoms in an alternating up-down fashion (at ∼21.7 THz at Γ), and in-plane bending of the C–H bonds (at ∼34.5 THz at Γ) with the six H atoms in the centre of the molecule bending in one direction, while the hydrogens at both ends bend in the other. For all these modes one can understand that intermolecular interactions play a non-negligible role, such that the phase shifts of the displacements in neighbouring unit cells have an impact on the interaction energy such that the phonon frequency becomes strongly q-dependent, which then causes a comparatively large group velocity. Similar considerations apply to the other acenes as shown in Section S4.4 (ESI†), which also contains a brief discussion of the densities of states per group velocity.

Phonon contribution to the heat capacity

While the above discussion focused on low-frequency phonons, there are also quantities for which higher-frequency phonons play a role. This is particularly true for thermodynamic properties at room temperature or at elevated temperatures, at which higher-lying phonon states also become occupied. An important example of such a thermodynamic property is the phonon contribution to the molar heat capacity at constant volume, c_V^m. This quantity (calculated from the full phonon band structures) is shown in Fig. 4(a) as a function of temperature for the studied OSCs together with the experimental values for the molar heat capacities at constant pressure, c_p^m, at room temperature. The simulations fully confirm the experimental trend with the absolute values being somewhat smaller in the calculations than in the measurements. This is not unexpected, as c_p^m is generally larger than c_V^m, since in c_V^m the effect of thermal expansion is not considered.⁸⁰

Fig. 4 (a) Molar heat capacity and (b) heat capacity per unit cell normalised by the classical Dulong–Petit limit (3Nk_B) as a function of temperature for the studied organic semiconductor crystals. The vertical dashed lines indicate the temperature of 300 K, while the horizontal line emphasises the common value in the normalised heat capacities of ∼0.33 observed at that temperature. Note that the upper limit of the shown temperature range already exceeds the melting point of benzene (at ∼279 K) and naphthalene (∼353 K) at ambient pressure.⁸¹ The open circles in (a) indicate the experimentally determined values for the molar constant-pressure heat capacity for benzene,^82,83 naphthalene,^83,84 anthracene,^83,85 tetracene,⁸⁶ and pentacene⁸⁶ in the corresponding colours. (c) Expanded region of panel (b) indicated by the dashed blue rectangle. (d) Density of states (DOS) normalised with the number of bands (3N) of the studied acenes: benzene (1A) to pentacene (5A).

For comparing the different acenes on an “equal footing” it is useful to convert the molar heat capacity, c_V^m, to the heat capacity per unit cell, C_V (by dividing c_V^m by Avogadro's number and multiplying it by the number of molecules per unit cell). In the classical Dulong–Petit limit (i.e., for the temperature T → ∞) C_V amounts to

Here, N is the number of atoms contained in the unit cell, and k_B is the Boltzmann constant. Notably, the normalised heat capacity per unit cell, C_V/(3Nk_B), ranges between 0 and 1 regardless of the number of atoms or molecules per unit cell. It is plotted in Fig. 4(b) as a function of temperature for the different acenes with an expanded version of the low-temperature region in Fig. 4(c).

Starting from C_V/(3Nk_B) = 0 at 0 K (in agreement with the third law of thermodynamics), the heat capacity rises continuously with T, concomitant with an increasing occupation of the higher-frequency states.^78,80 At room temperature, the heat capacity is still far from its classical limit (C_V/(3Nk_B) = 1), which is a consequence of the many high-frequency modes occurring in acenes (with the C–H stretching vibrations at ∼93 THz as the most extreme example). Interestingly, at room temperature, C_V/(3Nk_B) is virtually identical for all the systems and amounts to ∼0.33. This implies that the nature of the atoms contained in the molecules and especially the carbon-to-hydrogen ratio (which is 1.0 for benzene and 1.6 for pentacene) play essentially no role for that quantity. Consequently, all that really counts for the molar heat capacity is the total number of atoms per molecule (see Section S10 for more details, ESI†), such that the increase of the molar heat capacity with the length of the acenes at room temperature (see Fig. 4(a)) is merely a consequence of the larger numbers of atoms in the bigger molecules. The explanation for this lies in the (normalised) phonon densities of states shown in Fig. 4(d). They show that with increasing complexity of the molecules, the additional intramolecular vibrations are relatively uniformly distributed over the entire frequency range. Variations in the normalised heat capacity would occur only if the density of states in a specific frequency region differed notably for a given material. This is, in fact, the case in the low-frequency region, where there are twelve intermolecular modes per unit cell, independent of the length of the molecules (see above). Consequently, for the shorter acenes the weight of the intermolecular modes in that region increases relative to the entire DOS. For the longer acenes, this is only partially compensated by the decrease of the frequencies of the first intramolecular modes with molecular length as discussed above. As a result, at low temperatures (where mostly the low-frequency phonons count) the normalized heat capacity per unit cell is significantly increased, especially in benzene and also in naphthalene (see Fig. 4(c)).

Conclusions

This article uses state-of-the art density-functional-theory-based simulations to explain the relationship between the phonon band structures of acenes and the respective structures of the molecules and crystals. The discussion focuses on the low-frequency region (up to ∼7 THz), which is particularly relevant for heat and charge-transport processes. It is dominated by intermolecular vibrations (i.e., motions of the two molecules per unit cell as more or less rigid objects) and by intramolecular modes characterised by distortions of the entire molecular backbones. At first glance, the calculated phonon band structures appear to be highly complex with no indications for apparent structure–property relationships. An in-depth analysis of the individual bands, however, allows us to distinguish between phonons associated with intramolecular in-phase and out-of-phase bending and torsion modes, as well as intermolecular “translational and rotational rigid intermolecular modes” (TRIMMs and RRIMMs). The evolution of the intramolecular modes and of the TRIMMS with the acene length can be rationalized very well using simple classical models of bending beams and (torsional) spring–mass oscillators, and the details of the arrangements of the molecules in the different acene crystals typically play only a minor role. The simple models allow us to rationalise, why the frequencies of the considered intramolecular modes drop with acene length and why such a behaviour is observed only for certain TRIMMs. The situation becomes less clear for the RRIMMS, especially for those associated with rotations around axes other than the long molecular axis. This is primarily attributed to the system-dependent orientations of the respective rotational axes due to variations in the arrangements of the acene molecules in the unit cells.

As far as the group velocities of the acoustic phonons are concerned, we observe pronounced anisotropies with the highest group velocities reaching values of up to 5.33 km s⁻¹. For the longitudinal modes these anisotropies can again be rationalised based on classical models and the crystal structure. This also applies to the evolution of the maximum and minimum longitudinal group velocities. Unfortunately, the situation becomes much more complex for the (primarily) transverse modes. Considering the band dispersions and group velocities of the optical modes above 2 THz, one finds values of up to ∼2.8 km s⁻¹, which is only a factor of ∼1.9 smaller than the fastest acoustic phonons. Notably, even for frequencies up to 35 THz, bands with group velocities beyond 1 km s⁻¹ are observed, which can be rationalised by the nature of the involved vibrations. Finally, as an example of a thermodynamic quantity, we show how the phonon properties correlate with the phonon heat capacity as a function of temperature and chain length.

In conclusion, the presented results show how analogies to classical macroscopic oscillators can be exploited to estimate and explain the shifts in phonon frequencies and the associated change in physical properties in a series of related molecular crystals. These considerations are, by no means, restricted to the here-studied series of acenes but can also be extrapolated to other organic semiconductors. Thus, we expect the above observations and their explanations to serve as a foundation on which the future systematic investigation of phonon properties of more complex organic semiconductor crystals can be built.

Author contributions

T. K. and E. Z. conceptualised the present study of the structure–property relationships in phonon-derived properties of organic semiconductor crystals. T. K. conducted all the quantum–mechanical simulations as well as the formal analysis of the results, including coding of the necessary analysis software tools. T. K. wrote the original draft of the manuscript and prepared all the figures, which were revised by T. K. and E. Z., who was also responsible for the supervision of the project.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was funded in part by the Austrian Science Fund (FWF) [P33903-N]. We also acknowledge the Graz University of Technology for financial support through the Lead Project (LP-03) and the use of HPC resources provided by the ZID. The computational results have been in part achieved using the Vienna Scientific Cluster (VSC3). T. K. acknowledges additional funding as the recipient of a DOC Fellowship (25783) of the Austrian Academy of Sciences at the Institute of Solid State Physics and the Institute of Physical and Theoretical Chemistry.

Notes and references

1. H. Bronstein, C. B. Nielsen, B. C. Schroeder and I. McCulloch, Nat. Rev. Chem., 2020, 4, 66–77.
2. C. Wang, H. Dong, L. Jiang and W. Hu, Chem. Soc. Rev., 2018, 47, 422–500.
3. A. Salehi, X. Fu, D. H. Shin and F. So, Adv. Funct. Mater., 2019, 29, 1–21.
4. G. A. Chamberlain, Sol. Cells, 1983, 8, 47–83.
5. N. Yeh and P. Yeh, Renewable Sustainable Energy Rev., 2013, 21, 421–431.
6. B. Kumar, B. K. Kaushik and Y. S. Negi, Polym. Rev., 2014, 54, 33–111.
7. J. Nyman and G. M. Day, CrystEngComm, 2015, 17, 5154–5165.
8. K. T. Butler, K. Svane, G. Kieslich, A. K. Cheetham and A. Walsh, Phys. Rev. B, 2016, 94, 1–5.
9. E. Venuti, R. G. Della Valle, A. Brillante, M. Masino and A. Girlando, J. Am. Chem. Soc., 2002, 124, 2128–2129.
10. A. Brillante, R. G. D. Valle, L. Farina, A. Girlando, M. Masino and E. Venuti, Chem. Phys. Lett., 2002, 357, 32–36.
11. B. P. van Eijck, J. Comput. Chem., 2001, 22, 816–826.
12. A. J. Cruz-Cabeza, S. M. Reutzel-Edens and J. Bernstein, Chem. Soc. Rev., 2015, 44, 8619–8635.
13. D. Wang, L. Tang, M. Long and Z. Shuai, J. Phys. Chem. C, 2011, 115, 5940–5946.
14. J. C. Duda, P. E. Hopkins, Y. Shen and M. C. Gupta, Appl. Phys. Lett., 2013, 102, 251912.
15. J. C. Duda, P. E. Hopkins, Y. Shen and M. C. Gupta, Phys. Rev. Lett., 2013, 110, 1–5.
16. M. Chabinyc, Nat. Mater., 2014, 13, 119–121.
17. Y. Zhang, Y. J. Heo, M. Park and S. J. Park, Polymers, 2019, 11, 167.
18. O. Bubnova, Z. U. Khan, A. Malti, S. Braun, M. Fahlman, M. Berggren and X. Crispin, Nat. Mater., 2011, 10, 429–433.
19. V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J.-L. Brédas, Chem. Rev., 2007, 107, 926–952.
20. N. Karl, Synth. Met., 2003, 134, 1–9.
21. N. Tessler, Y. Preezant, N. Rappaport and Y. Roichman, Adv. Mater., 2009, 21, 2741–2761.
22. S. Fratini, D. Mayou and S. Ciuchi, Adv. Funct. Mater., 2016, 26, 2292–2315.
23. H. Oberhofer, K. Reuter and J. Blumberger, Chem. Rev., 2017, 117, 10319–10357.
24. A. Troisi and G. Orlandi, J. Phys. Chem. A, 2006, 110, 4065–4070.
25. A. Troisi and G. Orlandi, Phys. Rev. Lett., 2006, 96, 1–4.
26. S. Illig, A. S. Eggeman, A. Troisi, L. Jiang, C. Warwick, M. Nikolka, G. Schweicher, S. G. Yeates, Y. Henri Geerts, J. E. Anthony and H. Sirringhaus, Nat. Commun., 2016, 7, 1–10.
27. S. Fratini, S. Ciuchi, D. Mayou, G. T. De Laissardière and A. Troisi, Nat. Mater., 2017, 16, 998–1002.
28. T. F. Harrelson, V. Dantanarayana, X. Xie, C. Koshnick, D. Nai, R. Fair, S. A. Nuñez, A. K. Thomas, T. L. Murrey, M. A. Hickner, J. K. Grey, J. E. Anthony, E. D. Gomez, A. Troisi, R. Faller and A. J. Moulé, Mater. Horiz., 2019, 6, 182–191.
29. G. Schweicher, G. D’Avino, M. T. Ruggiero, D. J. Harkin, K. Broch, D. Venkateshvaran, G. Liu, A. Richard, C. Ruzie, J. Armstrong, A. R. Kennedy, K. Shankland, K. Takimiya, Y. H. Geerts, J. A. Zeitler, S. Fratini and H. Sirringhaus, Adv. Mater., 2019, 31, 1902407.
30. R. S. Sánchez-Carrera, P. Paramonov, G. M. Day, V. Coropceanu and J. L. Brédas, J. Am. Chem. Soc., 2010, 132, 14437–14446.
31. M. Masino, A. Girlando, A. Brillante, L. Farina, R. G. Delia Valle and E. Venuti, Macromol. Symp., 2004, 212, 375–380.
32. E. L. Bokhenkov, E. M. Rodina, E. F. Sheka and I. Natkaniec, Phys. Status Solidi B, 1978, 85, 331–342.
33. I. Natkaniec, E. L. Bokhenkov, B. Dorner, J. Kalus, G. A. MacKenzie, G. S. Pawley, U. Schmelzer and E. F. Sheka, J. Phys. C-Solid State Phys., 1980, 13, 4265–4283.
34. B. Dorner, E. L. Bokhenkov, S. L. Chaplot, J. Kalus, I. Natkaniec, G. S. Pawley, U. Schmelzer and E. F. Sheka, J. Phys. C-Solid State Phys., 1982, 15, 2353–2365.
35. V. K. Jindal, J. Kalus, E. L. Bokhenkov, S. L. Chaplot, B. Dorner, I. Natkaniec, G. S. Pawley and E. F. Sheka, J. Phys. C-Solid State Phys., 1982, 15, 7283–7294.
36. V. K. Jindal and J. Kalus, Phys. Status Solidi, 1986, 133, 89–99.
37. R. Bhandari and V. K. Jindal, J. Phys.: Condens. Matter, 1991, 3, 899–907.
38. U. Schmelzer, E. L. Bokhenkov, B. Dorner, J. Kalus, G. A. Mackenzie, I. Natkaniec, G. S. Pawley and E. F. Sheka, J. Phys. C-Solid State Phys., 1981, 14, 1025–1041.
39. J. F. J. Jordan, J. Kalus, U. Schmelzer and G. Eckold, Phys. Status Solidi, 1989, 155, 89–105.
40. T. Kamencek, S. Wieser, H. Kojima, N. Bedoya-Martínez, J. P. Dürholt, R. Schmid and E. Zojer, J. Chem. Theory Comput., 2020, 16, 2716–2735.
41. F. Brown-Altvater, T. Rangel and J. B. Neaton, Phys. Rev. B, 2016, 93, 60–63.
42. J. George, R. Wang, U. Englert and R. Dronskowski, J. Chem. Phys., 2017, 147, 074112.
43. N. Bedoya-Martínez, A. Giunchi, T. Salzillo, E. Venuti, R. G. Della Valle and E. Zojer, J. Chem. Theory Comput., 2018, 14, 4380–4390.
44. C. R. Groom, I. J. Bruno, M. P. Lightfoot and S. C. Ward, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 2016, 72, 171–179.
45. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
46. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
47. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
48. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
49. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
50. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
51. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
52. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
53. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
54. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
55. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
56. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
57. P. Vinet, J. R. Smith, J. Ferrante and J. H. Rose, Phys. Rev. B: Condens. Matter Mater. Phys., 1987, 35, 1945–1953.
58. G. P. Francis and M. C. Payne, J. Phys.: Condens. Matter, 1990, 2, 4395–4404.
59. T. Kamencek and E. Zojer, 2021, NOMAD Dataset Phonons Crystalline Acenes, DOI: 10.17172/NOMAD/2021.09.28-1.
60. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
61. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
62. G. J. Piermarini, A. D. Mighell, C. E. Weir and S. Block, Science, 1969, 165, 1250–1255.
63. S. C. Capelli, A. Albinati, S. A. Mason and B. T. M. M. Willis, J. Phys. Chem. A, 2006, 110, 11695–11703.
64. C. P. Brock and J. D. Dunitz, Acta Crystallogr., Sect. B: Struct. Sci., 1990, 46, 795–806.
65. D. Holmes, S. Kumaraswamy, A. J. Matzger and K. P. C. Vollhardt, Chem. – Eur. J., 1999, 5, 3399–3412.
66. R. B. Campbell, J. M. Robertson and J. Trotter, Acta Crystallogr., 1962, 15, 289–290.
67. C. C. Mattheus, A. B. Dros, J. Baas, A. Meetsma, J. L. De Boer and T. T. M. Palstra, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 2001, 57, 939–941.
68. G. E. Bacon, N. A. Curry, S. A. Wilson and R. Spence, Proc. R. Soc. London, Ser. A, 1964, 279, 98–110.
69. R. Guido, D. Valle and E. Venuti, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 206–212.
70. S. Azadi and R. E. Cohen, J. Chem. Phys., 2016, 145, 064501.
71. R. J. Bell, P. Dean and D. C. Hibbins-Butler, J. Phys. C-Solid State Phys., 1970, 3, 2111–2118.
72. J. T. Edwards and D. J. Thouless, J. Phys. C-Solid State Phys., 1972, 5, 807–820.
73. J. Canisius and J. L. van Hemmen, J. Phys. C-Solid State Phys., 1985, 18, 4873–4884.
74. T. Kamencek, N. Bedoya-Martínez and E. Zojer, Phys. Rev. Mater., 2019, 3, 116003.
75. M. Christensen, A. B. Abrahamsen, N. B. Christensen, F. Juranyi, N. H. Andersen, K. Lefmann, J. Andreasson, C. R. H. H. Bahl and B. B. Iversen, Nat. Mater., 2008, 7, 811–815.
76. W. Li, J. Carrete, G. K. H. H. Madsen and N. Mingo, Phys. Rev. B, 2016, 93, 1–5.
77. T. Zhu, K. Swaminathan-Gopalan, K. J. Cruse, K. Stephani and E. Ertekin, Adv. Funct. Mater., 2018, 28, 1–8.
78. C. Kittel, Introduction to Solid State Physics, Wiley, New York, 2004.
79. J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids, Oxford University Press, Oxford, 2001.
80. D. C. Wallace and H. Callen, Am. J. Phys., 2005, 40, 1718–1719.
81. D. R. Lide, CRC Handbook of Chemistry and Physics, Taylor & Francis, Boca Raton, Florida, 2009, 90th edn.
82. R. G. Ross, P. Andersson and G. Bäckström, Mol. Phys., 1979, 38, 377–385.
83. R. G. Ross, P. Andersson and G. Bäckström, Mol. Phys., 1979, 38, 527–533.
84. J. P. McCullough, H. L. Finke, J. F. Messerly, S. S. Todd, T. C. Kincheloe and G. Waddington, J. Phys. Chem., 1957, 61, 1105–1116.
85. M. Radomska and R. Radomski, Thermochim. Acta, 1980, 40, 405–414.
86. M. Fulem, V. Laštovka, M. Straka, K. Růžička and J. M. Shaw, J. Chem. Eng. Data, 2008, 53, 2175–2181.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1tc04708f

---