Modeling the electronic structure of organic materials: A solid-state physicist’s perspective

Density-functional theory (DFT) is one of the most successful methods for electronic-structure calculations. Its theoretical foundation lies upon the Hohenberg-Kohn theorems [71], defining the unique correlation between the potential and the density of a many-electron system, which can be used to represent all its ground-state properties. The most popular and efficient recipe to implement DFT is given by the Kohn-Sham (KS) equations [72], which map the many-body problem into an auxiliary system of non-interacting particles ruled by the Schrödinger equation

$$\left[-\dfrac{\nabla^{2}}{2}+v_{s}(\mathbf{r})\right]\phi_{i}(\mathbf{r})=\epsilon^{KS}_{i}\phi_{i}(\mathbf{r}).$$

(1)

The eigenfunctions of Eq. (1) are the states of the fictitious KS system of independent electrons that are used to express the electron density,

$$n(\mathbf{r})=\sum_{i=1}^{occ.}|\phi_{i}(\mathbf{r})|^{2},$$

(2)

and the eigenvalues represent their corresponding single-particle energies. Beside the kinetic-energy term, the Hamiltonian of Eq. (1) includes an effective potential which is given by the sum of three terms:

$$v_{s}(\mathbf{r})=v_{ext}(\mathbf{r})+v_{H}[n](\mathbf{r})+v_{xc}[n](\mathbf{r}).$$

(3)

The external potential, $v_{ext}$, describes the electron-nuclear Coulomb attraction, the Hartree potential, $v_{H}[n](\mathbf{r})$, the mean-field electron-electron coupling, and the exchange-correlation (xc) potential, $v_{xc}[n](\mathbf{r})$, includes the remaining interactions. While the form of the first two potentials is known, unfortunately, the same is not true for $v_{xc}[n](\mathbf{r})$. Hence, the accuracy of the results provided by the KS equations crucially depends on the approximation chosen for this quantity.

The local-density approximation (LDA), originally proposed by Kohn and Sham in their seminal paper from 1965 [72], is based on the homogeneous electron gas model and represents the lowest rank of the so-called Jacob’s ladder of density-functional approximations [73]. The generalized gradient approximation (GGA) [74] is the next rung, with $v_{xc}$ depending on both the density and its gradients. Higher-level approximations include, on the one hand, the so-called metaGGA functionals [75, 76], and, on the other hand, hybrid functionals, where $v_{xc}$ is enhanced by a fraction of Fock exchange [77]. Among the latter, global hybrids are distinguished from the range-separated ones, in which the long- and short-range parts of the Coulomb potential are augmented separately [78, 79]. Interestingly, the diffusion of these families of functionals has followed separate paths in chemistry and physics, reflecting the substantially different performance of the various approximations on (isolated) molecules and solids. While for inorganic crystals the LDA and GGA have been extensively adopted for decades and specific flavors of range-separated hybrid functionals such as HSE [80] have gained popularity only recently, in the study of molecules the use of hybrid functionals has become a standard since the turn of the turn of the century.

The definition of the basis set for the KS equations (Eq. 1) is another crucial aspect in the implementation of DFT applied to non-periodic and extended systems. The delocalized character of the electronic wave-functions in solids is prone to a representation in terms of pure plane-waves or to plane-waves augmented by local orbitals [81]. The solution of the KS problem in k-space typically exploits the crystal periodicity and enables the application of efficient parallelization schemes [82]. On the other hand, the non-periodic nature of the electronic states in molecules calls for a real-space description with localized basis sets. Gaussian functions, atom-centered orbitals, and grid representations are the most common options ensuring an optimal trade-off between accuracy and computational costs [83, 84, 85].

The efficient ab initio scheme provided by DFT is limited by construction to the description of ground-state properties. In order to access excited states, it is necessary to take other routes. One option is to adopt the time-dependent extension of DFT (TDDFT) [86]. This approach is routinely applied both in a perturbative flavor [87] as well as by propagating the KS equations in real time [88]: the former approach gives access to the single-particle composition of the excitations but is limited to the linear response; the latter grants access to the response of all orders [89, 90] even in presence of an external, time-dependent electric fields [91, 92], but does not provide any information about the nature of the excited states. In both flavors, TDDFT results are crucially dependent of the choice of the exchange-correlation functional. Calculations based on hybrid functionals are largely applied in the chemistry community to describe optical excitations of molecules [69], even in conjunction with embedding methods [93] to account for solvation effects. On the other hand, in spite of enormous efforts in the past decades to develop appropriate approximations [94], TDDFT is not considered reliable to describe excitations in solids, not even when composed of molecules [95, 96]. For this reason, its popularity in the solid-state physics community is relatively low.

The state-of-the-art approach to describe electronic and optical excitations in extended systems is many-body perturbation theory (MBPT) [97], including the $GW$ approximation for the self-energy [98] and the solution of the Bethe-Salpeter equation [99]. Modern implementations of MBPT [100, 101] are built on top of DFT in the usual frameworks outlined in Section 2.1 for both periodic (see, e.g., Refs. [102, 103]) and non-periodic systems [104, 105, 106]. In this framework, the single-particle Green’s function

$$G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\sum_{i}\dfrac{\phi_{i}^{*}(\mathbf{r})\phi_{i}(\mathbf{r}^{\prime})}{\omega-\epsilon^{KS}_{i}-i\eta}$$

(4)

is defined in terms of the non-interacting KS states and energies obtained from Eq. (1). The expression of the dynamically screened Coulomb interaction

$$W(\mathbf{r},\mathbf{r}^{\prime},\omega)=\int\text{d}^{3}r_{1}\,\varepsilon^{-1}(\mathbf{r},\mathbf{r}_{1},\omega)\;v_{C}(\mathbf{r}_{1},\mathbf{r}^{\prime})\;,$$

(5)

in addition to the bare Coulomb potential $v_{C}$, includes the frequency-dependent inverse dielectric function $\varepsilon^{-1}$, which is evaluated from KS states, too [97]. Eqs. (4) and (5) are the ingredients for the construction of the self-energy in the $GW$ approximation, $\Sigma=iGW$. In the original formulation by Lars Hedin [98], this equation was supposed to be solved self-consistently as part of a closed loop of equations. In practice, the perturbative “single-shot” approach $G_{0}W_{0}$ is typically adopted with good results especially for conventional semiconductors [107]. In this framework, the quasi-particle (QP) equation dressing the KS energies with the self-energy contribution reads [108]:

$$\epsilon_{i}^{QP}=\epsilon_{i}^{KS}+Z_{i}\langle\phi_{i}|\Re\Sigma(\epsilon_{i}^{KS})-v_{xc}|\phi_{i}\rangle,$$

(6)

where $Z_{i}$ is a renormalization factor compensating for the evaluation of $\Sigma$ at $\epsilon_{i}^{KS}$, rather than the correct $\epsilon_{i}^{QP}$, by means of a linear extrapolation.

Starting from the QP-corrected electronic structure, optical excitations are computed by solving the Bethe-Salpeter equation (BSE) [99], the equation of motion for the two-particle correlation function [109]. By construction, this formalism enables the description of excitons quantitatively accounting for electron-hole interactions. In practice, the problem is expressed in terms of an effective two-particle Schrödinger equation

$$\sum_{o^{\prime}u^{\prime}}\hat{H}^{BSE}_{ou,o^{\prime}u^{\prime}}A^{\lambda}_{o^{\prime}u^{\prime}}=E^{\lambda}A^{\lambda}_{ou},$$

(7)

where $o$ and $u$ stand for occupied and unoccupied states, respectively. In spin-unpolarized systems, the BSE Hamiltonian is expressed as:

$$\hat{H}^{BSE}=\hat{H}^{diag}+\hat{H}^{dir}+2\hat{H}^{x},$$

(8)

where the diagonal term, $\hat{H}^{diag}$, accounts for the energy differences between occupied and unoccupied states, $\hat{H}^{dir}$ corresponds to the direct Coulomb attraction between the positively-charged hole and the negatively-charge electron,

$$\hat{H}^{dir}=-\int\text{d}^{3}r\int\text{d}^{3}r^{\prime}\,\phi_{o}(\mathbf{r})\phi^{*}_{u}(\mathbf{r}^{\prime})\,W(\mathbf{r},\mathbf{r}^{\prime})\,\phi^{*}_{o^{\prime}}(\mathbf{r})\phi_{u^{\prime}}(\mathbf{r}^{\prime}),$$

(9)

and the exchange term $\hat{H}^{x}$ includes the exchange electron-hole repulsion,

$$\hat{H}^{x}=\int\text{d}^{3}r\int\text{d}^{3}r^{\prime}\,\phi_{o}(\mathbf{r})\phi^{*}_{u}(\mathbf{r})\,\bar{v}_{C}(\mathbf{r},\mathbf{r}^{\prime})\,\phi^{*}_{o^{\prime}}(\mathbf{r}^{\prime})\phi_{u^{\prime}}(\mathbf{r}^{\prime}).$$

(10)

Eq. (9) includes the screened Coulomb interaction (Eq. 5), with $\varepsilon$ evaluated at $\omega=0$ (static screening). In Eq. (10), $\bar{v}_{C}$ is the short-range part of the bare Coulomb potential accounting for local-field effects which are known to be particularly relevant in the optical response of organic materials [95, 96, 110, 111].

The eigenvalues of Eq. (7), $E^{\lambda}$, represent excitation energies and mark the resonances in the absorption spectra. In the context of inorganic semiconductors, where excitons are loosely bound and largely delocalized, exciton binding energies, $E_{b}$, are defined as the difference between excitation energies and the QP optical gap: $E_{b}=E^{\lambda}-E^{QP-opt}_{gap}$ [112]. In the spectra of these materials, excitons are identified by narrow resonances below the onset [113]. This definition turns out to be inappropriate for organic crystals where the optical absorption spectrum is usually formed by distinct peaks even beyond the absorption onset, which are present already in the (joint) density of states [96, 114, 115]. In this scenario, a more consistent definition of the exciton binding energy is given by $E_{b}=E^{\lambda}_{BSE}-E^{\lambda}_{IQPA}$, where the independent QP energies, $E^{\lambda}_{IQPA}$, are the eigenvalues of Eq. (8) with $\hat{H}^{BSE}\equiv\hat{H}^{diag}$.

The eigenvectors of Eq. (7), $A^{\lambda}$, contain information about the electronic composition and the spatial distribution of the excited states. They act as weighting factors in the transition dipole coefficients, defined as $\mathbf{t}^{\lambda}=\sum_{ou}A^{\lambda}_{ou}\langle o|\widehat{\mathbf{d}}|u\rangle$ for localized systems and in terms of the momentum operator,

$$\mathbf{t}^{\lambda}=\sum_{ou\mathbf{k}}A^{\lambda}_{ou\mathbf{k}}\frac{\langle o\mathbf{k}|\widehat{\mathbf{p}}|u\mathbf{k}\rangle}{\varepsilon^{QP}_{u\mathbf{k}}-\varepsilon^{QP}_{o\mathbf{k}}},$$

(11)

in periodic structures. The imaginary part of the macroscopic dielectric function

$$\text{Im}\,\varepsilon_{M}=\frac{8\pi^{2}}{\Omega}\sum\limits_{\lambda}|{\mathbf{t}^{\lambda}}|^{2}\delta(\omega-E^{\lambda}),$$

(12)

describes the optical absorption of the material within the unit cell volume $\Omega$. Given the vectorial character of $\mathbf{t}^{\lambda}$, $\varepsilon_{M}$ is a tensor with as many non-vanishing components as enabled by the crystal symmetries.

Furthermore, the eigenvectors of Eq. (7) are included in the definition of hole and electron densities, defined as

$$\rho_{h}^{\lambda}(\textbf{r})=\sum_{ou}|A_{ou}^{\lambda}|^{2}|\phi_{o}(\textbf{r})|^{2},$$

(13)

and

$$\rho_{e}^{\lambda}(\textbf{r})=\sum_{ou}|A_{ou}^{\lambda}|^{2}|\phi_{u}(\textbf{r})|^{2},$$

(14)

respectively, for transitions between occupied states $\phi_{o}$ and unoccupied states $\phi_{u}$. These quantities were originally introduced for isolated systems [116, 117]. In solids, Eqs. (13) and (14) retain the periodicity of the lattice [115]: as such, they highlight the fictitiously periodic distribution of the hole and the electron but cannot provide any information about the exciton localization. For this purpose, it is necessary to introduce the exciton wave function, a six-dimensional quantity expressed by

$$\Psi^{\lambda}(\mathbf{r}_{h},\mathbf{r}_{e})=\sum_{uo\mathbf{k}}A^{\lambda}_{uo\mathbf{k}}\phi^{*}_{o\mathbf{k}}(\mathbf{r}_{h})\phi_{u\mathbf{k}}(\mathbf{r}_{e}),$$

(15)

and typically visualized by the isosurface of its square modulus with the hole or the electron position fixed.

The description of electronic structure provided by DFT (and even MBPT) assumes the systems to be in vacuo. Especially for molecules, this is hardly a realistic scenario. Molecules are usually emedded in an environment which has a strong influence on their electronic and optical properties. In order to capture these effects, implicit methods have been developed in quantum chemistry. The polarizable continuum model (PCM) is one of the most popular ones and describes the surrounding medium by a single polarizable continuum interacting with the charge density of the molecule hosted in a cavity. The validity of this approximation depends on the system. Since it fails to capture intermolecular interactions between static dipoles, van der Waals forces, or chemical bonding, including charge transfer, its accuracy is limited to those scenarios in which these effects are negligible, unless corresponding (usually empirical) corrections are included [118].

From a mathematical point of view, the polarization response of the medium described by the PCM is included through a density-dependent interaction term added to the KS Hamiltonian, $\hat{H}^{KS}[n]$, formulated in Eq. (1): $\hat{H}^{KS}[n]\rightarrow\hat{H}^{KS}[n]+\hat{v}^{\mathrm{PCM}}[n]$. The interaction term, $\hat{v}^{\mathrm{PCM}}[n]$, is a local electrostatic potential generated by the polarization of the surrounding. It can be exactly represented as arising from a surface charge density $\sigma[n]$,

$\displaystyle v^{\mathrm{PCM}}[n](\textbf{r})=\int_{\cal C}\text{d}^{2}s\,v_{c}(\textbf{r}-\textbf{s})\sigma[n](\textbf{s}),$

(16)

where $v_{c}(\textbf{r})=1/|\textbf{r}|$ is the bare Coulomb interaction and $\cal{C}$ is the cavity surface, i.e., the interface between the molecule and the medium. In turn, $\sigma[n]$ is determined from the electrostatic potential created by the molecule’s charges on $\cal{C}$:

$\displaystyle\sigma[n](\textbf{s})=\int_{\cal{C}}\text{d}^{2}s^{\prime}\,{\cal Q}(\textbf{s},\textbf{s}^{\prime})\left[\sum_{j=1}^{N}Z_{j}v_{c}(\textbf{s}^{\prime}-\textbf{R}_{j})-\int\text{d}^{3}r\,v_{c}(\textbf{s}^{\prime}-\textbf{r})n(\textbf{r})\right],$

(17)

where $Z_{j}$ and $\textbf{R}_{j}$ are nuclear charges and positions, respectively. The PCM matrix, $\cal{Q}$, encodes all information about the response of dielectric environment. Within the flexible integral equation formalism [119], it is most generally given as

$\displaystyle{\cal Q}=\left[\left(2\pi{\cal I}-{\cal D}_{e}\right){\cal S}_{i}+{\cal S}_{e}\left(2\pi{\cal I}+{\cal D}_{i}^{*}\right)\right]^{-1}\left[{\cal S}_{e}{\cal S}_{i}^{-1}\left(2\pi{\cal I}+{\cal D}_{i}\right)+\left(2\pi{\cal I}+{\cal D}_{e}\right)\right],$

(18)

where $\cal{I}$ is the identity operator, and ${\cal S}_{i}$, ${\cal D}_{i}$ and ${\cal D}^{*}_{i}$, as well as ${\cal S}_{e}$ and ${\cal D}_{e}$ are integral operators, defined by their action on an function $u$ supported on ${\cal C}$ as [120]:

	$\displaystyle\left({\cal S}_{i}u\right)(\textbf{s})$	$\displaystyle=\int_{\cal C}\text{d}^{2}s^{\prime}\,v_{c}(\textbf{s}-\textbf{s}^{\prime})u(\textbf{s}^{\prime})$		(19a)
	$\displaystyle\left({\cal D}_{i}u\right)(\textbf{s})$	$\displaystyle=\int_{\cal C}\text{d}^{2}s^{\prime}\left[\hat{\textbf{n}}(\textbf{s}^{\prime})\cdot\nabla^{\prime}v_{c}(\textbf{s}-\textbf{s}^{\prime})\right]u(\textbf{s}^{\prime})$		(19b)
	$\displaystyle\left({\cal D}^{*}_{i}u\right)(\textbf{s})$	$\displaystyle=\int_{\cal C}\text{d}^{2}s^{\prime}\left[\hat{\textbf{n}}(\textbf{s})\cdot\nabla v_{c}(\textbf{s}-\textbf{s}^{\prime})\right]u(\textbf{s}^{\prime})$		(19c)
	$\displaystyle\left({\cal S}_{e}u\right)(\textbf{s})$	$\displaystyle=\int_{\cal C}\text{d}^{2}s^{\prime}\,W(\textbf{s},\textbf{s}^{\prime})u(\textbf{s}^{\prime})$		(19d)
	$\displaystyle\left({\cal D}_{e}u\right)(\textbf{s})$	$\displaystyle=\varepsilon\int_{\cal C}\text{d}^{2}s^{\prime}\left[\hat{\textbf{n}}(\textbf{s}^{\prime})\cdot\nabla^{\prime}W(\textbf{s},\textbf{s}^{\prime})\right]u(\textbf{s}^{\prime}),$		(19e)

where $\hat{\textbf{n}}(\textbf{s})$ is the cavity surface normal vector and $\varepsilon$ is the dielectric constant of the material directly on the exterior of ${\cal C}$. The term $W(\textbf{r},\textbf{r}^{\prime})$ appearing in ${\cal S}_{e}$ and ${\cal D}_{e}$ is the screened Coulomb interaction in the external medium, disregarding the presence of the cavity. If the so-described environment is homogeneous with dielectric constant $\varepsilon$ (a reasonable approximation for solvents or crystalline environments), it is simply given as $W(\textbf{r},\textbf{r}^{\prime})=\varepsilon^{-1}v_{c}(\textbf{r}-\textbf{r}^{\prime})$. In this case, Eq. (18) can be simplified to

$\displaystyle{\cal Q}=-{\cal S}_{i}^{-1}\left(2\pi\frac{\varepsilon+1}{\varepsilon-1}{\cal I}-{\cal D}_{i}\right)^{-1}\left(2\pi{\cal I}-{\cal D}_{i}\right),$

(20)

see Ref. [121]. When instead the medium is anisotropic, adjustments to the formalism are needed, as discussed in Ref. [122] and in Section 4 below.

In the first step, we evaluate the frontier states of the isolated ANT molecule from DFT adopting a range-separated hybrid functional and applying the $GW$ approximation in the perturbative approach $G_{0}W_{0}$ (see computational details in the Appendix). In Fig. 1a), the isosurfaces of the highest occupied and lowest unoccupied molecular orbitals (HOMO and LUMO, respectively) are visualized in a schematic energy-level diagram including the neighboring levels, HOMO-1 and LUMO+1. The picture provided by these results, yielding a fundamental gap of 6.94 eV, is in excellent agreement with the experimental references for the electronic structure of this molecule (6.9 eV, see Ref. [123]).

In the bulk phase, ANT crystallizes in a monoclinic structure, space group $P2{{}_{1}}/a$  [124], including two inequivalent molecules in the unit cell, see Fig. 1c), where both the real-space representation of the material and its Brillouin zone (BZ) are shown. The resulting band structure, reported in the bottom panel of Fig. 1c), exhibits the hallmarks of organic crystals. First and foremost, the electronic states are no longer described by isolated energy levels but by bands with a finite k-dispersion providing information about the carrier mobility along a certain direction. The relatively flat character of the bands of ANT compared to their counterparts in inorganic semiconductors is related to the molecular nature of the constituents. Since in the energy window displayed in Fig. 1c) the electronic states have $\pi$/$\pi^{*}$ character, the dispersion is larger along those directions where the delocalization of the corresponding electronic cloud is maximized. For this reason, in the electronic structure of longer members of the oligoacene family, such as tetracene and pentacene, the band dispersion is overall enhanced [10, 13, 52, 125]. Another typical characteristic of organic crystals in the band structure of ANT is the twofold splitting of each electronic level, due to the presence of two inequivalent molecules in the unit cell. Depending on the specific direction in reciprocal space, the “sub-states” are either energetically slit up or almost degenerate: the former scenario is realized when the corresponding wave function is spread over both molecules in the unit cell, thereby inducing electronic repulsion; in the latter case, the wave function is localized on only one inequivalent molecule in the unit cell, and thus the mutual repulsion is minimized [52]. As we will see in Subsection 3.2, this characteristic impacts the optical spectrum of the material giving rise to Davydov splitting [52, 126]. The computed band-gap, 4.30 eV, indicates a substantial renormalization with respect to result obtained for the isolated molecule (6.94 eV, see Fig. 1a) as a consequence of the screening exerted by the molecules forming the crystal [127]. Notice that the aforementioned value, 4.30 eV, overestimates the experimental band gap of the ANT crystal, ranging between 3.90 and 4.00 eV [128, 129, 130]. We can ascribe this discrepancy to the fact that the band structure reported in Fig. 1c) is obtained adopting the GGA in DFT, with the QP correction to the band gap included by means of an empirical scissors operator, see details in the Appendix. This approach is often used when experimental references are available (see, e.g., Refs. [11, 52, 131]), given the high numerical costs of $GW$ calculations on these systems.

In order to understand whether the discussed electronic properties of the ANT crystal are primarily induced by short- or long-range effects, we investigate a non-periodic cluster formed by the two inequivalent molecules extracted from the unit cell of the crystal (see Fig. 1b). Neglecting the lattice periodicity, the energy levels are obviously dispersionless, confirming the long-range nature of this property. However, the twofold splitting in the energy levels of the molecule is preserved and its magnitude is consistent with the result obtained in the crystal at the high-symmetry point $Z$ (see Fig. 1c). This finding suggests that the short-range interactions between the molecules included in the unit cell are responsible for the Davydov splitting, hinting that this effect can be reproduced even by non-periodic models. The isosurface plots of the orbitals shown in Fig. 1b) reveal the correspondence between the HOMO and the HOMO-1 (LUMO and LUMO+1) of the cluster and the HOMO (LUMO) of the isolated molecule. In both manifolds, the lower-energy state is localized on the same molecule of the cluster. The orbital segregation is consistent with the character of the corresponding states in the crystal at the $\Gamma$ point. The value of the fundamental gap in the cluster, 6.66 eV, is less than 0.3 eV smaller than the one computed for the isolated molecule. This result, in line with previous findings obtained at the same level of theory for longer oligoacenes [132], is unsurprising and confirms that the band-gap reduction seen in the crystal cannot be captured with such a minimal model. Yet, as discussed in Section 4, a convenient shortcut is available for an accurate evaluation of this quantity at affordable computational costs, without the need to simulating the full periodic crystal.

We now turn to the optical properties of ANT, considering again the isolated molecule, the crystal, and the bimolecular cluster extracted from its unit cell. The absorption spectrum calculated for gas-phase compound (Fig. 2a) is characterized by an optically active excitation at the onset. The energy of this resonance is in excellent agreement with the experimental value [133] and the corresponding excitation is given by the HOMO$\rightarrow$LUMO transition polarized along the short molecular axis. This characteristic, common to all oligoacenes [12], can be understood from the parity of the frontier orbitals shown in Fig. 1a). Consistently, the hole and electron densities (Eq. 13 and 14) exhibit the same spatial distribution as the HOMO and the LUMO, respectively. The orientation of the transition dipole moment of this excitation explains its lower oscillator strength with respect to the higher-energy one, polarized along the long molecular axis, which is visible in the spectrum above 5 eV (see Fig. 1a).

Moving now to the result obtained for the bimolecular cluster (Fig. 2b), we notice essentially the same spectral features discussed above for the isolated moiety, except that in this case, the lowest-energy peak is formed by two almost degenerate excitations. Their energy separation, of the order of 10 meV in this system, corresponds to the Davydov splitting. The electron and hole densities of the two excitations giving rise to the first peak [see vertical bars in Fig. 2b), left panel] are localized on the same molecule. The lowest-energy excitation stems from the transition between the HOMO and the LUMO+1, while the second one from HOMO-1$\rightarrow$LUMO. Their polarization is along the short molecular axis, as in the isolated molecule. Excitations polarized along the long molecular axis appear at higher energies, above 5 eV.

The absorption spectrum calculated for the crystal is plotted showing two diagonal components of the dielectric tensor (see Fig. 2c). A sharp resonance along $b$ dominates the onset followed by a very weak excitation along $c$, in agreement with earlier results obtained at the same level of theory [11]. The presence of two excitations with two polarizations is due to Davydov splitting. However, in the crystal, this effect is not only related to the presence of two inequivalent molecules in the simulation cell, as in the cluster model, but also to their mutual orientation with respect to the crystalline axes. Appropriate rotation of the dielectric tensor along the electric field directions enhances the spectral intensity of the two Davydov components [52, 134]. The spatial distribution of the electron-hole pair corresponding to the lowest-energy exciton of the ANT crystal is shown as the square modulus of exciton wave function, see Eq. (15). In Fig. 2c), the isosurface represents the correlated electron distribution with respect to the fixed hole position, marked in the plot by a red dot. In analogy with its longer oligoacene siblings [12, 52], lowest-energy exciton of the ANT crystal has Frenkel character. Although this result may suggest that the adopting a cluster model is sufficient to capture the exciton nature in this system [compare Fig. 2b) and c)], the localization of the electron-hole pair in the crystal is actually a consequence of the large the electron-hole coupling strength, as extensively discussed in previous work on oligoacenes [12, 52, 60]. Neglecting the proper treatment of these effects, as provided by the solution of the BSE, would lead to an incorrect description not only of the optical spectra but especially of the exciton distribution, which would erroneously retain the periodicity of the electronic wave functions.

As a first example, we calculate the band-gap of the ANT crystal simulating atomistically only a single moiety and accounting implicitly for the surrounding molecules with the aid of the PCM. In this analysis, instead of computing the fundamental gap as the difference between LUMO and HOMO energies, we assume the definition $E_{gap}=\text{IP}-\text{EA}$. IP and EA are the ionization potential and electron affinity of the molecule, respectively, calculated as the total energy difference between the neutral species and its cation (ANT⁺) and anion (ANT^-): $\text{IP}=E_{tot}^{ANT^{+}}-E_{tot}^{ANT}$ and $\text{EA}=E_{tot}^{ANT}-E_{tot}^{ANT^{-}}$. This corresponds to the so-called $\Delta$SCF approach, which we apply herein adopting the GGA for the exchange-correlation functional (see computational details in the Appendix).

We start from the gas-phase molecule, for which we obtain IP $=7.07$ eV and EA$=0.70$ eV, both in good agreement with corresponding experimental data: measurements reported for IP range from 7.40 to 7.45 eV [135, 136], while for EA, reference values are between 0.53 and 0.66 eV [137, 138, 139]. The fundamental gap of gas-phase ANT is therefore equal to 6.37 eV (see Fig. 3a, right), i.e., slightly smaller than the experimental prediction (6.74 – 6.92 eV) extracted from the values of EA and IP reported above and reported in Ref. [123]. Discrepancies between calculated and measured values can be ascribed to the adopted GGA for the exchange-correlation functional.

Next, we apply the same procedure embedding the ANT molecule, its anion, and its cation into a PCM cavity filled by a homogeneous medium with $\varepsilon=3.10$, corresponding to the static dielectric constant computed for the ANT crystal in the random-phase approximation (RPA). In this setup, the calculated IP decreases to 5.79 eV and the EA increases to 1.84 eV (see Fig. 3a, right). It is interesting to notice that both quantities vary by an almost identical amount in comparison to their values in gas phase. The fundamental gap becomes 3.95 eV, in excellent agreement with available measurements for transport gaps in ANT [128, 129, 130]. It is worth noting that this good match is enabled by the low dispersion of the bands in the ANT crystal (see Fig. 1c), which is in turn an consequence of the localized character of electronic states of the constituents. In the presence of highly dispersive bands, where the energy levels assume very different values at varying k-vector, such a good agreement with an effective model for the environment is not foreseeable. Band dispersion is the consequence of precisely the type of intermolecular interactions that are neglected in the PCM and that require explicit quantum-mechanical modelling to be captured.

To prove that the results obtained for ANT are not just a favorable coincidence, we apply the same procedure to other short members of the oligoacene family. The results reported in Table 1 and the shown agreement with experimental references confirm the ability of this approach to correctly reproduce the band gaps of the crystals. It is worth noting that the static dielectric functions used to feed the PCM in these calculations have been computed from the RPA for the corresponding crystalline phases. However, this is not the only way to retrieve these quantities. A cheap alternative is to evaluate them from linear-response calculations on the gas-phase moiety and from the Clausius-Mossotti equation, which connects the microscopic polarizability of the single unit with the macroscopic dielectric constant of the crystal [140, 141]. Otherwise, for the considered systems, which are well-known and extensively characterized, corresponding empirical values are available in the literature [129, 142].

The generally good accuracy of the adopted $\Delta$SCF approach in the prediction of IP and EA makes viable the usage of optimally tuned, Koopman-compliant hybrid functionals [143]. Parameters therein are adjusted in a non-empirical way until the HOMO energy matches the $\Delta$SCF-determined IP; the EA can be accounted for as well, demanding that it agree with either the LUMO energy of the neutral molecule or - more strictly compliant with Janak’s theorem [144] - the HOMO energy of the anion. Actually, the IP and EA obtained with $\Delta$SCF are relatively robust with respect to the replacement of DFT with Fock exchange, meaning that hybrid functionals do not necessarily improve very much the estimate of these values with respect to pure DFT. However, hybrid functionals - particularly the long-range-corrected ones with a large amount of Fock exchange in the asymptotic limit - are essential when computing optical excitations with TDDFT, enabling acceptably accurate results even for charge-transfer excitation energies, which pure DFT notoriously struggles with [145]. Such hybrid TDDFT approaches coupled with the PCM can yield reasonable estimates for bandgaps and even for exciton binding energies [140, 141], in spite of neglecting static intermolecular interactions, and artificially confining excitons, e.g., into a single unit cell. These approximations can be softened by replacing a single molecule by a cluster, thus describing a larger portion of the system atomistically and quantum-mechanically. In this case, some pitfalls with respect to symmetry have to be avoided, which is generally lower for clusters than for both isolated molecules and the organic crystals [146].

Using the generic prescription of Eq. (18) for building the PCM matrix, one can model even more complex dielectric environments than an isotropic solvation cavity. The only requirement is an expression for the screened Coulomb interaction $W$. In Ref. [122], we extended the PCM to describe implicitly an arbitrary stack of layered materials intrinsically characterized by a high degree of anisotropy, introducing LayerPCM. This scheme reproduces the increasingly popular scenario of organic molecules adsorbed on two-dimensional transition-metal dichalcogenides (TMDC) [147, 148, 149, 150, 151, 152, 153]. In this case, the screened interaction in the region above the slab can be expressed as

$\displaystyle W(\textbf{r},\textbf{r}^{\prime})=v_{c}(\textbf{r}-\textbf{r}^{\prime})+W^{p}(\textbf{r},\textbf{r}^{\prime}),$

(21)

where the so-called polarization contribution, $W^{p}$, is related to the polarization of the substrate. For a semi-infinite substrate, $W^{p}$ can be written as the electrostatic potential of a single image charge within the substrate region. For a stacked substrate, it is more convenient to calculate it numerically as

$\displaystyle W^{p}(\textbf{r},\textbf{r}^{\prime})=\int_{0}^{\infty}\text{d}q_{\parallel}\,J_{0}(q_{\parallel}\Delta r_{\parallel})\frac{1-\varepsilon(q_{\parallel})}{1+\varepsilon(q_{\parallel})}e^{-q_{\parallel}(z^{\prime}+z-2z_{0})},$

(22)

using a Gauss-Laguerre quadrature [154]. In Eq. (22), $\Delta r_{\parallel}=[(x-x^{\prime})^{2}+(y-y^{\prime})^{2}]^{1/2}$, $J_{0}$ is the zeroth-order Bessel function of first kind [155], and $\varepsilon(q_{\parallel})$ is the effective nonlocal dielectric constant of the substrate, which can be determined using a transfer matrix formalism [122]. An alternative to this Bessel-function-based approach for obtaining $W^{p}$ is the recursive calculation of image charges [156].

We apply LayerPCM to model the hybrid interface formed by an ANT molecule adsorbed on a MoS₂ monolayer (see Fig. 3b). Three parameters are required to model the TMDC (or any other layered material) within LayerPCM: the layer thickness $t$, as well as the perpendicular and parallel components of the static dielectric constant, $\varepsilon_{\perp}$ and $\varepsilon_{\parallel}$, respectively. We employ a self-consistent scheme to obtain these quantities from first principles. The MoS₂ monolayer is simulated atomistically with periodic boundary conditions in all three directions, inserting a large amount of vacuum in the perpendicular direction to decouple the replicas. The dielectric constants are determined through linear-response calculations within the RPA. It should be noted that corresponding results are effective values, $\varepsilon_{\perp,\mathrm{eff}}$ and $\varepsilon_{\parallel,\mathrm{eff}}$, describing simultaneously the TMDC and the vacuum layer. In order to extract the dielectric constants of the material excluding the vacuum, we take an initial guess for $t$ and perform the following transformations:

	$\displaystyle\varepsilon_{\parallel}$	$\displaystyle=\left(\varepsilon_{\parallel}^{\mathrm{eff}}-1\right)\frac{c}{t}+1$		(23a)
	$\displaystyle\varepsilon_{\perp}$	$\displaystyle=\left[\left(\varepsilon_{\perp,\mathrm{eff}}^{-1}-1\right)\frac{c}{t}+1\right]^{-1},$		(23b)

where $c$ is the perpendicular lattice constant of the adopted supercell, i.e. the sum of the thickness of the TMDC and the height of the vacuum. The different formulas adopted for the in-plane and out-of-plane directions is due to the different boundary conditions for electric fields with parallel and perpendicular polarization with respect to the dielectric interfaces [157]. Equipped with these values, $t$ is obtained by connecting smoothly the plane-averaged semi-local $v_{xc}$ at short range with the image potential $V_{\mathrm{im}}$ of the MoS₂ layer at long range, corresponding to the correct asymptotic tail of the exact exchange-correlation potential [158]. This approach was employed in the past to determine mirror planes for metal surfaces [159]. For the dielectric monolayer, the image potential follows as a special case of Eq. (22) with $\textbf{r}=\textbf{r}^{\prime}$ and $\varepsilon(q_{\parallel})$ corresponding to a single layer:

	$\displaystyle V_{\mathrm{im}}(z)$	$\displaystyle=-\frac{1}{2}[(\varepsilon_{\perp}\varepsilon_{\parallel})^{1/2}+1][(\varepsilon_{\perp}\varepsilon_{\parallel})^{1/2}-1]\times$
		$\displaystyle\times\int_{0}^{\infty}\text{d}q_{\parallel}\,\frac{e^{-2q_{\parallel}z}}{1+\varepsilon_{\perp}\varepsilon_{\parallel}+2(\varepsilon_{\perp}\varepsilon_{\parallel})^{1/2}\coth(q_{\parallel}(\varepsilon_{\parallel}/\varepsilon_{\perp})^{1/2}t)}.$		(24)

The value obtained for $t$ by smoothly matching $v_{xc}$ and $V_{\mathrm{im}}$ is reinserted into Eq. (23) to calculate new dielectric constants. This procedure is iterated self-consistently. For MoS₂, we obtain $\varepsilon_{\parallel}=16.64$, $\varepsilon_{\perp}=11.25$, and $t=5.4~{}\mathrm{\AA}$. Owing to the prescription in Eq. (23b), $\varepsilon_{\perp}$ is extremely sensitive to changes in $t$. Indeed, values for $\varepsilon_{\perp}$ differing by a factor of 2 have been reported in Ref. [160] using a different definition of $t$. While such a wide range of results may seem alarming, it can be understood by realizing that the dielectric constant does not single-handedly describe the polarization response of the medium, but the geometry plays a crucial role, too. This is particularly true when the dimensions of the medium are small. Very different combinations of $\varepsilon_{\perp}$ and $t$ can thus give rise to similar physical scenarios.

Having defined this scheme, we evaluate the band gap of the physisorbed ANT molecule, where the underlying layer is accounted for via its dielectric function only. As in the crystalline environment, the presence of the TMDC reduces significantly the band gap of the moiety. This is a well-known effect that has been reported in a number of first-principles studies on hybrid interfaces formed by conjugated molecules adsorbed on TMDC monolayers [161, 162, 163, 164, 165, 166]. Comparing the results obtained for isolated ANT and for ANT adsorbed on a MoS₂ monolayer, we notice a reduction of IP and an increase of EA by about 0.8 and 1 eV respectively, which lead to a fundamental gap of 4.57 eV for the physisorbed molecule (see Fig. 3b, right panel). It is worth noting that the difference between IP and EA implies that already the relatively small ANT ions cannot be assimilated to point charges, for which both values would be equal. The band-gap renormalization for the molecule adsorbed on the TMDC is lower than in the crystal. Comparing the results reported in Fig. 3, panels a) and b), this behavior can be ascribed to more significant energetic readjustement of the HOMO level (set equal to -IP, according to Janak’s theorem) in the two scenarios. In contrast, the LUMO level, estimated as -EA, differs energetically by 150 meV in the implicit simulation of the crystalline environment and of the monolayer substrate.

To analyze more in depth the electronic properties of ANT physisorbed on the MoS₂ monolayer, we perform a fully atomistic benchmark calculation of the whole hybrid system simulating a supercell corresponding to 4$\times$4 the unit cell of the TMDC and with the molecule adsorbed flat at an approximate distance of 3.3 Å from the substrate. The electronic structure of the system, evaluated for consistency with the previous runs at the GGA level of DFT (PBE functional, see computational details in the Appendix), features a type-II band alignment between the organic and inorganic components (see Fig. 4). The HOMO of the molecule is situated within the energy gap of MoS₂, while the LUMO is higher up in the conduction band. We note in passing the clear signatures of mixing between the HOMO-1 of ANT and the wave functions of MoS₂ at the high-symmetry point $M$ in the valence region, which was found to be a hybridization hotspot between TMDCs and polycyclic aromatic hydrocarbons [167]. It is well known that DFT - particularly with semi-local functionals - is not very good at predicting energy level alignments: it notoriously underestimates band gaps as well as molecular IPs and, furthermore, does not capture image-charge induced level renormalization [168, 169]. Thus, we correct the frontier energy levels of the hybrid system as well as of the freestanding monolayer for reference with many-body perturbation theory, using the single-shot $G_{0}W_{0}$ method (Fig. 4b). As expected, this leads to a band-gap opening; the fundamental gap of MoS₂ is increased from 1.70 to 2.76 eV in the hybrid system, and to 2.77 eV in the freestanding configuration, in agreement with previously reported results on a similar level of theory [170, 171, 172, 173]. The presence of the molecule does not lead to a significant band-gap renormalization in the TMDC. The most important change brought about by the QP correction is the inverted order of the occupied energy levels. The severe underestimation of the ionization energy of ANT induced by DFT in the GGA is corrected by a large downshift of the HOMO, while the valence-band maximum of MoS₂ barely changes, resulting in a change from a type-II to type-I alignment with the molecular levels hugging those of the TMDC. Similar observations have been made for other hybrid interfaces [163, 165]. Hence, DFT alone yields unreliable results for the alignment of occupied states; for the unoccupied ones, however, the situation is not as dramatic, considering that the $GW$ corrections to the conduction-band minimum and to the LUMO differ only by 0.27 eV, compared to 0.89 eV between the valence-band maximum and the HOMO.

Interestingly, using the $\Delta$SCF+LayerPCM approach described above for calculating the renormalized band gap of ANT, we obtain results in excellent agreement with the $GW$-corrected ones, even adopting a GGA functional for DFT (Fig. 4b). Taking as a point of comparison the frontier levels of isolated MoS₂ computed from $G_{0}W_{0}$, $\Delta$SCF+LayerPCM yields a significantly better estimate for the level alignment of the hybrid interface than DFT alone. The only salient difference between the fully atomistic $GW$ calculation on the hybrid system and the effective treatment is given by a rigid shift of the MoS₂ levels, which is brought about by charge transfer [164]. Since the Fermi energy of ANT is higher than that of MoS₂, some electronic population is transferred from the organic to the inorganic side of the interface. The resulting partial charges - positive on ANT, negative on MoS₂ - slightly realign the levels of the subsystems with their corresponding electrostatic energy. This effect is obviously not captured when describing the two subsystems individually. However, the amount of charge transferred is based on the level alignment obtained from DFT in the GGA, as we apply only the perturbative $G_{0}W_{0}$ approach and not the self-consistent flavor of $GW$. The staggered line-up predicted by GGA can be assumed to lead to an overestimation of the difference between the individual Fermi energies and, thus, of charge transfer. Ongoing work to perturbatively include charge-transfer effects based on the energy levels predicted by calculations on the individual the subsystems is expected to further improve the accuracy of this method.

4T:TCNQ crystallizes in a triclinic lattice with unit cell containing a donor/acceptor pair with basal planes facing each other [181] (see Fig. 5a). The complex arrangement of the molecules with respect to the crystal structure and their relation with the BZ is schematically visualized in Fig. 5b), where the directions connecting $\Gamma$ with the other high-symmetry points are highlighted to ease the reading of the band-structure plot (Fig. 5c). By inspecting this result, we immediately notice highest-occupied and lowest-unoccupied states sticking out with respect to the manifold of valence and conduction bands, respectively. Interestingly, the character of these states varies at different k-points in the BZ. At the zone edges (high-symmetry points $Z$, $N$, $M$, and $R$), the wave functions are spatially segregated over the donor (valence states) and acceptor molecules (conduction states), respectively, see Fig. 5d). Notably, both valence-band maximum and conduction-band minimum are located at $M$. Elsewhere in the BZ, the electronic states are hybridized between the donor and the acceptor, in analogy with the HOMO and the LUMO of the corresponding molecular complexes [41, 44].

This counter-intuitive behavior of the frontier states in the co-crystal can be explained by the long-range nature of its electronic wave functions [115]. To this end, we introduce a tight-binding (TB) model assuming as a basis the molecular orbitals of the isolated bimolecular donor/acceptor cluster extracted from the crystalline unit cell. We assimilate the system to a quantum-mechanical two-level model including the highest-occupied and lowest-unoccupied energy levels, with corresponding wave functions expressed as

$$\psi_{\lambda\textbf{k}}(\textbf{r})=\sum_{m=H,L}C^{\lambda}_{m}(\textbf{k})\tilde{\varphi}_{m\textbf{k}}(\textbf{r}),$$

(25)

where $C^{\lambda}_{m}(\textbf{k})$ are the orthonormal expansion coefficients for either level $\lambda$, and

$$\tilde{\varphi}_{m\textbf{k}}(\textbf{r})=\sum_{j}\varphi_{m}(\textbf{r}+\textbf{R}_{j})e^{-i\textbf{k}\cdot\textbf{R}_{j}}.$$

(26)

In Eq. (26), $\varphi_{m}(\textbf{r})$ are the HOMO ($m=H$) and LUMO ($m=L$) of the isolated donor/acceptor complex, and $\textbf{R}_{j}$ are the lattice vectors. From Eqs. (25) and (26), it is evident that $\psi_{\lambda\textbf{k}}(\textbf{r})$ has the same periodicity of the co-crystal and thus fulfills Bloch’s theorem. It is worth highlighting the k-point dependence of the expansion coefficients which brings about the k-dependent character of the electronic states within each band seen in Fig. 5c). The wave functions and energies associated to the valence and conduction bands are obtained by diagonalizing the following problem:

$$\hat{H}^{TB}\begin{pmatrix}C_{H}^{\lambda}(\textbf{k})\\ C_{L}^{\lambda}(\textbf{k})\end{pmatrix}=\begin{pmatrix}A_{HH}(\textbf{k})&A_{HL}(\textbf{k})\\ A_{HL}(\textbf{k})&A_{LL}(\textbf{k})\end{pmatrix}\begin{pmatrix}C_{H}^{\lambda}(\textbf{k})\\ C_{L}^{\lambda}(\textbf{k})\end{pmatrix}=E_{\lambda}(\textbf{k})\begin{pmatrix}C_{H}^{\lambda}(\textbf{k})\\ C_{L}^{\lambda}(\textbf{k})\end{pmatrix},$$

(27)

where the nearest-neighbor (NN) approximation is assumed and the overlap integral is $\sum_{j}e^{-i\textbf{k}\cdot\textbf{R}_{j}}\int\varphi_{n}(\textbf{r})\varphi_{m}(\textbf{r}+\textbf{R}_{j})d^{3}r\approx\delta_{mn}$. Both assumptions are justified by the overall localized nature of the electronic states in molecular crystals. The matrix elements on the left-hand-side of Eq. (27) are calculated as

$$A_{nm}(\textbf{k})=A_{mn}(\textbf{k})=\sum_{j\in NN}e^{-i\textbf{k}\cdot\textbf{R}_{j}}\int\text{d}^{3}r\,\varphi_{n}^{*}(\textbf{r})\hat{H}^{TB}\varphi_{m}(\textbf{r}+\textbf{R}_{j})=t^{(0)}_{nm}+2\sum_{j\in NN}t_{nm}^{(j)}\cos(\textbf{k}\cdot\textbf{R}_{j}),$$

(28)

where $t_{nm}^{(j)}=\int\text{d}^{3}r\,\varphi_{n}^{*}(\textbf{r})\hat{H}^{TB}\varphi_{m}(\textbf{r}\pm\textbf{R}_{j})=\int\text{d}^{3}r\,\varphi_{n}^{*}(\textbf{r}\pm\textbf{R}_{j})\hat{H}^{TB}\varphi_{m}(\textbf{r})$ for $n,m\in\{H,L\}$ are the on-site ($\textbf{R}_{j}=0$) and hopping ($\textbf{R}_{j}=\textbf{R}_{a},\textbf{R}_{b},\textbf{R}_{c}$) integrals. The eigenvectors of Eq. (25) for the wave functions of the valence (VB) and conduction bands (CB) are expressed as

$$C^{VB}_{H}(\textbf{k})=w(\textbf{k}),\;\;\;\;\;C^{VB}_{L}(\textbf{k})=-w(\textbf{k})\mu(\textbf{k})$$

(29)

and

$$C^{CB}_{H}(\textbf{k})=w(\textbf{k})\mu(\textbf{k}),\;\;\;\;\;C^{CB}_{L}(\textbf{k})=w(\textbf{k}),$$

(30)

respectively, where

$$\mu(\textbf{k})=\frac{A_{HL}(\textbf{k})}{E_{(-)}(\textbf{k})+\sqrt{E^{2}_{(-)}(\textbf{k})+A^{2}_{HL}(\textbf{k})}},$$

(31)

with

$$E_{(-)}(\textbf{k})=A_{LL}(\textbf{k})-A_{HH}(\textbf{k})$$

(32)

and

$$w(\textbf{k})=\frac{1}{\sqrt{1+\mu^{2}(\textbf{k})}}.$$

(33)

Plugging Eqs. (29) and (30) into Eq.(25), we get:

$$\psi_{VB,\textbf{k}}(\textbf{r})=\sum_{j}e^{-i\textbf{k}\cdot\textbf{R}_{j}}w(\textbf{k})\left[\varphi_{H}(\textbf{r}+\textbf{R}_{j})-\mu(\textbf{k})\varphi_{L}(\textbf{r}+\textbf{R}_{j})\right],$$

(34)

and

$$\psi_{CB,\textbf{k}}(\textbf{r})=\sum_{j}e^{-i\textbf{k}\cdot\textbf{R}_{j}}w(\textbf{k})\left[\mu(\textbf{k})\varphi_{H}(\textbf{r}+\textbf{R}_{j})+\varphi_{L}(\textbf{r}+\textbf{R}_{j})\right].$$

(35)

The energies associated to these states are calculated from Eq. (27) as

$$E_{\lambda}(\textbf{k})=\sum_{j}\sum_{n,m}C^{\lambda}_{n}(\textbf{k})C^{\lambda}_{m}(\textbf{k})t^{(j)}_{nm}\cos(\textbf{k}\cdot\textbf{R}_{j})\;\;\;\;(\lambda=VB,CB).$$

(36)

With this formula, we can finally plot the energies of the VB and CB fit from the TB model (Eq. 36) with those output by the DFT calculations on the co-crystal, see Fig. 6, where we adopt a different k-path with respect to Fig. 5c) to ease visualization. The agreement between the two results is very good and confirms the validity of the model, which, therefore, can be used to rationalize the varying spatial distribution of the VB and CB wave functions in the BZ. To this end, we introduce the so-called segregation factor,

$$\mathcal{S}(\mathbf{k})=1-w^{2}(\textbf{k})\left(1-|\mu(\textbf{k})|\right)^{2}=\frac{2\mu(\textbf{k})}{1+\mu^{2}(\textbf{k})}\in[0,1],$$

(37)

and we plot it in Fig. 6 for a direct visualization of the band character. Large values of $\mathcal{S}(\mathbf{k})$, found in the vicinity of the high-symmetry points $M$, $R$ as well as $Z$ and $N$, correspond to spatial segregation of the wave function on either molecule. On the other hand, the small magnitude of $\mathcal{S}(\mathbf{k})$ close to $X$, $\Gamma$, $Y$, and $L$ is indicative of orbital delocalization across the donor/acceptor interface.

We can rationalize these results going back to the equations introduced above. From Eq. (31), it appears that $\mu(\mathbf{k})$ is maximized when $A_{HL}(\textbf{k})$ is large, i.e., when there is non-negligible mixing between the HOMO and the LUMO in the donor/acceptor complex. This happens along the crystallographic orientations that are substantially collinear to the lattice vectors aligned with the stacking direction of the molecules (see Fig. 5b). Since the off-diagonal elements of $\hat{H}^{TB}$ are non-negligible, the expansion coefficients of the wave-functions in Eqs. (29)-(30) have similar magnitudes. In this scenario, long-range interactions prevail over local couplings between donor and acceptor moieties, and induce the spatial segregation of the frontier states seen in Fig. 5d). On the other hand, along the directions that are approximately parallel to the molecular planes and thus orthogonal to the stacking direction, $A_{HL}(\textbf{k})\rightarrow 0$ and no mixing between HOMO and LUMO takes place. In this case, $\hat{H}^{TB}$ becomes diagonal, and one $C^{\lambda}_{m}(\textbf{k})$ coefficient dominates over the others. Consequently, $\mu(\mathbf{k})\rightarrow 0$ implying $\mathcal{S}(\mathbf{k})\rightarrow 0$: local couplings prevail over the long-range interactions, and the hybridized character of the frontier states is preserved in the periodic wave functions (Fig. 5d).

With the gained understanding of the electronic properties of the 4T:TCNQ co-crystal, we now analyze its optical excitations. Similar to ANT (Section 3 and Ref. [11]) and to organic crystals in general [7, 10, 52, 96, 110, 111, 13], the absorption spectrum of 4T:TCNQ is dominated by a strong resonance at the onset (see Fig. 7), on which we focus in the following. The pronounced anisotropy of the co-crystal is apparent in both the real and imaginary part of the dielectric tensor computed from the BSE [Fig. 7a), left and middle panels]. All the six inequivalent components contribute to the lowest-energy excitation. The one with the largest spectral strength is along $zz$, which is roughly aligned to the stacking direction of the molecules in the unit cell [see Fig. 5a)-b)]. This behavior is consistent with the first optical excitation of the isolated donor/acceptor complex [115], visualized in the spectrum in Fig. 7a), right panel. The main contribution to the lowest-energy excitation in the co-crystal comes from the transition between the VB maximum and the CB minimum, again in line with the behavior of the cluster where the first excited state corresponds to the HOMO$\rightarrow$LUMO transition [115]. However, as discussed above, the wave functions of the frontier states are segregated on the donor and acceptor molecules, respectively (see Fig. 5d). As a result, the first exciton in the co-crystal has charge-transfer character, as indicated in Fig. 7b) by the plots of the hole and electron densities, opposite to the Frenkel nature of its counterpart in the isolated complex (see Fig. 7c). In light of this behavior, the strong intensity of first excitation of the co-crystal may appear puzzling. To solve this conundrum, we make use again of the TB model introduced in Section 5.1.

Starting from the definition of Bloch’s states (Eq. 26), we peel off the periodic part

$$u_{\lambda}(\textbf{k},\textbf{r})=\sum_{j}\sum_{m=H,L}C^{\lambda}_{m}(\textbf{k})\varphi_{m}(\textbf{r}+\textbf{R}_{j})e^{-i\textbf{k}\cdot(\textbf{r}+\textbf{R}_{j})}$$

(38)

and use it to express the matrix elements for the optical transitions as

$$\langle u_{CB}(\textbf{k},\textbf{r})|\hat{i\nabla_{\textbf{k}}}|u_{VB}(\textbf{k},\textbf{r})\rangle=\sum_{m=H,L}C_{m}^{CB}(\textbf{k})(i\nabla_{\textbf{k}})C_{m}^{VB}({\textbf{k}})+\sum_{m,m^{\prime}=H,L}C_{m^{\prime}}^{CB}({\textbf{k}})C_{m}^{VB}({\textbf{k}})\textbf{d}_{m^{\prime}m},$$

(39)

where

$$\textbf{d}_{m^{\prime}m}=\int\text{d}^{3}r\,\varphi_{m^{\prime}}^{*}(\textbf{r})\,\hat{\textbf{r}}\,\varphi_{m}(\textbf{r})$$

(40)

is the transition dipole moment between the (localized) molecular orbitals of the donor/acceptor cluster extracted from the unit cell of the co-crystal. This term is not the only one entering the expression of the matrix elements. The first term on the right-hand side of Eq. (39), including the k-derivative of the expansion coefficients introduced in Eqs. (29)-(30), enables long-range couplings between the periodic states even when their spatial overlap is small. Thanks to this contribution, charge-transfer excitations like the lowest-energy one in the spectrum of the 4T:TCNQ co-crystal and of its siblings with partially or fully fluorinated acceptors [115] exhibit finite oscillator strength. Notice that, in contrast, charge-transfer excitations in non-periodic donor/acceptor organic interfaces are characterized by very weak spectral intensity [179, 180], due to the negligible spatial overlap of the wave functions involved in the transition.

The results reported in Section 3 for anthracene molecules and crystals are obtained through a two-step procedure. First, all structures are optimized by minimization of the interatomic forces until the threshold of 10^-3 eV/Å. In the optimization of the molecular clusters, the positions of C atoms are clamped in order to preserve the herringbone angle between the two molecules assigned in the input structure to mimic their arrangement in the crystal. Otherwise, the dimer relaxes to a configuration with the basal planes of the two molecules facing each other. These calculations are performed with the all-electron code FHI-aims [84] adopting tight integration grids, TIER2 basis sets [189], and the generalized-gradient approximation in the Perdew-Burke-Ernzerhof (PBE) parameterization [190] for the exchange-correlation potential. Van der Waals interactions are included via the pairwise Tkatchenko-Scheffler scheme [191].

In the second step, electronic and optical properties are computed for the obtained geometries. For the finite systems (isolated molecules and clusters), the MOLGW code [106] is employed with Gaussian-type cc-pVTZ basis sets [192] in conjunction with the frozen-core approximation and the resolution-of-identity approximation [193]. The range-separated hybrid functional CAM-B3LYP [194] is adopted in the DFT step as a starting point for the subsequent $G_{0}W_{0}$ and BSE runs. Calculations on the crystal are performed using the all-electron full-potential code exciting, implementing the family of linearized augmented plane-wave plus local orbitals (LAPW+LO) methods [81]. Muffin-tin radii of 1.05 bohr and 0.7 bohr are chosen for carbon and hydrogen, respectively, along with a plane-wave cutoff G_max=4.7 bohr^-1. The GGA in the PBE parameterization [190] is used for the exchange-correlation potential. In the DFT calculations, the BZ is sampled by a $8\times 6\times 5$ k-mesh, while in the BSE step [103] a $\Gamma$-shifted $8\times 6\times 5$ k-point mesh is adopted. The QP correction is added to the KS gap through a scissors operator equal to $\delta=$2.27 eV, extimated by aligning the lowest-energy excitation computed for the ANT crystal with available experimental data [195]. In the band-structure plot shown in Fig. 1c), the binding energy of the first exciton, $E_{b}=$1.30 eV, output by our BSE calculation, has been added to $\delta$, in analogy with the procedure adopted in Ref. [52]. In the solution of the BSE, the screened Coulomb potential $W$ is computed using 200 empty bands and the energy threshold for the local field effects is set to 1.0 Ha. In the construction and diagonalization of the BSE Hamiltonian, 10 occupied and 10 unoccupied bands are considered.

DFT calculations coupled with the PCM are performed with a locally modified copy of the Octopus code (v9.2) [196]. Nuclear potentials and core electrons are approximated by norm-conserving SG15 pseudopotentials [197], the exchange-correlation potential by PBE [190]. Neutral and singly-ionized molecules are treated in spin-restricted and unrestricted frameworks, respectively. KS wave functions are represented on a real-space grid generated by uniform sampling with spacing 0.21 Å of a union of spheres of radius 4.5 Å centered at each atomic position. The PCM cavity is generated in a similar, yet not identical fashion [198], using smaller spheres with atom-dependent radii. Here, we take the van der Waals radius of the respective element, scaled by a factor of 1.1.

For the atomistic calculation of the isolated TMDC monolayer and the ANT@MoS₂ interface as well as for the determination of the dielectric constants fed into the PCM, we exploit the interface between the Quantum Espresso suite (v6.8) [199] and the Yambo code (v5.0.2) [102, 200]. With the former, we obtain the DFT starting point for many-body calculations performed with the latter. We again employ SG15 pseudopotentials [197] for the cores and PBE for the exchange-correlation functional [190], neglecting spin-orbit coupling. Structural relaxations are done with a wave-function cutoff of 60 Ry and including the pairwise Tkatschenko-Scheffler dispersion correction [191], reducing the components of all forces to below 10^-3 Ha/bohr. At this step, we use k-point samplings of $8\times 8\times 1$, $2\times 2\times 1$, and $4\times 6\times 4$ for MoS₂, ANT@MoS₂, and the acene crystals, respectively. To simulate the ANT@MoS₂ interface, we choose a lattice constant of 20 Å in the non-periodic direction including a sufficiently large vacuum region to decouple the replicas, and we include a dipole correction. We use Wannier interpolation [201] as implemented in the Wannier90 code (v3.1.0) [202] to determine the density of states of ANT@MoS₂ on a $100\times 100\times 1$ k-grid, allowing for the application of only minimal artificial broadening in the electronic structure, which endows the obtained linewidths with physical meaningfulness. The corresponding band structure is unfolded to the BZ of the MoS₂ unit cell [167, 203]. The DFT calculations for the generation of the required large number of unoccupied orbitals in the many-body simulations are done with a wavefunction cutoff of 40 Ry, sticking with the same $4\times 6\times 4$ k-mesh for the acenes, but increasing it to $24\times 24\times 1$ and $6\times 6\times 1$ for MoS₂ and ANT@MoS₂, respectively, which is necessary to obtain converged $GW$ results. We calculate the dielectric functions of the acenes by means of linear-response calculation on the level of the RPA, featuring 300 conduction bands and a $\textbf{G},\textbf{G}^{\prime}$ sum kinetic cutoff of 3 Ry. The values used in the PCM correspond to the trace of the $\omega=0$ dielectric tensor, averaging $xx$, $yy$, and $zz$ components. For MoS₂, we take 100 conduction bands and a cutoff of 5 Ry to calculate the (effective) anisotropic dielectric functions within the RPA.

The $GW$ calculations on MoS₂ and ANT@MoS₂ are performed with the Yambo code [102, 200] adopting the plasmon-pole approximation for the dynamically screened interaction $W$. 300 (MoS₂) and 500 (ANT@MoS₂) conduction bands are included for both $W$ and the correlation part of the self-energy, $\Sigma_{c}$ and employ the sum-over-states terminator [204]. The $\textbf{G},\textbf{G}^{\prime}$ sum is cut off at kinetic energies of 5 Ry. A Coulomb cutoff [205] along the out-of-plane direction is applied, truncating the interaction at 19.57 Å $\lesssim$ ($c=20$ Å). Coulomb integrals over the BZ for the smallest $\sim$300 G vectors are calculated with the random integration method, replacing the sum over a uniformly sampled Brillouin zone by a 10⁶ q-point Monte-Carlo integration. Here, we neglect the anisotropy of $W$ for small q, choosing the corresponding polarization vector as $\frac{(1,1,1)}{\sqrt{3}}$, representing an orientational average.

The results for the 4T:TCNQ co-crystal presented in Section 5 have been obtained from DFT using the code Quantum Espresso [199] and from MBPT using Yambo [102, 200]. In the former step, a plane-wave basis set cutoff of 50 Ry (200 Ry) for the wave-functions (electron density) and norm conserving pseudopotentials [206] are employed. The PBE functional [190] is adopted in conjunction with the Tkatchenko-Scheffler pairwise scheme [191] to account for van der Waals interactions. A 4$\times$4$\times$4 k-grid is used to sample the BZ. Self-consistent calculations are carried out on optimized geometries obtained by minimizing interatomic forces with a threshold of 10^-4 Ry/Å. BSE calculations are performed on top of the DFT electronic structure with the quasi-particle correction added in the form of a scissors operator of 1.41 eV to the conduction bands based on available experimental spectra [44]. The BSE Hamiltonian is constructed including 20 occupied and 40 unoccupied states, and by 4$\times$4$\times$4 k-point grid; it is diagonalized using the Haydock-Lanczos algorithm [207].

Calculations on the 4T:TCNQ cluster are performed with MOLGW [106] adopting the same computational parameters employed for ANT, except for the Gaussian basis set, which is here reduced to aug-cc-pVDZ.