*Article* **Prediction of Strong Transversal s(TE) Exciton–Polaritons in C60 Thin Crystalline Films**

**Vito Despoja 1,2,\* and Leonardo Maruši´c 3**


**Abstract:** If an exciton and a photon can change each other's properties, indicating that the regime of their strong bond is achieved, it usually happens in standard microcavity devices, where the large overlap between the 'confined' cavity photons and the 2D excitons enable the hybridization and the band gap opening in the parabolic photonic branch (as clear evidence of the strong exciton–photon coupling). Here, we show that the strong light–matter coupling can occur beyond the microcavity device setup, i.e., between the 'free' s(TE) photons and excitons. The s(TE) exciton–polariton is a polarization mode, which (contrary to the p(TM) mode) appears only as a coexistence of a photon and an exciton, i.e., it vanishes in the non-retarded limit (*<sup>c</sup>* → <sup>∞</sup>). We show that a thin fullerene C<sup>60</sup> crystalline film (consisting of *N* C<sup>60</sup> single layers) deposited on an Al2O<sup>3</sup> dielectric surface supports strong evanescent s(TE)-polarized exciton–polariton. The calculated Rabi splitting is more than Ω = 500 meV for *N* = 10, with a tendency to increase with *N*, indicating a very strong photonic character of the exciton–polariton.

**Keywords:** eksciton-polaritons; molecular crystals; optical conductivity; photonics

#### **1. Introduction**

The interaction between photons and polarization modes can result in the formation of hybrid photon polarization modes, called polaritons [1]. Very common platforms for studying strong light–matter interactions are the gapped systems, such as semiconductors [2,3] or molecules [4], placed in microcavity devices, where the cavity exciton–polaritons are formed. The quantum nature of a cavity exciton–polariton manifests in the form of the Bose– Einstein condensation, which has recently been experimentally detected [2,5,6]. The cavity exciton–polaritons are routinely observed in bulk [7–10] and quantum well systems [3,11], e.g., devised from GaAs [3]. Two-dimensional (2D) materials, such as semiconducting monolayers, thin heterostructures, and films, are even more attractive than their bulk counterparts, due to the reduced Coulomb screening and the corresponding large exciton binding energies [12–19] that enable the formation of well-defined exciton–polaritons even at room temperatures [20]. The first 2D exciton–polaritons were obtained in a monolayer of a transition metal dichalcogenide (TMD) MoS2, where the Rabi splitting between the exciton and the cavity photon of ∼50 meV was observed [21]. Further photoluminescence studies showed clear anti-crossing behavior and splitting of the exciton–polariton in other 2D TMD cavity devices, e.g., in MoSe<sup>2</sup> [22], WS<sup>2</sup> [23], WSe<sup>2</sup> [24,25], and in the MoSe2-WSe<sup>2</sup> heterostructure [26]. In addition, the real-space imaging of the exciton–polaritons has been done by means of near-field scanning optical microscopy for WSe<sup>2</sup> thin films [27]. Finally, a remarkable Rabi splitting of 440 meV was theoretically predicted in the monolayer hexagonal boron nitride cavity device [28], suggesting extraordinary strong light–matter interaction in 2D heterostructures.

The strongest exciton–photon coupling is achieved in the organic dye molecule thin films placed in a microcavity [29]. For example, Rabi splitting of <sup>Ω</sup> ≤ 450 meV [30–32],

**Citation:** Despoja, V.; Maruši´c, L. Prediction of Strong Transversal s(TE) Exciton–Polaritons in C<sup>60</sup> Thin Crystalline Films. *Int. J. Mol. Sci.* **2022**, *23*, 6943. https://doi.org/ 10.3390/ijms23136943

Academic Editor: Ana María Díez-Pascual

Received: 24 May 2022 Accepted: 20 June 2022 Published: 22 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

0.7 eV [33], and even more than 1 eV [34,35] have been detected when various organic dye molecules were placed in a planar microcavity. The theoretical approach to such systems is mostly based on the two-level or boson–boson Hamiltonian model with an arbitrary coupling constant [4,36]. Using graphene [37] or perovskite [38]-layered heterostructures, one can obtain tunable microcavity devices of high performance, which can be applied as photonic detectors or emitters [38,39], but also as platforms for studying the exciton– photon coupling.

All these studies use the same concept: an exciton in a semiconducting nanostructure hybridizes with a photon 'confined' in a metallic microcavity. Such cavity photons and excitons are expected to interact stronger as the photon is more confined (i.e., the overlap between the exciton and photon is larger) and as the exciton oscillator strength [28] is larger. In this paper, we change the concept and explore the coupling between the 'free' photons and the excitons in the 2D nanostructures. The coupling between the 'free' photons and the polarization modes (such as plasmons, phonons, or excitons at surfaces or in 2D nanostructures) is a well-known, widely explored phenomenon [40–47]. The inherent property of all these eigenmodes (called plasmon polaritons, phonon polaritons, or exciton– polaritons) is their evanescent character, i.e., they are well-defined eigenmodes with the electromagnetic field (wave function) strongly localized at the interface or within the 2D nanostructure. Moreover, these modes usually have p or transverse-magnetic (TM) polarization, i.e., the electric field (and, therefore, the currents as well) has a component parallel to the direction of propagation. The most relevant point is that the p(TM) polarized plasmon polariton, phonon polariton, or exciton–polariton branches reduce to the plasmon, phonon, or exciton branches in the non-retarded (*<sup>c</sup>* → <sup>∞</sup>) limit. On the other hand, the electric field and the current of s or transverse-electric (TE) electromagnetic eigenmodes are perpendicular to the direction of their propagation. An especially attractive aspect of these polarization modes is that they **do not** exist in the non-retarded (*<sup>c</sup>* → <sup>∞</sup>) limit, i.e., they appear only as the coexistence of a photon and an exciton. The extent of the photon's participation in the s(TE) exciton–polariton is determined from the bending of the horizontal exciton branch (*ω*ex) at the exciton–photon crossing (*ω*ex = *Q*ex*c*/ √ *ǫ*, where *Q*ex is the photon wave vector at the exciton–photon crossing point), which we call the Rabi splitting Ω, to keep the terminology compatible with the cavity systems. Even though the s(TE) surface or 2D polaritons do exist [40] for some conditions, there is still no experimental evidence of such modes. However, the hybridization between the s(TE) Bloch surface waves (BSWs) (i.e., the photons confined between a truncated photonic crystal and a semi-infinite dielectric), and the excitons has been experimentally demonstrated in both inorganic (quantum well and TMD monolayer) and organic systems [48–51].

We show that very strong s(TE) exciton–polaritons may occur in layered van der Waals (vdW) heterostructures. The prototypical layered nanostructure we investigate in this paper is a thin film of the FCC fullerite (crystalline fullerene) cut along the (111) planes so that it formed several (*N*) molecular (C60) layers. The crystalline C<sup>60</sup> films were also deposited on a dielectric Al2O<sup>3</sup> surface to make the simulation more realistic. The epitaxial growth of the C<sup>60</sup> thin films of various thicknesses on various metallic or dielectric substrates under ambient conditions and in high vacuum was studied in references [52–58]. Some experimental studies even show that the crystalline growth of the C<sup>60</sup> thin films on pentacene buffer layers is exclusively (111)-oriented [55]. Theoretical, molecular dynamic simulations of the C<sup>60</sup> multilayer epitaxial growth and stability on various substrates were investigated in references [58–61]. These experimental/theoretical studies suggest that our model system is indeed highly realistic.

In this paper, the light–matter interaction was studied using our quantum electrodynamic Bethe–Salpeter equation approach (QE-BSE) developed in references [28,47]. This approach describes both excitons and photons by bosonic propagators *σ* and Γ, respectively, derived from the first principles. The C<sup>60</sup> optical conductivity *σ* was calculated using *ab initio* G0W0-BSE method [47,62], and the free proton propagator Γ was derived by solving Maxwell's equation at the vacuum/dielectric interface [63,64]. The dielectric surface was

described by the local dielectric function *ǫM*(*ω*), also determined from the first principles. The exciton–photon coupling was achieved by dressing the free-photon propagator Γ with excitons at the random phase approximation (RPA) level. We studied the s(TE)-polarized exciton–polariton in the C<sup>60</sup> thin crystalline film as a function of the number *N* of the C<sup>60</sup> single layers. We obtained a very strong hybridization between the exciton in the C<sup>60</sup> thin film and the s(TE) free photons. The hybridization increased with *N*, and for *N* = 10, we achieved the Rabi splitting Ω even larger than 1000 meV and 500 meV for self-standing and supported C<sup>60</sup> films, respectively.

The paper is organized as follows. In Section 2, we present the geometry of the system and the derivation of the optical conductivity *σ* ˜ of the C<sup>60</sup> single layer using the G0W0-BSE approach with the solution of Dyson's equation for the electric field propagator E = <sup>Γ</sup> + <sup>Γ</sup>*<sup>σ</sup>* ˜ E. In Section 3, we present the spectra of the electromagnetic modes *S* = ReE, as well as the dispersion relations and intensities of the s(TE)-polarized exciton–polaritons in the C<sup>60</sup> thin films of different thicknesses *N*, in vacuum, and at a Al2O<sup>3</sup> surface. The conclusions are presented in Section 4.

#### **2. Theoretical Formulation**

We assume the the C<sup>60</sup> molecules, upon deposition on the crystal surface, self-assemble in a regular FCC structure (the most stable bulk structure of crystalline fullerene) forming a (111) surface, as shown in Figure 1a. The FCC crystal lattice constant is taken to be *a*3D = 14 Å [61], and the separation between the layers is fixed to be ∆ = *a*3D/ √ 3 = 8.1 Å. Each crystal plane forms a 2D hexagonal Bravais lattice with the lattice constant *a*2D = *a*3D/ √ 2 = 9.9 Å, as shown in Figure 1b. The C<sup>60</sup> films, occupying *z* > 0 halfspace, are immersed in a dielectric medium described by a dielectric constant *ǫ*0. The dielectric response of the substrate, occupying *z* < 0 half-space, is approximated by a local macroscopic dielectric function *ǫ*M(*ω*).

**Figure 1.** (**a**) C<sup>60</sup> molecules upon deposition on the surface self-assemble in a regular FCC structure forming a (111) surface. The C<sup>60</sup> layers, occupying *z* > 0 half-space, are immersed in a dielectric medium described by a dielectric constant *ǫ*0. The dielectric response of the substrate, occupying *z* < 0 half-space, is approximated by a local macroscopic dielectric function *ǫ*M(*ω*). The FCC lattice constant is *a*3D = 14 Å so that the separation between the layers is ∆ = *a*3D/ √ 3 = 8.1 Å. (**b**) Each crystal plane forms a 2D-hexagonal Bravais lattice with lattice constant *a*2D = *a*3D/ √ 2 = 9.9 Å.

### *2.1. Calculation of Electric Field Propagator* E

The quantity we used to extract the information from the electromagnetic modes in the C<sup>60</sup> films deposited on a dielectric surface was the electric field propagator E*µν*, which provides the electric field produced by an external oscillating point dipole **p**<sup>0</sup> *e* <sup>−</sup>*iω<sup>t</sup>* placed at point **r** ′ as [63,64]

$$E\_{\mu}(\mathbf{r},\omega) = \sum\_{\nu=\mathbf{x},\mathbf{y},\mathbf{z}} \mathcal{E}\_{\mu\nu}(\mathbf{r},\mathbf{r}',\omega) p\_{\nu}^{0}. \tag{1}$$

The propagator E is the solution of Dyson's equation [28,47,64,65]

$$\mathcal{E}\_{\mu\nu}(\mathbf{r}, \mathbf{r}', \omega) = \Gamma\_{\mu\nu}(\mathbf{r}, \mathbf{r}', \omega) + \sum\_{\mathbf{a}, \beta = \mathbf{x}, \mathbf{y}, \mathbf{z}} \int d\mathbf{r}\_1 \int d\mathbf{r}\_2 \,\Gamma\_{\mu\alpha}(\mathbf{r}, \mathbf{r}\_1, \omega) \sigma\_{\mathbf{a}\beta}(\mathbf{r}\_1, \mathbf{r}\_2, \omega) \mathcal{E}\_{\beta\nu}(\mathbf{r}\_2, \mathbf{r}', \omega), \tag{2}$$

where the integration is performed over the entire space, *σ* is the nonlocal conductivity tensor of the deposited C<sup>60</sup> thin film, and Γ is the electric field propagator in the absence of the C<sup>60</sup> film, i.e., when *σ* = 0 [64,65]. The propagator Γ also includes the electromagnetic field scattering at the medium/substrate interface. If each molecule is approximated as a point polarizable dipole, then the optical conductivity of the C<sup>60</sup> film can be written as

$$\sigma\_{\mu\nu}(\mathbf{r}, \mathbf{r}', \omega) = \sum\_{i=1}^{N} \sum\_{\mathbf{R}\_{\parallel}} \sigma\_i^{\mu\nu}(\omega) \delta(\rho - \mathbf{R}\_{\parallel}) \delta(z - z\_i) \delta(\rho' - \mathbf{R}\_{\parallel}) \delta(z' - z\_i), \tag{3}$$

where *σ µν i* (*ω*) is the optical conductivity tensor of a single molecule in the *i*-th molecular layer. This approximation is fully justified in the optical limit 2*πc*/*ω*light >> *R*M, where *R*<sup>M</sup> is the radius of a C<sup>60</sup> molecule. Note that although all the molecules are equal, we distinguished between their conductivities in different layers *σ<sup>i</sup>* ; *i* = 1, *N*, due to the different influences of the substrates on a molecule in a different layer. The 2D Bravais lattice translation vectors spanning the molecular crystal are

$$\mathbf{R}\_{\parallel} = n\mathbf{a}\_1 + m\mathbf{a}\_2; \ n, m \in \mathbb{Z},\tag{4}$$

where **a**<sup>1</sup> and **a**<sup>1</sup> are the primitive vectors, as illustrated in Figure 1b. The molecular layers occupy the planes

$$z\_i = z\_0.z\_0 + \Delta \, z\_0 + 2\Delta \, \dots \, \, z\_0 + (N-1)\Delta \, \dots$$

where *N* is the number of molecular layers. Due to the planar translational invariance of the substrate, the Fourier transform of the propagator Γ is

$$
\Gamma\_{\mu\nu}(\mathbf{r}, \mathbf{r}', \omega) = \int \frac{d\mathbf{Q}}{(2\pi)^2} \Gamma\_{\mu\nu}(\mathbf{Q}, \omega, z, z') e^{i\mathbf{Q}(\mathbf{p} - \mathbf{p}')},\tag{5}
$$

where **Q** = (*Qx*, *Qx*) are the two-dimensional wave vectors. The propagator E should also include the effects of the electromagnetic field Bragg scattering at the 2D crystal lattice, so that its Fourier transform is

$$\mathcal{E}\_{\mu\nu}(\mathbf{r}, \mathbf{r}', \omega) = \sum\_{\mathbf{g}\_{\parallel}} \int \frac{d\mathbf{Q}}{(2\pi)^2} \mathcal{E}\_{\mathbf{g}\_{\parallel}}^{\mu\nu}(\mathbf{Q}, \omega, z, z') e^{i\mathbf{Q}\cdot\mathbf{p}} e^{-i(\mathbf{Q} + \mathbf{g}\_{\parallel})\mathbf{p}'},\tag{6}$$

where

$$\mathbf{g}\_{\parallel} = n\mathbf{b}\_1 + m\mathbf{b}\_2; \ n, m \in \mathbb{Z} \tag{7}$$

are 2D reciprocal vectors, while **b**<sup>1</sup> and **b**<sup>1</sup> are primitive reciprocal vectors. After inserting (3), (5), and (6) in (2), it becomes an equation in the (**Q**, *ω*, *z*) space

$$\mathcal{E}\_{\mathbf{g}\_{\parallel}}^{\mu\nu}(\mathbf{Q},\omega,z,z') = \Gamma\_{\mu\nu}(\mathbf{Q},\omega,z,z') + $$

$$\sum\_{\mathbf{x},\mathbf{f}=\mathbf{x},y,z} \sum\_{i} \sum\_{\mathbf{g}'\_{\parallel}} \Gamma\_{\mu\mathbf{a}}(\mathbf{Q},\omega,z,z\_{i}) \mathcal{e}\_{i}^{a\beta}(\omega) \mathcal{E}\_{\mathbf{g}\_{\parallel}+\mathbf{g}'\_{\parallel}}^{\beta\nu}(\mathbf{Q}-\mathbf{g}'\_{\parallel},\omega,z\_{i\nu}z'), \tag{8}$$

where the surface optical conductivity is

$$\mathfrak{d}\_{\rm i}^{a\beta}(\omega) = \frac{1}{S\_{\rm fcc}} \sigma\_{\rm i}^{a\beta}(\omega),\tag{9}$$

and *S*fcc = (**a**<sup>1</sup> × **a**2)**z**ˆ = *a* 2 2D √ 3/2 is the area of the 2D unit cell. If we neglect the electromagnetic field Bragg scattering, by introducing **<sup>g</sup>**<sup>k</sup> = **<sup>g</sup>** ′ <sup>k</sup> <sup>=</sup> 0 in Equation (8), and inserting *z* = *z<sup>i</sup>* and *z* ′ = *z<sup>j</sup>* , the equation becomes the matrix tensor equation:

$$\mathcal{E}\_{\mu\nu}(\mathbf{Q},\omega,z\_{i\prime}z\_{j}) = \Gamma\_{\mu\nu}(\mathbf{Q},\omega,z\_{i\prime}z\_{j}) + \sum\_{a,\beta=x,y,z} \sum\_{k} \Gamma\_{\mu\mu}(\mathbf{Q},\omega,z\_{i\prime}z\_{k}) \psi\_{k}^{\mathbf{a}\beta}(\omega) \mathcal{E}\_{\beta\nu}(\mathbf{Q},\omega,z\_{k\prime}z\_{j}), \tag{10}$$

where E *µν*(*z<sup>i</sup>* , *zj*) is the electric field propagator within (*i* = *j*) or between (*i* 6= *j*) the C<sup>60</sup> layers. The electrical field propagator in the absence of the C<sup>60</sup> film can be written as

$$
\Gamma = \Gamma^0 + \Gamma^{sc},
\tag{11}
$$

where the propagator of the 'free' electric field (or free photons propagator) is [63,64]

$$\Gamma^{0}(\mathbf{Q},\omega,z,z') = -\frac{4\pi i}{\varepsilon\_{0}\omega}\delta(z-z')\mathbf{z}\cdot\mathbf{z} - \frac{2\pi\omega}{\beta\_{0}c^{2}}e^{i\beta\_{0}|z-z'|}\sum\_{q=s,p}\mathbf{e}\_{q}^{0}\cdot\mathbf{e}\_{q}^{0}.\tag{12}$$

The propagator of the scattered electric field in the region *z*, *z* ′ > 0 is [63]

$$\mathbf{T}^{\rm sc}(\mathbf{Q},\omega,z,z') = -\frac{2\pi\omega}{\beta\_0 c^2} e^{i\beta\_0(z+z')} \sum\_{q=s,p} r\_q \cdot \mathbf{e}\_q^+ \cdot \mathbf{e}\_q^-. \tag{13}$$

Here, the unit vectors of the s(TE)-polarized electromagnetic field are

$$\mathbf{e}\_{\mathbf{s}}^{0,\pm} = \mathbf{Q}\_0 \times \mathbf{z}.$$

and the unit vectors of the p(TM) polarized electromagnetic field are

$$\mathbf{e}\_p^{0,\pm} = \frac{c}{\omega\sqrt{\epsilon\_0}} [\alpha\_{0,\pm}\beta\_0\mathbf{Q}\_0 + Q\mathbf{z}]\_{\prime}$$

where *α*<sup>0</sup> = −*sgn*(*z* − *z* ′ ), *α*<sup>±</sup> = ∓1, and **Q**<sup>0</sup> and **z** are the unit vectors in the **Q** and *z* directions, respectively. The reflection coefficients of the s(TE) and p(TM) polarized electromagnetic waves at the medium/substrate interface are

$$r\_{\mathbf{s}} = \frac{\beta\_0 - \beta\_{\mathbf{M}}}{\beta\_0 + \beta\_{\mathbf{M}}}$$

and

$$r\_p = \frac{\beta\_0 \epsilon\_\mathbf{M} - \beta\_\mathbf{M} \epsilon\_0}{\beta\_0 \epsilon\_\mathbf{M} + \beta\_\mathbf{M} \epsilon\_0}.$$

respectively, and the complex wave vectors in perpendicular (*z*) direction are

$$
\beta\_{0,\mathbf{M}} = \sqrt{\frac{\omega^2}{c^2} \epsilon\_{0,\mathbf{M}}(\omega) - |\mathbf{Q}|^2}.
$$

The *β*<sup>0</sup> and *β*<sup>M</sup> determine the character of the electromagnetic modes at the medium/dielectric substrate interface. To simplify the interpretation, we assume that the dielectric medium is vacuum, i.e., *ǫ*<sup>0</sup> = 1, and that the dielectric function of the substrate is constant *ǫ*M(*ω*) = *ǫ*M, which is a plausible approximation for many wide band gap insulators. In the vacuum, for *ω* > *Qc*, *β*<sup>0</sup> is a real number; therefore, the electromagnetic modes have a radiative character, and for *ω* < *Qc*, *β*<sup>0</sup> is an imaginary number so that the electromagnetic modes have evanescent character. The two regions are separated by the so-called 'light-line' *ω* = *Qc*, as illustrated by the magenta line in Figure 2a. In analogy to that, the two regions in the substrate are separated by the *ω* = *Qc*/ √ *ǫ*<sup>M</sup> line, as illustrated by the green line in Figure 2a. Since *ǫ*<sup>M</sup> > 1, the slope of the light-line in the substrate is smaller than in the vacuum, so in the gap *Qc*/ √ *ǫ*<sup>M</sup> < *ω* < *Qc*, the light propagates freely into the substrate but has an evanescent character in the vacuum region, as illustrated in Figure 2a. Therefore, the exciton–polariton mode *ω*ex-pol is expected to appear in the fully evanescent region *ω* < *Qc*/ √ *ǫ*M, since in that region it cannot be irradiated into the surrounding media. The evanescent character of the electric field produced by the exciton–polariton in the C<sup>60</sup> film is illustrated in Figure 2b.

**Figure 2.** (**a**) The character of the electromagnetic modes at the dielectric/vacuum (*ǫ*M/*ǫ*0) interface. In the region *ω* > *Qc*, the electromagnetic modes are entirely radiative (both in vacuum and in the dielectric), in the region *Qc*/ √ *ǫ*<sup>M</sup> < *ω* < *Qc*, they are radiative in the dielectric and evanescent in vacuum, and in the *ω* < *Qc*/ √ *ǫ*<sup>M</sup> region, they have fully evanescent character. In the latter region, the photon and molecular exciton (*ω*ex) hybridize, and an exciton–polariton (*ω*ex-pol) occurs. The measure of the coupling strength between the exciton *ω*ex and the photon is given by Rabi splitting Ω. (**b**) The evanescent electric field *Eµ*(*z*) produced by an exciton–polariton in the C<sup>60</sup> film.

We chose the electromagnetic modes to propagate in the **Q**<sup>0</sup> = **y** directions. In this case, the Dyson Equation (10) decouples into the independent matrix and the matrix tensor equations for the s(TE) and p(TM) polarizations, respectively. Here, we explore the s(TE)-polarized electromagnetic modes, which satisfy the matrix equation:

$$\mathcal{E}\_{\text{xx}}(|\mathbf{Q}|\mathbf{y},\omega,z\_{i\prime}z\_{j}) = \Gamma\_{\text{xx}}(|\mathbf{Q}|\mathbf{y},\omega,z\_{i\prime}z\_{j}) + \tag{14}$$

$$\sum\_{k} \Gamma\_{\text{xx}}(|\mathbf{Q}|\mathbf{y},\omega,z\_{i\prime}z\_{k}) \tilde{\sigma}\_{k}^{\text{xx}}(\omega) \mathcal{E}\_{\text{xx}}(|\mathbf{Q}|\mathbf{y},\omega,z\_{k\prime}z\_{j}),$$

where after using (11)–(13)

$$\Gamma\_{\rm xx}(|\mathbf{Q}|\mathbf{y},\omega,z\_i,z\_j) = -\frac{2\pi\omega}{\beta\_0 c^2} \left[ e^{i\beta\_0|z\_i-z\_j|} + r\_s e^{i\beta\_0(z\_i+z\_j)} \right].\tag{15}$$

The first term in (15) represents the incident electromagnetic field, while the second term represents the one reflected at the vacuum/substrate boundary. In the electrostatic or non-retarded limit *<sup>c</sup>* → <sup>∞</sup> and

$$\lim\_{\omega \to \infty} \Gamma\_{\text{xx}}(|\mathbf{Q}| \mathbf{y}, \omega, z\_i, z\_j) = 0,\tag{16}$$

so that all the properties presented below are a direct consequence of the binding between the s(TE)-polarized light (photons) and the molecular excitons, and they vanish, i.e., do not exist in the electrostatic limit.

#### *2.2. Calculation of the Optical Conductivity of a Single Molecule*

First, we explain the calculation of the molecular conductivity *σ µν i* (*ω*) in a 'selfstanding molecule'(*z*<sup>0</sup> → <sup>∞</sup>), and then we extend that to the case when the molecule is close to a dielectric surface, i.e., when *z*<sup>0</sup> is finite, and chosen to be a characteristic vdW distance. The basic ingredients we need to calculate the molecular conductivity *σµν*(*ω*) are the molecular orbitals *φn*(**r**) and the energies *En*, which can be obtained by solving the Kohn–Sham equation self-consistently. We assume that the molecules are periodically repeating so that they form a simple cubic (sc) Bravais lattice with the unit cell *a* and volume Ωsc = *a* 3 . The unit cell *a* is chosen so that there is no intermolecular overlap. This allows the molecular states |*φn*i to be calculated at the <sup>Γ</sup> point only. It should be emphasized that the sc lattice and the unit cell *a* are not related to the previously described FCC lattice with the unit cell '*a*3D'. The purpose of the sc lattice is only to determine the molecular states |*φn*i at the <sup>Γ</sup> point using the plane-wave DFT code. From now on, the conductivity of the 3D molecular crystal will be denoted as *σ* 3D *µν* (*ω*).

The nonlocal optical conductivity tensor in the 3D molecular crystal is [62,64–66]

$$\sigma\_{\mu\nu}^{\text{3D}}(\mathbf{r}, \mathbf{r}', \omega) = \frac{2i}{\omega} \sum\_{nm} \sum\_{n'm'} \mathcal{K}\_{n \to n'}^{m \leftarrow m'}(\omega) j\_{nm}^{\mu}(\mathbf{r}) [j\_{n'm'}^{\upsilon}(\mathbf{r}')]^\*,\tag{17}$$

where

$$j\_{nm}^{a}(\mathbf{r}) = \frac{e\hbar}{2im} \left\{ \phi\_n^\*(\mathbf{r}) \partial\_a \phi\_m(\mathbf{r}) - [\partial\_a \phi\_n^\*(\mathbf{r})] \phi\_m(\mathbf{r}) \right\} \tag{18}$$

represents the current produced by the transition between the molecular states |*φn*i → |*φm*i. Considering that the Bloch wave functions at the Γ point *φ<sup>n</sup>* are periodic functions, tensor (17) can be expanded in the Fourier series

$$
\sigma\_{\mu\nu}^{\text{3D}}(\mathbf{r}, \mathbf{r}', \omega) = \frac{1}{\Omega\_{\text{sc}}} \sum\_{\mathbf{G} \mathbf{G}'} e^{i \mathbf{G} \mathbf{r}} e^{-i \mathbf{G}' \mathbf{r}'} \sigma\_{\mathbf{G} \mathbf{G}'}^{\mu \nu}(\omega), \tag{19}
$$

where the Fourier coefficients are

$$
\sigma\_{\mathbf{G}\mathbf{G}'}^{\mu\nu,\mathfrak{3D}}(\omega) = \frac{i}{\omega} \frac{2}{\Omega\_{\mathbf{sc}}} \sum\_{nm} \sum\_{n'm'} j\_{nm}^{\mu}(\mathbf{G}) \mathcal{K}\_{n \to n'}^{m \leftarrow m'}(\omega) [j\_{n'm'}^{\upsilon}(\mathbf{G}')]^\*,\tag{20}
$$

and the current vertices are

$$j\_{nm}^{\mathfrak{a}}(\mathbf{G}) \, \, = \int\_{\Omega\_{\mathbf{sc}}} d\mathbf{r} e^{-i\mathbf{G}\mathbf{r}} \, j\_{nm}^{\mathfrak{a}}(\mathbf{r}). \tag{21}$$

The four-point polarizability K can be obtained by solving the Bethe–Salpeter (BS) equation [28,47]

$$\mathcal{K}\_{\mathfrak{n}\to\mathfrak{n}'}^{\mathfrak{m}\leftarrow\mathfrak{m}'}(\omega) = \mathcal{L}\_{\mathfrak{n}\to\mathfrak{n}'}^{\mathfrak{m}\leftarrow\mathfrak{m}'}(\omega) \, + \sum\_{\mathfrak{n}\_1\mathfrak{m}\_1} \sum\_{\mathfrak{n}\_2\mathfrak{m}\_2} \mathcal{L}\_{\mathfrak{n}\to\mathfrak{n}\_1}^{\mathfrak{m}\leftarrow\mathfrak{m}\_1}(\omega) \, \Xi\_{\mathfrak{n}\_1\to\mathfrak{n}\_2}^{\mathfrak{m}\_1\leftarrow\mathfrak{m}\_2} \mathcal{L}\_{\mathfrak{n}\_2\to\mathfrak{n}'}^{\mathfrak{m}\leftarrow\mathfrak{m}'}(\omega),\tag{22}$$

where the time-ordered electron–hole propagator is defined as

$$\mathcal{L}\_{n \to n'}^{m \gets m'} (\omega) \;= \delta\_{nn'} \delta\_{mm'} \int\_{-\infty}^{\infty} \frac{d\omega'}{2\pi i} \mathcal{G}\_{\mathbb{n}}(\omega') \mathcal{G}\_{\mathbb{m}}(\omega + \omega'). \tag{23}$$

The propagator (or Green's function) of an electron or a hole in a molecular state |*φn*i is

$$G\_{\rm n}(\omega) = \frac{1}{\omega - E\_{\rm n} + E\_{\rm n}^{\rm XC} - \Sigma\_{\rm n}^{\rm X} - \Sigma\_{\rm n}^{\rm C,0}(\omega)},\tag{24}$$

where the exchange of self-energy is

$$
\Sigma\_n^{\chi} = -\frac{1}{\Omega\_{\rm sc}} \sum\_m f\_m \sum\_{\mathbf{G} \mathbf{G}'} \rho\_{nm}^\*(\mathbf{G}) V\_{\mathbf{G} \mathbf{G}'}^{\mathbf{C}}(\mathbf{Q}) \rho\_{nm}(\mathbf{G}'). \tag{25}
$$

The correlation self-energy in the G0W<sup>0</sup> approximation is

$$\Gamma\_{\rm m}^{\rm C,0}(\omega) = \frac{1}{\Omega\_{\rm sc}} \sum\_{\rm m} \sum\_{\rm G\mathcal{G}'} \rho\_{\rm nm}^\*(\mathbf{G}) \rho\_{\rm nm}(\mathbf{G}') \left\{ (1 - f\_{\rm m}) \, \Gamma\_{\rm G\mathcal{G}'}^0(\omega - \mathbf{E}\_{\rm m}) - f\_{\rm m} \, \Gamma\_{\rm G\mathcal{G}'}^0(\mathbf{E}\_{\rm m} - \omega) \right\}, \tag{26}$$

where the correlation propagator is defined as

$$
\Gamma^{0}\_{\mathbf{G}\mathbf{G}'} (\omega) = \int\_0^\infty d\omega' \frac{S^0\_{\mathbf{G}\mathbf{G}'} (\omega')}{\omega - \omega' + i\delta'} \tag{27}
$$

and *f<sup>n</sup>* = *θ*(*E<sup>F</sup>* − *En*) is the Fermi–Dirac distribution at *T* = 0. The spectrum of the electronic excitation in a self-standing molecule is

$$S^{0}\_{\mathbf{G}\mathbf{G}'}(\omega) = -\frac{1}{\pi}Im\mathcal{W}^{0}\_{\mathbf{G}\mathbf{G}'}(\omega). \tag{28}$$

To avoid double-counting, we excluded the DFT exchange–correlation contribution *E* XC *n* from the KS energy *E<sup>n</sup>* in Equation (24). In the quasi-particle (QP) approximation, the electrons and holes have energies *E* QP *<sup>n</sup>* , which are the real poles of Green's function (24), i.e., they satisfy the equation

$$E\_n - E\_n^{\chi\mathcal{C}} + \Sigma\_n^{\chi} + \text{Re}\Sigma\_n^{\mathcal{C},0}(E\_n^{\mathcal{Q}\mathcal{P}}) = E\_n^{\mathcal{Q}\mathcal{P}}.\tag{29}$$

The electron/hole Green functions can now be approximated as

$$G\_n^{\rm QP}(\omega) = \frac{1 - f\_n^{\rm QP}}{\omega - E\_n^{\rm QP} + i\delta} + \frac{f\_n^{\rm QP}}{\omega - E\_n^{\rm QP} - i\delta},\tag{30}$$

and after they are used in (23), the time-ordered electron–hole propagator becomes

$$\mathcal{L}\_{n \to n'}^{m \gets m'} (\omega) \ = \frac{f\_n^{\text{QP}} - f\_m^{\text{QP}}}{\omega + E\_n^{\text{QP}} - E\_m^{\text{QP}} + i \delta s n (E\_m^{\text{QP}} - E\_n^{\text{QP}})} \delta\_{nm'} \delta\_{mm'} \tag{31}$$

where *δ* = 0 <sup>+</sup>. The 'time-ordered' screened Coulomb interaction is the solution of Dyson's equation

$$\mathcal{W}^{0}\_{\mathbf{G}\mathbf{G}'}(\omega) = V^{\mathbb{C}}\_{\mathbf{G}\mathbf{G}'} + \sum\_{\mathbf{G}\_1\mathbf{G}\_2} V^{\mathbb{C}}\_{\mathbf{G}\mathbf{G}\_1} \chi^{0}\_{\mathbf{G}\_1\mathbf{G}\_2}(\omega) \mathcal{W}^{0}\_{\mathbf{G}\_2\mathbf{G}'}(\omega),\tag{32}$$

where the matrix of the 'time-ordered' irreducible polarizability is

$$\chi^{0}\_{\mathbf{G}\mathbf{G}'}(\omega) = \frac{2}{\Omega\_{\mathbf{sc}}} \sum\_{nm} \frac{(f\_n - f\_m)\rho\_{nm}(\mathbf{G})\,\rho^\*\_{nm}(\mathbf{G}')}{\hbar\omega + E\_n - E\_m + i\delta\text{sgn}(E\_m - E\_n)} \tag{33}$$

and the charge vertices are

$$\rho\_{\rm mm}(\mathbf{G}) = \int\_{\Omega\_{\rm sc}} d\mathbf{r} \,\phi\_{n}^{\*}(\mathbf{r})e^{-i\mathbf{G}\mathbf{r}}\phi\_{m}(\mathbf{r}).\tag{34}$$

Since we calculate the optical conductivity of a single isolated benzene molecule, we have to exclude the effect on its polarizability due to the interaction with the surrounding molecules in the sc lattice. This is accomplished by using the truncated Coulomb interaction [67]

$$V\_{\mathbb{C}}(\mathbf{r} - \mathbf{r}') = \frac{\Theta(|\mathbf{r} - \mathbf{r}'| - R\_{\mathbb{C}})}{|\mathbf{r} - \mathbf{r}'|},\tag{35}$$

where Θ is the Heaviside step function, and *R*<sup>C</sup> is the range of the Coulomb interactions, i.e., the radial cutoff. The Coulomb interaction matrix to be used in (32) is then

$$V\_{\mathbf{G}\mathbf{G}'}^{\mathbb{C}} = \frac{1}{\Omega} \int\_{\Omega} d\mathbf{r} \int\_{\Omega} d\mathbf{r}' e^{-i\mathbf{G}\mathbf{r}} V\_{\mathbf{C}}(\mathbf{r}, \mathbf{r}') e^{i\mathbf{G}'\mathbf{r}'} = \frac{4\pi}{|\mathbf{G}|^{2}} [1 - \cos|\mathbf{G}|\mathcal{R}\_{\mathbf{C}}| \delta\_{\mathbf{G}\mathbf{G}'} . \tag{36}$$

The Bethe–Salpeter kernel is

$$
\Xi\_{n \to n'}^{m \leftarrow m'} = \Xi\_{n \to n'}^{H, m \leftarrow m'} - \frac{1}{2} \Xi\_{n \to n'}^{F, m \leftarrow m'} \tag{37}
$$

where the BS–Hartree kernel is

$$\Sigma\_{n \to n'}^{\text{H,m} \leftarrow m'} = \frac{1}{\Omega\_{\text{sc}}} \sum\_{\mathbf{G}\_1 \mathbf{G}\_2} \rho\_{nm}^\*(\mathbf{G}\_1) \ V\_{\mathbf{G}\_1 \mathbf{G}\_2}^{\mathbb{C}} \rho\_{n'm'}(\mathbf{G}\_2),\tag{38}$$

and the BS–Fock kernel is

$$\Xi\_{n \to n'}^{\mathbf{F}, m \leftarrow m'} = \frac{1}{\Omega\_{\text{sc}}} \sum\_{\mathbf{G}\_1 \mathbf{G}\_2} \rho\_{m'}^\*(\mathbf{G}\_1) \mathcal{W}\_{\mathbf{G}\_1 \mathbf{G}\_2}^0(\omega = 0) \rho\_{mm'}(\mathbf{G}\_2). \tag{39}$$

Here, the index '0' in *W*, *S*, Γ, and Σ <sup>C</sup> emphasizes that we consider the screened interaction in a self-standing molecule. Finally, considering that the interaction *V*<sup>C</sup> prevents the correlations between the conductivities in the adjacent cells, the conductivity of an isolated molecule is equal to the conductivity per unit cell

$$
\sigma\_{\mu\nu}(\omega) = \int\_{\Omega\_{\text{sc}}} d\mathbf{r} \int\_{\Omega\_{\text{sc}}} d\mathbf{r}' \sigma\_{\mu\nu}^{\text{3D}}(\mathbf{r}, \mathbf{r}', \omega). \tag{40}
$$

After using (19) in (40), the optical conductivity of a single molecule becomes

$$
\sigma\_{\mu\nu}(\omega) = \Omega\_{\text{sc}} \sigma\_{\mathbf{G}=0\mathbf{G}'=0}^{\mu\nu,\text{3D}}(\omega). \tag{41}
$$

After combining Equations (9), (20), and (41), we determine the explicit expression for the surface optical conductivity

$$\tilde{\sigma}\_{\mu\nu}(\omega) = \frac{2i}{\omega \,\mathrm{S}\_{\mathrm{fcc}}} \sum\_{nm} \sum\_{n'm'} j\_{nm}^{\mu}(\mathbf{G}=\mathbf{0}) \mathcal{K}\_{n\to n'}^{m\leftarrow m'}(\omega) [j\_{n'm'}^{\nu}(\mathbf{G}'=\mathbf{0})]^{\*},\tag{42}$$

which can be used in Dyson's equation for the electric field propagator (10). It is important to note that the dimension of the conductivity (42) is exactly the quantum of conductance *G*<sup>0</sup> = <sup>2</sup>*π<sup>e</sup>* 2 *h* , already the standardized unit for describing the optical conductivity in 2D crystals [28,47,65]. Accordingly, the *σ*˜*µν*(*ω*) represents the optical conductivity of one (e.g., *i*-th) molecular layer.

#### *2.3. Optical Conductivity in a Molecule Physisorbed at a Dielectric Surface*

We assume that a fullerene molecule, centered at *z* = *z<sup>i</sup>* , is physisorbed at the supporting crystal occupying the *z* < 0 half-space, as illustrated in Figure 3.

**Figure 3.** The fullerene molecule C<sup>60</sup> centered at *z* = *z<sup>i</sup>* is physisorbed at the supporting crystal occupying the half-space *z* < 0, with the dielectric properties approximated by the macroscopic dielectric function *ǫ*M(*ω*).

We further assume that the bonding between the molecule and the supporting crystal has the vdW character, which implies a small electronic overlap between the molecule and the substrate and, therefore, a small impact on the short-range electronic correlations to the molecular QP and optical properties. More precisely, the processes involving the electronic hopping between the molecule and the supporting crystal are neglected. However, we shall retain the processes of scattering an electron or a hole or an excited electron–hole pair by the potential induced at the crystal surface, ∆*V*.

If we have two valence electrons at points **r** and **r** ′ in a self-standing molecule, they interact via the bare Coulomb truncated potential *V*<sup>C</sup> (35), but they also polarize the molecule itself so that the total interaction between them is given by the potential *W*<sup>0</sup> , obtained as the solution of Equation (32). After the polarizable surface is brought close to the molecule, the electrons at **r** and **r** ′ polarize the surface as well so that the interaction between them, neglecting the polarization of the molecule, becomes

$$V\_{\mathbb{S}}(\mathbf{r}, \mathbf{r}', \omega) \;=\; V\_{\mathbb{C}}(\mathbf{r}, \mathbf{r}') + \Delta V\_{\mathbb{S}}(\mathbf{r}, \mathbf{r}', \omega), \tag{43}$$

where ∆*V*S(**r**,**r** ′ , *ω*) is the induced dynamic Coulomb potential coming from the excitations of the electronic modes or phonons at the surface. The total interaction between the electrons (including the polarization of the molecule) is then the solution of Dyson's equation

$$\boldsymbol{W\_{\mathbf{G}\mathbf{G}'}}(\omega) = \boldsymbol{V\_{\mathbf{G}\mathbf{G}'}}(\omega) + \sum\_{\mathbf{G\_1}\mathbf{G\_2}} \boldsymbol{V\_{\mathbf{G\_1}}^{\mathbf{S}}}(\omega) \boldsymbol{\chi\_{\mathbf{G\_1}\mathbf{G\_2}}^{0}}(\omega) \boldsymbol{W\_{\mathbf{G\_2}\mathbf{G}'}^{\mathbf{S}}}(\omega) \tag{44}$$

where

$$
\Delta V\_{\mathbf{G}\mathbf{G}'}^{\mathbf{S}}(\omega) = \frac{1}{\Omega\_{\mathbf{sc}}} \int\_{\Omega\_{\mathbf{sc}}} d\mathbf{r} \int\_{\Omega\_{\mathbf{sc}}} d\mathbf{r}' e^{-i\mathbf{G}\mathbf{r}} \Delta V\_{\mathbf{S}}(\mathbf{r}, \mathbf{r}', \omega) e^{i\mathbf{G}'\mathbf{r}'}.\tag{45}
$$

Here, the integration is constrained within the unit cell Ωsc centered at **r**<sup>0</sup> = (0, 0, *zi*) (to avoid interaction with the neighboring molecules, via ∆*V*S), as shown in Figure 3. The induced interaction in the region *z* > 0 can be written as [62]

$$
\Delta V(\mathbf{r}, \mathbf{r}', \omega) = \int \frac{d\mathbf{Q}}{(2\pi)^2} v\_{\mathbf{Q}} D(\mathbf{Q}, \omega) e^{i\mathbf{Q}(\mathbf{P} - \mathbf{P}')} e^{-Q(z + z')} \,\tag{46}
$$

where *v<sup>Q</sup>* = 2*π* |**Q**| . Since the supporting crystal dielectric response is approximated by the local 3D macroscopic dielectric function *ǫ*M(*ω*), the surface excitation propagator can be approximated as [68]

$$D(\mathbf{Q}, \omega) \approx D(\omega) = \frac{1 - \epsilon\_{\mathbf{M}}(\omega)}{1 + \epsilon\_{\mathbf{M}}(\omega)}.\tag{47}$$

After (46) and (47) are used in (45), we have

$$
\Delta V\_{\mathbf{G}\mathbf{G}'}^{\mathbf{S}}(\omega) = \frac{D(\omega)}{\Omega\_{\text{sc}}} \int \frac{d\mathbf{Q}}{(2\pi)^2} v\_{\mathbf{Q}} e^{-2Qz\_i} \mathbf{F}\_{\mathbf{G}}(\mathbf{Q}) \mathbf{F}\_{\mathbf{G}'}^\*(\mathbf{Q}), \tag{48}
$$

where the form factors are defined as

$$F\_{\mathbf{G}}(\mathbf{Q}) = 8(-1)^{\eta\_z} \frac{\sin[(Q\_x - G\_x)\frac{a}{2}] \sin[(Q\_y - G\_y)\frac{a}{2}] \sinh[Q\_z^a]}{(Q\_x - G\_x)(Q\_y - G\_y)(Q + iG\_z)}.$$

The reciprocal vectors of the sc lattice are **G** = (*Gx*, *Gy*, *Gz*), where *G<sup>x</sup>* = 2*πn<sup>x</sup> a* , *G<sup>y</sup>* = 2*πn<sup>y</sup> a* , *G<sup>z</sup>* = 2*πn<sup>z</sup> a* and *<sup>n</sup>x*, *<sup>n</sup>y*, *<sup>n</sup><sup>z</sup>* <sup>∈</sup> <sup>Z</sup>. After we determine the 'bare' potential (which includes the polarization of the surface) *V*<sup>S</sup> and the 'total' potential (which includes the polarization of the surface and of the molecule) *W*S, the calculation of the QP and optical properties of a molecule near the dielectric surface is equal to the procedure described in Section 2.2, except for that in the BS–Hartree kernel (38), we need to replace

$$V^{\mathbb{C}}\_{\mathbf{G}\_1\mathbf{G}\_2} \to V^{\mathbb{S}}\_{\mathbf{G}\_1\mathbf{G}\_2}(\omega)\_{\prime\prime}$$

in the BS–Fock kernel (39)

$$\mathcal{W}^{0}\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}(\omega=0) \to \mathcal{W}^{\mathbf{S}}\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}(\omega=0) \prime$$

and the spectrum of the electronic excitations (28) used to calculate the correlation selfenergy (26) is

$$S^{0}\_{\mathbf{G}\mathbf{G}'}(\omega) \to S^{\mathbf{S}}\_{\mathbf{G}\mathbf{G}'}(\omega) = -\frac{1}{\pi}Im\mathcal{W}^{\mathbf{S}}\_{\mathbf{G}\mathbf{G}'}(\omega). \tag{49}$$

#### *2.4. Computational Details*

The fullerene KS orbitals *φn*(**r**) and energies *E<sup>n</sup>* were calculated using the planewave self-consistent field DFT code (PWSCF) within the QUANTUM ESPRESSO (QE) package [69]. The core-electron interaction was approximated by the norm-conserving pseudopotentials [70,71] so that the number of occupied valance states was 120. The exchange-correlation (XC) potential was approximated by the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) functional [72]. The plane-wave cut-off energy was 60 Ry. The molecules were arranged in the simple cubic Bravais lattice of the unit cell *a* = 18 Å with one molecule per unit cell. Since there was no intermolecular overlap, the ground state electronic density was calculated at the Γ point only. The geometries were fully relaxed, with all forces .0.02 eV/Å. The RPA 'time-ordered' screened Coulomb interactions *<sup>W</sup>*0,S (Equations (32) and (44)) were calculated using the energy cut-off 2 Ry (∼27 eV), and the band summations '(*n*, *m*)' in the irreducible polarizability (33) were performed over 240 molecular states. The exchange self-energy (25) was calculated using the energy cut-off 8 Ry (∼109 eV) and the correlation self-energy (26) was determined (according to the cut-off in *W*0) using the energy cut-off 2 Ry (∼27 eV); the band summation '*m*' was performed using 240 molecular orbitals. The BS–Hartree kernel (38) and the 'bare' BS–Fock kernel (the Equation (39), derived using the bare interaction *V*C), were calculated using the energy cutoff 8 Ry (∼109 eV); the induced Fock kernel (the Equation (39), derived using the induced interaction *W*0,S − *V*C), was calculated using the energy cut-off 2 Ry (∼27 eV). During the evaluation of the BSE–HF kernels, we used 42 occupied (HOMO − 41, . . . , HOMO) and 42 unoccupied (LUMO, . . . , LUMO + 41) molecular states, i.e., the dimensions of the BSE–HF kernel matrix was 2 × 42 × 42 = 3528. To achieve the accurate (experimental) exciton energy, the calculation was performed beyond the Tamm–Dancoff approximation. The damping parameters *δ* used in (31) and in (33) were 50 meV and 200 meV, respectively. We

assume that the dielectric medium was vacuum (i.e., *ǫ*<sup>0</sup> = 1), and that the substrate was the aluminum-oxide Al2O3, described by the macroscopic dielectric function

$$\epsilon\_{\mathbf{M}}(\omega) = 1/\epsilon\_{\mathbf{G}=0\mathbf{G}'=0}^{-1}(\mathbf{q} \approx 0, \omega),\tag{50}$$

where the dielectric matrix is *<sup>ǫ</sup>*<sup>ˆ</sup> <sup>=</sup> ˆI <sup>−</sup> *<sup>V</sup>*<sup>ˆ</sup> *<sup>χ</sup>*<sup>ˆ</sup> 0 . The irreducible polarizability *χ* 0 is determined using Equation (33) for <sup>Ω</sup>sc → <sup>Ω</sup>, *<sup>n</sup>* → (*n*, **<sup>k</sup>**) and *<sup>m</sup>* → (*m*, **<sup>k</sup>** + **<sup>q</sup>**). Here **<sup>k</sup>**, **<sup>q</sup>**, and **G** are the 3D wave vector, the transfer wave vector, and the reciprocal lattice vector, respectively, corresponding to the bulk Al2O<sup>3</sup> crystal. The bare Coulomb interaction is *<sup>V</sup>***GG**′(**q**) = <sup>4</sup>*<sup>π</sup>* |**q+G**| <sup>2</sup> *δ***GG**′ . The ground state electronic density of the bulk Al2O<sup>3</sup> is calculated using 9 × 9 × 3 K-mesh, the plane-wave cut-off energy is 50 Ry, and the Bravais lattices are hexagonal (12 Al and 18 O atoms in the unit cell) with the lattice constants *a*Al2O<sup>3</sup> = 4.76 Å and *c*Al2O<sup>3</sup> = 12.99 Å. The Al2O<sup>3</sup> irreducible polarizability *χ* 0 is calculated using the 21 × 21 × 7 *k*-point mesh and the band summations (*n*, *m*) are performed over 120 bands. The damping parameter is *δ* = 100 meV and the temperature is *T* = 10 meV. For the optically small wave vectors **q** ≈ 0, the crystal local field effects are negligible, i.e., the crystal local field effects cut-off energy is set to zero. Using this approach, the Re*ǫ*<sup>M</sup> is almost constant (Re*ǫ*<sup>M</sup> <sup>≈</sup> 3) for low frequencies (*<sup>ω</sup>* <sup>&</sup>lt; 3 eV), i.e., in the IR and visible range, while Im*ǫ*<sup>M</sup> is equal to zero up to the band gap energy (*Eg*∼6 eV). Therefore, Al2O<sup>3</sup> is a good choice for the substrate in the visible and near-UV frequency range, since its electronic excitations are above that range, and its IR active SO phonons (at *ω*SO < 100 meV) [73] are far below the C<sup>60</sup> excitons, which means that in the frequency range of interest, there is no dissipation of the electromagnetic energy in the substrate (it is transparent). In addition to that, the dielectric function is mostly constant, but in this calculation, we used the fully dynamical and complex *ǫ*M(*ω*). The integration in (48) was performed over the two-dimensional wave vectors **Q** = (*Qx*, *Qx*) using a 121 × 121 rectangular mesh and the cut-off wave vector *Q*<sup>C</sup> = 0.5 a.u. For the radial cut-off in the truncated interaction (36), we used *R*<sup>C</sup> = *a*/2 = 9 Å.

#### **3. Results and Discussion**

Strong exciton–photon hybridization usually occurs due to the large overlap between excitons of large oscillatory strengths and confined photons. This is experimentally achieved by placing nanostructures or molecules of large oscillatory strength in a metallic cavity. Our goal here was to explore the electromagnetic eigenmodes, which are a mixture of 'free' s(TE) photons and excitons, i.e., a 'free' traveling photon (not confined between metallic walls) captured by an exciton and 'trapped' in a molecular nanostructure. This concept makes it possible to clearly estimate the contribution of the 'free' photon in the hybrid exciton–photon mode.

First, we determined the QP and optical properties of a self-standing layer of the C<sup>60</sup> molecules, which corresponded with the properties of the gas-phase molecule, since at the *G*0*W*0-BSE stage of the calculations, the molecules were not coupled. The calculated G0W<sup>0</sup> HOMO-LUMO gap was *E<sup>g</sup>* = 4.66 eV, which is in good agreement with the experimental value 5 eV [74–77]. When a molecular layer is deposited on an Al2O<sup>3</sup> dielectric surface, the band gap reduces to *E<sup>g</sup>* = 4.26 eV. Here, the separation was chosen to be *z*<sup>0</sup> = 6.5 Å which corresponds to the characteristic vdW atom–atom separation of 3 Å (Note that *z*<sup>0</sup> is defined as the distance between the substrate and the center of the molecule, as denoted in Figure 1, so the C<sup>60</sup> molecule radius (3.57 Å) has to be added to the vdW atom–atom distance). For comparison, we determined the HOMO-LUMO gap for the molecule at a silver (Ag) surface also described in terms of the *ab initio* macroscopic dielectric function. For the same separation (*z*<sup>0</sup> = 6.5 Å), the gap was *E<sup>g</sup>* = 3.81 eV; compared with the insulator surface, the reduction is twice as strong. The image theory estimation of the HOMO-LUMO gap at the metallic surface is *Eg*∼3.56 eV which, as expected, overestimates the G0W<sup>0</sup> result. Figure 4a shows the calculated G0W0-BSE optical conductivities *σ*˜*x*(*ω*) in the selfstanding molecular layer (black solid) and in the molecular layer deposited at the Al2O<sup>3</sup>

dielectric surface (cyan dashed). For comparison, the experimental optical absorption in the fullerene C<sup>60</sup> [78] is presented by the red circles. This experimental study, as well as some others [79,80], show three broad excitation bands at *ωex*1∼ 3.9 eV, *ωex*2∼4.9 eV and *ωex*3∼6 eV, which is in very good agreement with our results. It should be emphasized that, since the band gap is *E<sup>g</sup>* = 4.66 eV, only the excitation band *ωex*<sup>1</sup> can be considered as an exciton by definition, and its binding energy is *E<sup>g</sup>* − *h*¯ *ωex*1∼0.76 eV. In the literature, the strong maximum *ωex*<sup>3</sup> is usually referred to as the *π* plasmon [81,82]. All three excitons are the result of the electronic transitions within the C<sup>60</sup> *π*-complex [83]. When the molecule is deposited at the Al2O<sup>3</sup> surface, the excitation band barely changes at all, due to the well-known cancellation effect: the substrate weakens the interaction between the excited electron and the hole, which reduces the exciton binding energy and, therefore, cancels the gap reduction. This phenomenon was studied in detail in references [62,66,84]. Since the influence of the dielectric surface on the molecular optical conductivity is weak, we shall further approximate *σ*˜ *x i* (*ω*) ≈ *σ*˜ *x* (*ω*), where *σ*˜ *x* (*ω*) is the optical conductivity in the self-standing molecular single layer, i.e., for *<sup>z</sup>*<sup>0</sup> → <sup>∞</sup>.

**Figure 4.** (**a**) The optical conductivities *σ*˜*xx*(*ω*) in the C<sup>60</sup> single layer in vacuum (black solid), in the C<sup>60</sup> single layer at the Al2O<sup>3</sup> dielectric surface (cyan dashed) (where *z*<sup>0</sup> = 6.5 Å), and the experimental optical absorption in the gas phase fullerene C<sup>60</sup> (red circles). (**b**) The lower and upper exciton–polariton branches, LPB and UPB, respectively, in the self-standing C<sup>60</sup> film (blue dots) and the C<sup>60</sup> film deposited at the Al2O<sup>3</sup> surface (red dots). The LPB corresponds to the dispersion relation of the exciton–polariton *ω*ex1-pol(*Q*) appearing in Equation (53). The number of C<sup>60</sup> layers is *N* = 6. (**c**) The spectra of the s(TE)-polarized electromagnetic modes *S*(*Qex*<sup>1</sup> , *ω*) in the self-standing C<sup>60</sup> films for *Q*<sup>0</sup> ex1 <sup>=</sup> 0.02 nm−<sup>1</sup> (blue solid) and in the C<sup>60</sup> films at the Al2O<sup>3</sup> surface for *Q*ex1 = 0.035 nm−<sup>1</sup> (red dashed). The number of the single layers is *N* = 0, 1, 2, 3, . . . , 10, where the case *N* = 0 corresponds to the spectrum of the free photons in vacuum or at the vacuum/Al2O<sup>3</sup> interface (the photon continuum). All spectra for the supported films are multiplied by factor 2.

This means that the impact of the dielectric substrate on the electromagnetic modes in the C<sup>60</sup> films is reduced to the propagator Γ sc, which appears in Dyson's Equation (14). The spectra of the s(TE)-polarized electromagnetic modes in the C<sup>60</sup> films will be analyzed using the real part of the propagator E in the topmost molecular layer

$$S(\mathbf{Q}\_{\prime}\boldsymbol{\omega}) = \operatorname{Re}\mathcal{E}\_{\text{xx}}(|\mathbf{Q}|\mathbf{y}\_{\prime}\boldsymbol{\omega}\_{\prime}\boldsymbol{z}\_{\text{N}}\boldsymbol{z}\_{\text{N}}).\tag{51}$$

As already mentioned, the spacing parameters were chosen to be *z*<sup>0</sup> = 6.5 Å and ∆ = 8.1 Å. In the near field spectroscopy experiment, the incident photon of the wavelength

$$
\lambda\_{\rm ex1} = \frac{2\pi c}{\omega\_{\rm ex1}\sqrt{\epsilon(\omega\_{\rm ex1})}} \tag{52}
$$

couples to polarization modes in some sub-wavelength nanostructure (e.g., AFM tip) and is scattered (diffracted) to all *λ* > *λ*ex1 and *λ* < *λ*ex1 so that it could excite all radiative electromagnetic modes, as well as all evanescent ones at a fixed frequency *ω*ex1. The Rabi splitting (the measure of the free-photon participation in the exciton–polariton mode) is defined as the difference between the exciton *ωex*<sup>1</sup> and the exciton–polariton *ω*ex1-pol

$$
\hbar \Omega = \hbar \omega\_{\rm ex1} - \hbar \omega\_{\rm ex1-pol}(Q\_{\rm ex1}) \tag{53}
$$

where *Qex*<sup>1</sup> = 2*π*/*λ*ex1 is the wave vector at which the exciton *ω*ex1 and the photon cross.

Figure 4b shows the dispersion relations of the exciton–polaritons *ω*ex1-pol(*Q*) in the self-standing C<sup>60</sup> film (blue dots), and in the C<sup>60</sup> film deposited at a Al2O<sup>3</sup> surface (red dots). The number of C<sup>60</sup> single layers in the film was chosen to be *N* = 6. The energies of the exciton–polaritons *ω*ex1-pol(*Q*) correspond to the maxima in the spectral function *S*(*Q*, *ω*), appearing below the exciton energy *ω*ex1 [see Figure 4c]. For the selfstanding film, *ǫ*(*ω*) = 1, so that the crossing wave vector, according to the Equation (52), was *Q*<sup>0</sup> ex1 = 0.02 nm−<sup>1</sup> , and the Rabi splitting, according to the Equation (53), was *h*¯ Ω<sup>0</sup> = 650 meV. The red dots in Figure 4b show that the presence of the substrate reduces the bending of the exciton–polariton dispersion relation significantly, therefore reducing the corresponding Rabi splitting. In this case, *ǫ*(*ω*) = *ǫ*M(*ω*), which gives the crossing wave vector *Q*<sup>S</sup> ex1 = 0.035 nm−<sup>1</sup> and, thus, the Rabi splitting ¯*h*Ω<sup>S</sup> = 103 meV.

Figure 4c shows the spectra of the s(TE)-polarized electromagnetic modes *S*(*Qex*1, *ω*) in the self-standing C<sup>60</sup> films for *Q*<sup>0</sup> ex1 = 0.02 nm−<sup>1</sup> (blue solid) and in the C<sup>60</sup> films deposited at the Al2O<sup>3</sup> surface for *Q<sup>S</sup>* ex1 = 0.035 nm−<sup>1</sup> (red dashed). The number of layers is set to be *N* = 0, 1, 2, 3, . . . , 10, where the case *N* = 0 corresponds to the spectra of the free photons in vacuum or at the vacuum/Al2O<sup>3</sup> interface (photons continuum). Due to the lower intensity, all spectra for the supported films were multiplied by factor 2. We can clearly see the exciton–polariton peaks separating from the photon continuum as the number of layers increases. In the self-standing films, the exciton–polariton is already well separated for *N* = 2, while in the supported films, it occurs for *N* = 4. For *N* = 10, the giant light–matter coupling was achieved providing large Rabi splitting Ω<sup>0</sup> = 1334 meV and Ω*<sup>S</sup>* = 670 meV in the self-standing and the supported films, respectively. Here we limited our study to the exciton–polaritons in the quasi-2D nanostructures so that the maximum number of the molecular layers was limited to *N* = 10, still satisfying the sub-wavelength limit *<sup>λ</sup>*ex1 <sup>&</sup>gt; (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>)∆. For a further increase of the number of layers *<sup>N</sup>*, the Rabi splitting continued to increase, but the experimental limitations in the realization of such a system raise the question of the plausibility of such a result.

We need to emphasize that the anti-crossing behavior and the Rabi splitting into the lower polariton band (LPB) and the upper polariton band (UPB) is a well-defined concept only when describing the interaction between the well-defined eigenmodes, such as excitons and cavity photons [33] or BSWs [48]. In our case, the only well-defined modes (boson) were the excitons *ω*ex1-*ω*ex3. On the other hand, the photons appeared as a continuum of eigenmodes above the light lines *Qc* or *Qc*/ √ *ǫ*M. This means that it only made sense to observe the deformation of the exciton dispersion *ω*ex1, characterized by the exciton bending parameter Ω. However, in Figure 5, it seems that the edge of the photon continuum *Qc* or *Qc*/ √ *ǫ*<sup>M</sup> behaves similar to a well-defined eigenmode, so in Figure 4b,c we also introduced and denoted the LPB and UPB, exactly as they appear in Figure 5b,e. The LPB corresponds to the exciton–polariton dispersion relation *ω*ex1-pol(*Q*) appearing in (53).

To present a more extensive picture of the electromagnetic modes in the C<sup>60</sup> films, including the hybridizations of the photon and the higher excitons *ω*ex2 and *ω*ex3, Figure 5a–c show the spectral intensity of the s(TE)-polarized electromagnetic modes *S*(*Q*, *ω*) in the self-standing C<sup>60</sup> films for *N* = 3, 6, and 9, respectively. Figure 5d–f show the same for the C<sup>60</sup> films deposited at the Al2O<sup>3</sup> surface. It is immediately obvious that the most intensive electromagnetic modes occur in the evanescent regions *ω* < *Qc* and *ω* < *Qc*/*ǫ*M(*ω*), for the self-standing and supported films, respectively. These modes produce an electric field, which spreads within the film and exponentially decays outside of the film [as illustrated in Figure 2b], which is one of the inherent properties of the 2D exciton–polaritons. Moreover, as the field is more confined, it is stronger, so that the polariton modes confined at a sub-wavelength scale are very interesting for many applications. The free photon continua can be seen as the low-intensity patterns in the regions *ω* > *Qc*, and *ω* > *Qc*/*ǫ*M(*ω*) in the self-standing and supported films, respectively, with their intensity weakening with the number of layers *N*. If the electromagnetic eigenmode, for example in the C<sup>60</sup> film at the Al2O3, would occur in the region *ω* > *Qc*/*ǫ*M(*ω*), it would be irradiated (decay radiatively) into the Al2O<sup>3</sup> crystal, so that these modes appear only as weak resonances. Such weak resonances can be seen as weak horizontal patterns (parallel with the excitons *ω*ex1-*ω*ex3) entering the radiative region.

**Figure 5.** The spectral intensity of the s(TE)-polarized electromagnetic modes *S*(*Q*, *ω*) in the selfstanding C<sup>60</sup> films for (**a**) *N* = 3, (**b**) *N* = 6, and (**c**) *N* = 9. Figures (**d**–**f**) show the same for the C<sup>60</sup> films deposited at the Al2O<sup>3</sup> surface. In both cases, the strong electromagnetic modes (*ω*ex1-pol, *ω*ex2-pol, and *ω*ex3-pol) occur in the evanescent regions *ω* < *Qc* and *ω* < *Qc*/*ǫ*M(*ω*).

In Figure 5a, we can see the beginning of the hybridization between the free photon and the three C<sup>60</sup> excitons (*ω*ex1-*ω*ex3), which visibly intersect the photon line *ω* = *Qc*. In Figure 5b,c, the photon line is already significantly deformed between the excitons and pushed far into the evanescent region, which means that the photon *Qc* and the three excitons *ω*ex1, *ω*ex2, and *ω*ex3 are strongly coupled and converted into three exciton– polaritons *ω*ex1-pol, *ω*ex2-pol, and *ω*ex3-pol. In Figure 5d–f, the exciton–polaritons are less

intense and pushed into the evanescent region *ω* < *Qc*/*ǫ*M(*ω*), the use of the substrate obviously weakens the intensity of the exciton–polaritons. However, in Figure 5e,f, we can see the formation of the three strong exciton–polariton modes *ω*ex1-pol–*ω*ex3-pol, for larger numbers of layers *N* = 6 and 9, respectively. Finally, this confirms that the strong binding between the transverse s(TE) photon and the molecular excitons is a realistic physical phenomenon that can be achieved experimentally. The use of a metal substrate is not recommended because a metallic surface will strongly quench the exciton–polariton modes. To support this argument, we performed the same calculation using a silver (Ag) substrate and obtained much weaker exciton–polariton modes.

Finally, we noticed that the hybridization between the exciton and the photon increases weakly if just the exciton oscillatory strength (e.g., S) increases, but strongly if the number of spatially separated molecular layers N increases. For example, the exciton–photon binding is much stronger in the N-separated layers of the oscillatory strength S than in one layer of the oscillatory strength N× S. This suggests that the photon-exciton coupling can be increased simply by increasing the number of layers N in a multilayered van der Waals heterostructure.

#### **4. Conclusions**

We showed that the 2D-layered heterostructures, consisting of a larger number of exciton-active single layers or 2D crystals (e.g., *N* > 5), can support evanescent s(TE) polarized exciton–polaritons with strong photonic character. We obtained giant Rabi splitting of more than Ω<sup>0</sup> = 1000 meV and Ω*<sup>S</sup>* = 500 meV in the self-standing and supported C<sup>60</sup> films for *N* = 10, respectively. This investigation has fundamental and practical contributions. We predict the existence of the evanescent s(TE) polarization modes with significant photonic character, which vanish exactly without the photon admixture (for *<sup>c</sup>* → <sup>∞</sup>). Unlike the well-known p(TM) polarization modes (e.g., the plasmon polariton for *<sup>c</sup>* → <sup>∞</sup> collapses into a plasmon), this is a new fundamental contribution. We also demonstrated that exciton–photon coupling can be manipulated simply by changing the number of single layers (*N*) in a vdW-layered heterostructure. Moreover, due to the fact that the vdW heterostructure with a thickness of just a few molecular (or atomic) layers supports the confined photons, it can be easily implemented in the photonic integrated circuits or photonics chips. The disadvantage of these types of photonic modes is that they do not couple directly to the free photons (external light). However, once excited, these modes can be easily manipulated (since they are trapped in the nanostructure). For example, by patterning the vdW nanoribbons on a dielectric wafer (patterning the photonics circuits), the direction of the photon propagation can easily be modified. Moreover, the exciton–polaritons can be easily switched (at the contact) into evanescent modes in another nanostructure. Finally, these layered vdW heterostructures can be applied in photonic devices, such as light sources (LED), telecommunications (as waveguides or optical cables), or chemical and biological sensing.

**Author Contributions:** Conceptualization, V.D.; methodology, V.D.; software, V.D. and L.M.; validation, V.D. and L.M.; formal analysis, V.D.; investigation, V.D. and L.M.; resources, V.D. and L.M.; data curation, V.D.; writing—original draft preparation, V.D.; writing—review and editing, V.D. and L.M.; visualization, V.D.; supervision, V.D.; project administration, V.D.; funding acquisition, V.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Croatian Science Foundation (grant no. IP-2020-02-5556) and the European Regional Development Fund for the "QuantiXLie Centre of Excellence" (grant KK.01.1.1.01.0004).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge support from the Donostia International Physic Center (DIPC) computing center, which provided the computational resources.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

