Loss spectroscopy of molecular solids: combining experiment and theory

The nature of the lowest-energy electronic excitations in prototypical molecular solids is studied here in detail by combining electron energy loss spectroscopy (EELS) experiments and state-of-the-art many-body calculations based on the Bethe–Salpeter equation. From a detailed comparison of the spectra in picene, coronene and tetracene we generally find a good agreement between theory and experiment, with an upshift of the main features of the calculated spectrum of 0.1–0.2 eV, which can be considered the error bar of the calculation. We focus on the anisotropy of the spectra, which illustrates the complexity of this class of materials, showing a high sensitivity with respect to the three-dimensional packing of the molecular units in the crystal. The differences between the measured and the calculated spectra are explained in terms of the small differences between the crystal structures of the measured samples and the structural model used in the calculations. Finally, we discuss the role played by the different electron–hole interactions in the spectra. We thus demonstrate that the combination of highly accurate experimental EELS and theoretical analysis is a powerful tool to elucidate and understand the electronic properties of molecular solids.

between theory and experiment, with an upshift of the main features of the calculated spectrum of 0.1-0.2 eV, which can be considered the error bar of the calculation. We focus on the anisotropy of the spectra, which illustrates the complexity of this class of materials, showing a high sensitivity with respect to the three-dimensional packing of the molecular units in the crystal. The differences between the measured and the calculated spectra are explained in terms of the small differences between the crystal structures of the measured samples and the structural model used in the calculations. Finally, we discuss the role played by the different electron-hole interactions in the spectra. We thus demonstrate that the combination of highly accurate experimental EELS and theoretical analysis is a powerful tool to elucidate and understand the electronic properties of molecular solids.

Introduction
Organic molecular solids-especially those that consist of fused aromatic rings only and do not contain heteroatoms or carry substituents-have become a focus of research in recent decades and have established a very attractive field at the crossover of solid-state physics and chemistry [1]. The wide variety of these compounds, the advantages of the relatively low cost and the possibility of modifying them using the methods of synthetic organic chemistry in a practically unlimited fashion have aroused high expectations for the development of new materials. In particular, their potential application in organic electronic devices has motivated many investigations in the past. These have for instance been exploited in organic field effect transistors in view of fundamental as well as applied aspects [2]. Moreover, organic semiconductors are also of interest for the manufacture of organic photovoltaic cells, organic light emitting diodes or organic spintronics [3][4][5].
Organic molecular materials, like other carbon-based materials, are well-known prototypes of electrical insulators and semiconductors. Even before the discovery of graphene many other carbon-based systems, such as nanotubes, fullerenes and graphite have played an important role in several fields, including superconductivity. In particular, the superconducting (alkali-metal doped) fullerides have attracted much attention, and rather high transition temperatures could be realized (e.g. T c of 18 K in K 3 C 60 [6] or T c = 38 K in Cs 3 C 60 [7,8]). Their exceptional electronic properties are attributed to the delocalized π -electrons and to their molecular structure (e.g. to dynamical Jahn-Teller effects [9]). Thus, doping of such materials with π-electron networks might bring out novel physical properties, such as superconductivity, as well in other compounds.
However, in the case of organic superconductors, no new system with high T c 's similar to those of the fullerides has been discovered in the past decade. In 2010 the field was finally renewed with the discovery of superconductivity in alkali-doped picene with a T c up to 18 K in case of K 3 picene [10]. After that, superconductivity was also reported in other alkalimetal intercalated polycyclic aromatic hydrocarbons, such as phenanthrene (T c = 5 K) [11], coronene (T c = 15 K) [12] and 1,2;8,9-dibenzopentacene (T c = 33 K) [13]. Especially in the latter case, T c is higher than in any other organic superconductor besides the alkali-metal doped fullerides. Therefore, these small-molecule-based materials open up new avenues for research into superconductivity and it becomes more and more evident that organic materials serve as a fascinating field for material science and solid-state physics.
Understanding the physical properties of new materials requires the investigation of their elementary excitations. The scope of the present work is to unveil the character and the nature of those excitations. To this end, electron energy loss spectroscopy (EELS) is a powerful tool to access the electronic excitations of materials. It measures the loss function, −Im[1/ M (q, ω)], where M is the macroscopic dielectric function, and q and ω are the transferred momentum and energy, respectively. While in optics only the dipole limit q → 0 of M is accessible, in EELS a range of momentum transfers q well beyond the optical limit can be investigated, providing detailed information on the character and the spatial localization of the electronic excitations [14][15][16].
Organic molecular solids are complex materials in which the properties of the crystal as a whole show a strong sensitivity with respect both to the individual properties of their building blocks and to their geometrical arrangement in the three-dimensional packing. This results, for instance, in a high level of anisotropy of the EELS and optical spectra. Moreover, molecular solids often can have many polymorphs obtained by changing the preparation conditions and undergo several structural phase transitions as a function of temperature or pressure. In view of this inherent complexity, a combined experimental and theoretical effort hence seems to be a very promising strategy to provide deeper insight into the fundamental low-energy electronic excitations of these materials. Here we present a detailed comparison of EELS spectra of several prototypical organic solids obtained using state-of-the-art experimental tools and first-principles many-body calculations. We focus on the complementarity aspects of theory and experiment and show how their combination can add a fruitful piece of clarification about the electronic properties of this kind of materials.

Molecular solids
The building blocks of all the compounds studied here are hydrocarbon molecules made of fused benzene rings. While a zigzag or armchair manner (see figure 1). Coronene, instead, is made of benzene rings joined in a circle. In the condensed phase all the investigated hydrocarbon molecules crystallize in a layered structure, often called herringbone crystal structure [17], with two molecules per unit cell. The crystal structure is either monoclinic (as for picene and coronene) or triclinic (as for tetracene).
In the solid, the molecular units are interacting through weak van der Waals forces. This implies that the band widths are generally small and the materials are insulators. Moreover the electronic properties of the solid are mainly dictated by the properties of the isolated molecules: the lowest-energy excitations are due to transitions between bands which originate from π and π * molecular states. This is the reason why excitations in these systems traditionally have been described using models based on molecular orbitals [1,18]. Molecular solids thus offer a unique playground to understand the nature of localized excitons and their interactions in the materials thanks to simple models.

Sample preparation
Single crystals of several hydrocarbons (tetracene, picene and coronene) with high quality were obtained via directional sublimation or physical vapour transport (details of the used techniques can be found in [21,22]). Thin films were cut from these crystals using an ultramicrotome equipped with a diamond knife.
For the present work, large thin films of organic compounds have also been produced by thermal evaporation under high vacuum onto single crystalline substrates (e.g. KBr) kept at room temperature (for details of this procedure and an overview over the parameters used see [23][24][25]). These thin films were floated off in distilled water, mounted onto standard electron microscopy grids and transferred into the spectrometer. In all cases the thickness of the films was about 100 nm.  [19,20]. Prior to the EELS measurements the films were characterized in situ using electron diffraction, whereby on the one hand a check to identify the crystal structure can be performed and on the other hand the orientation of single crystalline samples can be obtained. As an example for such an electron diffraction profile we show in figure 2 the profiles for solid picene and coronene measured in the two main crystallographic directions a * and b * . The measured diffraction profiles are in good correspondence to the published crystal structures of each compound [19,20,26]. In table 1 we summarize the structural information that can be directly read from our diffraction data.
The experimental structures from [19,20,26] have been used for the simulations in the various compounds. All calculations are done considering the atoms frozen in their crystallographic positions, thus neglecting the coupling with atomic vibrations. On the contrary, the effect of temperature-induced structural changes and of atomic motions can be monitored experimentally as a function of the temperature of the measurement.

Electron energy loss spectroscopy
All electron diffraction studies and loss function measurements were carried out using the 172 keV spectrometer described in detail elsewhere [27]. We note that at this high primary beam energy only singlet excitations are possible. The energy and momentum resolution were chosen to be 85 meV and 0.03 Å −1 , respectively. We have measured the loss function −Im[1/ M (q, ω)] for a momentum transfer q parallel to the film surface. In order to properly normalize the data taken along different directions in the reciprocal space, we carried out a Kramers-Kronig analysis (KKA). The raw data first have been corrected by subtracting contributions of multiple scattering processes and elimination of the contribution of the direct beam. The normalization necessary within the KKA has been done using the sum rule for all valence excitations. The results of the KKA show that data for different directions can be properly normalized in the energy range between 9-10 eV to reveal the relative intensities of the corresponding excitations.
The loss function can be written in terms of the real and imaginary parts of the dielectric function, 1 and 2 , respectively, Peaks in the measured spectra can thus derive from zeroes of 1 (when 2 is not too large), which correspond to collective plasmon excitations, or from structures in 2 , which are due to valence-conduction interband transitions. In particular, peaks of 2 that are located at energies smaller than the quasiparticle (QP) band gap correspond to bound excitons. The binding energy of the excitons is then defined as the difference between the band gap and their peak position in the spectra. For each vanishing momentum transfer q, the macroscopic dielectric function M (q, ω) can be calculated by solving the many-body Bethe-Salpeter equation (BSE) [28] M (q, ω) = 1 + lim where A λ and E λ are eigenstates and eigenvalues of the excitonic Hamiltonian written in a one-particle electron-hole (e-h) transition basis |n = |kvc = φ vk (r)φ ck (r) (where v (c) runs over valence (conduction) bands and k is in the first Brillouin zone). In equation (3) E nk are the QP energies,v is a modified Coulomb interaction in which the long-range componentv(G = 0) is set to 0 in reciprocal space, and W is the statically screened Coulomb interaction, calculated in the random phase approximation (RPA). From the eigenvectors A λ of the excitonic Hamiltonian one also obtains the exciton wavefunction as where r h and r e are the positions of the hole and the electron, respectively. Whenv and W are set to 0 in equation (3), the eigenstates A λ become diagonal in the e-h transition basis and equation (2) reduces to Fermi's golden rule in an independentparticle picture. The spectrum is thus given by a sum of independent transitions from valence to conduction states E vk → E ck , which can be calculated either using Kohn-Sham results from the local-density approximation (LDA) of density-functional theory or using many-body perturbation theory in the GW approximation [29] in which the self-energy is approximated by the product of the single particle Green's function G and the screened interaction W.
v and W enter the BSE kernel in equation (3) as e-h exchange and direct Coulomb interactions, respectively. The former is repulsive and is responsible for crystal local-field effects (LFE), the latter is attractive and contains excitonic effects.
Alternatively, the loss function can be calculated using time-dependent density functional theory (TDDFT) in the linear-response framework [30,31], where the polarizability χ can be obtained from a Dyson-like equation where χ 0 is the Kohn-Sham non-interacting polarizability, v the Coulomb interaction and f xc the exchange-correlation TDDFT kernel. When f xc = 0 the RPA is retrieved. In the present work we have also adopted the adiabatic LDA (ALDA) for f xc . From χ one obtains the microscopic dielectric function from and the macroscopic dielectric function M through [32,33] where G and G are reciprocal-lattice vectors and −1 G,G (q, ω) is the Fourier transform to reciprocal space of −1 (r 1 , r 2 , ω).

Calculation details
Following the computational approach adopted analogously in [16,23,34], all the calculations have been done within a plane-wave basis framework using norm-conserving Troullier-Martins pseudopotentials in the LDA with the Perdew-Zunger parameterization. Convergence has been reached with an energy cutoff of 40 Ha in the plane-wave basis set. For picene and tetracene we have used a 6 × 6 × 4 Monkhorst-Pack grid of k points and a 4 × 4 × 4 one for coronene.
GW corrections have been obtained from the standard first-order perturbative G 0 W 0 scheme with a plasmon-pole approximation. For self-energy calculations we have used 7000 plane waves in the expansion of the wavefunctions, 450, 600 and 550 empty bands in picene, coronene and tetracene, respectively. In the calculation of the RPA screened Coulomb interaction W = −1 v, the inverse of the dielectric function −1 was a matrix of size 600, 500, and 600 G vectors for picene, coronene and tetracene, respectively. BSE calculations have been done including the full static W G,G matrix and not only its diagonal, and considering 53 occupied and 29 empty bands for picene, which become 39 and 22 for coronene and 35 and 36 for tetracene. This set of bands gives spectra converged in the energy range shown in the following figures. GW corrections have been taken into account explicitly without using a scissor correction. Spectra have been calculated using the Haydock recursive technique. A broadening of 0.05 eV has been used.

Electron-hole interactions
From the theoretical point of view, it is straightforward to break down a calculated spectrum into the different contributions that give rise to it. This allows one to understand the role played by the different e-h interactions and make the connection with the relevant physical properties. We consider here the representative case of picene in the long wavelength limit q → 0 along the a * axis, but a similar analysis could be carried out for the other materials as well.
In an independent-particle picture, in which no e-h interactions are at play, the calculated spectrum is the sum of vertical e-h transitions. In a first approach the transitions are calculated in LDA (green line in figure 3). In the low-energy range here considered the different peaks in the loss function are due to interband transitions between bands derived from molecular energy levels. The spectra for 2 , which one would obtain from optical absorption spectroscopy, are thus pretty similar.
We then switch on the e-h interactions, starting from the e-h exchange repulsive interactionv. This is responsible of crystal LFE, which are the consequence of the spatial charge fluctuations that are induced on a microscopic scale by the external field. These induced local fields are stronger when the electronic system is polarizable and inhomogeneous. This is the case for picene along a * . When local fields are included in the calculation, the spectra get damped (as the result of a depolarization effect) and shifted to higher energy (compare blue and green lines in figure 3).
While Kohn-Sham LDA eigenvalues are not meant to describe electronic excitation (addition or removal of electrons), the GW approximation of many-body perturbation theory generally gives QP band structure in good agreement with experiment [28]. In the considered molecular solids, the main effect of GW corrections with respect to LDA is the band gap opening, from 2.4 eV to 4.1 eV in picene [23] (in coronene, going from LDA to GW, the band gap opens from 2.4 eV to 4.0 eV; in tetracene from 1.2 eV to 2.9 eV). The effect of the band gap opening associated to the GW corrections is thus a blueshift of the spectra (compare blue and black lines in figure 3).
Finally, we consider the effect of the e-h direct attractive interaction W . The main consequence of W is a redshift of the spectra, which partially compensates the blueshift due to GW corrections (compare black and red lines in figure 3). The final BSE spectrum is pretty similar to the one obtained in LDA (compare red and blue lines in figure 3), where both GW corrections and excitonic effects are absent. Hence the mutual cancellation of these two contributions could allow one to use LDA spectra (with the inclusion of LFE) to make qualitative comparisons with experiments (see yellow curve in the left and right panel of figure 3), in particular for the higher energy peaks at 4.5-6.0 eV. However, the onset in LDA is too low and the interpretation of the lowest energy peaks is different in the BSE calculation. Their energy is smaller than the QP band gap (4.1 eV in picene). Thus they correspond to bound excitons.
In the TDDFT framework, instead, the kernel f xc should contain the two partially compensating effects of the GW band-gap opening and the e-h attraction that is responsible for the formation of bound excitons. However, it has been shown that TDDFT kernels that do not have a long range 1/q 2 contribution in the long-wavelength limit are not able to capture the exciton formation in solids [37]. In figure 3 we in fact see that the ALDA spectrum is very similar to the LDA with LFE included (compare purple and blue lines), with the onset located at much too low energy (2.6 eV). Moreover it is missing the multiple bound-exciton structure inside the band gap (see the right panel of figure 3), for which a dynamical long-range kernel would be necessary [38][39][40][41][42]. The best overall agreement between the calculation and the experiment, especially for those bound excitons, is obtained with the BSE (see right panel of figure 3), with an upshift of the calculated spectrum of ∼ 0.1-0.2 eV, which can be considered as the error bar of the calculation. Therefore the BSE will be the framework used in the following for the analysis of the bound excitons in the spectra.

Anisotropy in the EELS spectra: picene, coronene and tetracene
A crucial aspect of the dielectric properties of all the organic compounds is their strong anisotropy. This is the result of the localization of the wavefunctions, which have a strong molecular character. For instance, in a picene molecule, the lowest energy π → π * transitions are visible only for q directions in the plane of the molecule, while they are dipole forbidden in the direction perpendicular to the plane. In turn, in the condensed phase the resulting spectra are highly sensitive to the arrangement of the molecules with respect to the crystal axes. In picene the a * and b * axes have a large component in the direction perpendicular to the planes of the molecules. Therefore π → π * transitions for momentum transfers q in the a * b * plane have lower intensities than for q along the c * axis. This also implies that a slight change in the orientation of the molecules has a large impact on the intensity of the peaks.
From the experimental point of view, in order to investigate the change in the spectra between the different q directions, it is essential to perform measurements on single crystalline samples. This is illustrated for picene in figure 4 (left panel) and for coronene in figure 5 (left panel). In these two figures we show a comparison of the loss functions measured in  polycrystalline films as well as in single crystals in an energy range of 2-7 eV. These data are taken with a small momentum transfer q of 0.1 Å −1 . In the spectra taken from polycrystalline films the spectral features from different directions are averaged out. In particular here we find that the q-direction dependence of the experimental spectra of the single crystals is larger in picene than in coronene.
In both measured directions one can observe maxima in the loss functions of the two compounds which can be seen as a signature of the energetically sharp and well-defined molecular electronic levels of the two molecules. In fact we can ascribe the spectral maxima to excitations between the energetically close-lying highest occupied and lowest unoccupied electronic states of the two crystals. Interestingly, the main features of the loss spectrum of solid coronene are in good agreement with previous EELS measurements in the gas phase [43,44] and optical absorption data [45,46].
Considering the onset of the experimental spectra, we find an optical gap of 2.8 eV in coronene and of 3.15 eV in picene. The excitation onset is followed by additional well-separated features. Unfortunately, for coronene different values for the transport energy gap from 3.29 up to 3.62 eV were previously published [47,48]. Consequently, from these experimental results only the lowest excitation that is observed can safely be attributed to a singlet exciton.
The main features of the spectra are well reproduced by the calculations, which can therefore safely be used to interpret the experiments. Besides the a * , and b * directions, as in the experiments, BSE calculations are performed also for q along the c * -axis (see figure 4 (right panel) for picene and figure 5 (right panel) for coronene). In the calculated spectra the onset is at 3.3 eV for picene and 2.9 eV for coronene. This implies a large exciton binding energy in both cases: 0.8 eV for picene and 1.1 eV for coronene. The position of the main peaks is in agreement with experiment within 0.1-0.2 eV. Also in the calculations coronene turns out to be less anisotropic than picene, even though in the theoretical results for both compounds the changes in peak intensities between the different q directions seem to be enhanced with respect to experiment.
In figures 6 we compare the experimental and theoretical BSE spectra of the loss function for q directions along the a * -and b * -axis of tetracene (see solid lines). The overall agreement between the two is good. Along the a * -axis theory underestimates the position of the peaks by ∼ 0.1 eV. However, the shape of the spectrum along the a * is as in the experiment, with a first prominent peak followed by two smaller peaks and a shoulder at higher energy. These additional structures in the past have been interpreted as vibrational satellites of the main excitonic peak (see e.g. [49]). Here we see instead that they arise from other bound excitons (the calculated band gap is 2.9 eV). In fact, they have mainly an electronic origin, since in the theoretical model the coupling with vibrational modes is absent.
We have calculated also the absorption spectra (see dashed lines in the right panel of figure 6), given by the imaginary part of the dielectric function 2 in the dipole limit.
By comparing the EELS and the absorption spectra (solid versus dashed lines) we find that the two are pretty similar. This implies that in the low energy part of the loss function spectrum the screening effect from the real part of the dielectric function 1 is not affecting much the position of the structures seen in the loss function of tetracene (see equation (1)).
Comparing the spectra along a * and b * we find, similarly as for picene and coronene, that the anisotropy is enhanced in the calculations with respect to experiment. In particular, in the measured spectra there are three peaks both along a * and b * , while in the calculated spectra along b * only two peaks are present. The first peak in the experimental spectrum along b * , located at a position corresponding to the main peak along a * , is in fact absent in the calculations.
Possible explanations for these differences will be investigated more deeply in section 4.4, after the discussion on the nature of these lowest-energy excitons in section 4.3.

Exciton character
We can use the results of the diagonalization of the excitonic Hamiltonian equation (3) to analyse the lowest-energy excitons. Their nature is the result of the competition between e-h exchange, the band structure (both tend to delocalize the exciton) and the e-h direct attractive term (which tends to localize the exciton) [16]. In the case of picene we have shown in [16] that the lowest-energy exciton is mainly a Frenkel (FR) exciton, localized on a single molecule. This is contrast to pentacene, where the lowest-energy exciton is mainly a delocalized charge-transfer (CT) exciton [16,[50][51][52]. We now address this issue in the case of tetracene and coronene.
In a first approximation we take into account only transitions between highest occupied-lowest unoccupied molecular orbitals (HOMO-LUMO) in which we set to zero the band dispersion and we neglect the overlap between wavefunctions localized on different molecules. In this situation the excitonic Hamiltonian equation (3) can be described using a simple two-level tight-binding model [16]. This is very helpful to understand the nature of the excitons, since the excitonic Hamiltonian has thus only two possible solutions: a FR or a CT exciton. In the former case, the e-h pair is localized on the same molecule, while in the latter it is delocalized over different molecular units. Due to the presence of two inequivalent molecules in the unit cell, each FR and CT state and their combinations split in two components (this effect is known as Davydov splitting [18]) that are symmetric (+) and antisymmetric (−) with respect to the exchange of the e-h pair between inequivalent molecules. The energy of FR excitons is Here is the HOMO-LUMO gap, W is the on-site term of the direct e-h interaction W in equation (3), while I and J are the excitation transfer integrals stemming from the exchange e-h interactionv in equation (3) and are related to scattering processes of an e-h pair between two equivalent and inequivalent molecules respectively. The energy of CT excitons are given on the other hand by: whereW are the attractive inter-site terms of the direct e-h interaction W in equation (3). We note that the symmetric and antisymmetric CT excitons are degenerate, while for the FR excitons the Davydov splitting stems from the J contribution from the exchange e-h interactionv. Comparing the FR and CT solutions (8)-(9), we find that the condition for the CT state to be the lowest-energy excitation is The on-site W term is always larger than the inter-site termW. However, when the size of the molecular units is not too small, the two terms can become comparable. If the terms I − |J | are large enough, the CT exciton can be the lowest-energy solution. This is what happens for pentacene, but not for picene [16]. When the finite width of the bands is taken into account, the two CT and FR solutions are mixed together: in a real system with a non-zero band dispersion, the lowest excitation is always a mixture of the two excitons. The result of the ab initio BSE calculation (including contributions from all the bands besides the HOMO-LUMO transition) for tetracene and coronene is shown in figure 7 where we plot the electronic distribution of the excitonic wavefunction equation λ (r h , r e ) (4) of the lowest excited state for a fixed position of the hole r h (represented by the blue ball in figure 7). From out calculations we thus see that in both tetracene and coronene the first exciton is delocalized up to the first nearest-neighbour molecular units, as it occurs also in pentacene. Also in these cases therefore we find a CT exciton and not a FR exciton.
Comparing tetracene and pentacene one finds that the exciton is more localized on tetracene than in pentacene [51]. This is a consequence of the smaller size of the tetracene molecule with respect to pentacene (four benzene rings instead of five) and can be understood in terms of our simplified model. In fact, as the number of benzene rings decreases, the difference between W andW increases, and the FR and CT states become more separate in energy. As a consequence, the hopping processes related to the electronic band dispersion are less efficient in mixing the two excitonic states and, in turn, the excitonic wave function results more localized. Interestingly, comparing the excitons of coronene and picene we find that, although the binding energy is larger in coronene, the exciton in coronene is more delocalized than in picene, where it is pure FR exciton. This shows that the knowledge of the exciton binding energy is not sufficient  to understand the exciton character. The localization nature of the exciton in fact is set by the energy difference between FR and CT states, which is due to the competition between direct and exchange e-h interaction (see equation (10)), and by the efficiency of the hopping terms in mixing FR and CT states.

Loss map in the a * b * plane of tetracene
In order to investigate in more detail the effect of the anisotropy in the spectra of molecular solids, we have obtained the full map of the loss function in the a * b * plane of tetracene in correspondence with the lowest energy excitations (compare figures 8 and 9 (left panel)). In this way it is possible to monitor the evolution of the spectra by rotating the q direction from a * to b * . In figure 9 (right panel) we also show the same map for the calculated absorption spectrum: it has a similar behaviour as the loss spectrum map. Both in experiment and theory the spectra evolve continuously between the two q directions, fully displaying the anisotropic character of the response function in tetracene. The maps of the loss function in the a * b * plane clearly illustrate the high sensitivity of the loss spectra of molecular solids with respect to their crystal structure.
The anisotropy of the spectra is related to the Davydov splitting that we have discussed in the previous section and can be understood in terms of a molecular picture. In this framework, the dipole matrix elements for the exciton ± λ can be generically written in terms of the dipole matrix elements of the isolated molecule, which we callx andŷ for the HOMO-LUMO transition for a polarization direction parallel to the short (x) and the long (y) axis of tetracene [16] where α A and β A are the projections of r along the x and y axes of the molecule A (and, analogously, α B and β B for the second molecule B in the unit cell). In the present case, due to the symmetry of the tetracene molecule, theŷ dipole is zero. Therefore, the dipole matrix element defined in equation (11) simplifies into p ± = ± λ |r|0 = α Ax ± α Bx = p A ± p B . We see that the dipoles associated to the two Davydov components are perpendicular each other (p A + p B is perpendicular to p A − p B ). Moreover, in tetracene the angle between b * and the x axis is approximately the same for the two molecules, i.e. |α A | |α B |. Thus one of the two Davydov components has its dipole moment mainly oriented along b * (i.e. it is visible along b * ). Since a * and b * are perpendicular each other in the crystal structure of [26], the other Davydov component is polarized mainly along a * . From this simple analysis, we conclude that the two structures at 2.30 eV and 2.57 eV, which set the onset of the calculated spectrum along the a * -and b * -axis, respectively, are the two Davydov components of the lowest exciton. We also understand that in all the directions that are intermediate between a * and b * , the weight of the two Davydov components changes according to the projection of the polarization direction with respect to the directions of the two dipoles p + and p − .
In the EELS experiment the crystal structure was not fully accessible: in particular, it was not possible to precisely determine the mutual orientation of the molecules. In any case, we know that in the experimental sample the angle between a * and b * is slightly different from 90 • (see table 1). These are the first indications to explain the main discrepancies that resulted in the comparison between calculated and measured spectra in section 4.2. From our previous analysis, we in fact understand that in the experiment in the spectrum along the b * -axis the two lowest peaks arise from the two Davydov components, while for a momentum transfer parallel to a * axis we only observe the lowest Davydov feature. The same conclusion was obtained also for the similar case of pentacene [21]. In the theoretical spectra in figure 9 we see that the first Davydov component, missing along b * , starts to appear already at an angle of 10 • from the b * direction.

Temperature effects
Finally, to explain the differences in the relative intensities of the peaks between experiment and theory, it is important to consider also the effect of temperature. In the calculations the atoms are kept frozen in the crystal structure determined from diffraction data from [26] that are measured at a given temperature (155 K in the present case). In the EELS experiment, instead, the temperature modifies the crystal structure of the samples and affects the coupling with the vibrational modes of the crystal. In figure 10 we compare the loss spectra measured along a * and b * as a function of temperature. While the temperature is not affecting much the position of the peaks (less than 0.1 eV from 20 to 350 K, as it can be seen from figure 11), the shape of  the spectra changes considerably. By increasing the temperature, the various spectral features become broader and the differences between the intensities of the peaks become smaller. This analysis thus provides a second indication to explain the main differences between theory and experiment.

Conclusion and perspective
In conclusion, we have performed a detailed comparison of the loss spectra of three prototypical molecular solids, obtained using state-of-the-art experimental tools and theoretical methods from many-body perturbation theory. We find a good agreement between theory and experiment, with a close correspondence in the position and the shape of the peaks. However, the calculated spectrum is upshifted by 0.1-0.2 eV, which can be considered as the error bar of the calculation. The main differences between the measured and the calculated spectra have been explained in terms of the small differences between the crystal structures of the measured samples and the structural model used in the calculations. The additional role of temperature should be taken into account in order to describe accurately the evolution of the intensities of the peaks in the spectra. Nevertheless, the exciton nature can be safely described by theoretical models that neglect the effect of temperature, which can be easily monitored experimentally. The combination of highly accurate experimental and theoretical analysis is thus a powerful tool to understand the electronic properties of this class of materials.