<figure><img src="image_1.png" /></figure>

HAL
open science

# Optical properties of quasi-two-dimensional objects from time-dependent density functional theory: Longitudinal versus transverse dielectric functions

Stefano Mazzei, Christine Giorgetti

► To cite this version:

Stefano Mazzei, Christine Giorgetti. Optical properties of quasi-two-dimensional objects from time-dependent density functional theory: Longitudinal versus transverse dielectric functions. Physical Review B, 2023, 107 (16), pp.165412. 10.1103/PhysRevB.107.165412. hal-04072130

HAL Id: hal-04072130

https://hal.science/hal-04072130

Submitted on 17 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Optical properties of quasi-two-dimensional objects from time-dependent density functional
theory: Longitudinal versus transverse dielectric functions

Stefano Mazzei and Christine Giorgetti

Laboratoire des Solides Irradiés, Ecole polytechnique, CNRS, CEA/DRF/IRAMIS, Institut Polytechnique de Paris, F-91128 Palaiseau, France
and European Theoretical Spectroscopy Facility (ETSF)

<figure><img src="image_3.png" /></figure>

(Received 26 September 2022; revised 1 March 2023; accepted 15 March 2023; published 17 April 2023)

Comprehension of the electronic properties of nano-objects is a key to defining dedicated properties, which can be adjusted by changing their size. Beyond confinement effects, the presence of interfaces, i.e., places where there is an abrupt change of electronic density, should also play a role. Time-dependent density functional theory (TD-DFT) is a state-of-the-art *ab initio* formalism in which this effect is accounted for through the so-called local field effects. In an earlier paper [S. Mazzei and C. Giorgetti, Phys. Rev. B 106, 035431 (2022)], we showed that the framework inherited from three-dimensional crystals could not provide reliable absorption spectra. In the present paper, we propose to calculate the macroscopic average of the dielectric tensor of a quasi-two-dimensional (2D) object from the response function of the density to the total macroscopic potential in order to avoid use of the so-called Adler and Wiser formula. We evidence that the inclusion of interfaces in the thickness of the slab causes the response function for the out-of-plane component to move sharply from the bulk absorption resonance to the plasmon one. This shows that the longitudinal-longitudinal contraction of the dielectric tensor is no longer equal to the transverse-transverse one in a quasi-2D object for out-of-plane perturbation. Nevertheless, we also show that the macroscopic average of the dielectric tensor of an ultrathin slab calculated within the longitudinal formalism of TD-DFT depicts the properties of the transverse reflectance and transmittance spectra of a thin slab.

DOI: 10.1103/PhysRevB.107.165412

# I. INTRODUCTION

Nanomaterials are expected to exhibit novel properties as a consequence of electronic confinement. In particular, varying their size or shape should enable us to define dedicated optical properties [1] due to the modification of the screening effects, leading, for example, to an adjustable band gap [2–5]. Understanding these properties at the atomic scale is key to achieving this goal, and *ab initio* theoretical spectroscopy [6] provides state-of-the-art frameworks in which it can be accomplished. In particular, time-dependent density functional theory (TD-DFT) [7–9] offers an efficient scheme based on the solution of the Dyson equation, which relates the density response functions to the total and external potentials. The kernel of this equation allows one to account for the exchange and correlation phenomena at the origin of the excitonic effects [10], and for the fluctuations at the atomic scale of the induced density, called the local field effects (LFEs), governed by the Hartree potential.

Powerful numerical codes have been developed for three-dimensional (3D) periodic crystals, exploiting the reciprocal space and plane-wave basis-set framework. To treat isolated nano-objects in this periodic formalism, one uses the supercell approach, where the object is embedded in vacuum. To avoid spurious interactions between replicas, one can extend the size of the supercell. It has been shown that the results behave like the effective medium theory [11–13], but this physical picture

is dubious since it conceives the properties of the isolated nano-object as an averaging with the surrounding vacuum. Moreover, it has also been shown that for the out-of-plane component in 2D geometry, the normalization of the response function with the size of the supercell prevents us from extracting a reliable macroscopic dielectric function [13]. Another approach, based on the use of a cutoff of the Coulomb interaction [14,15], has been proposed. Beyond the question of the vacuum still contained in the volume defined by the cutoff parameter, no out-of-plane component can be extracted for 2D material in the optical limit. Finally, a method called Selected-G was developed to overcome these drawbacks, with a modified expression for the Coulomb potential called the slab potential to isolate the replicas for quasi-2D objects, without any dependence on the size of the supercell [13,16]. It has been successfully applied to calculate the absorption spectrum as well as the second-harmonic generation in the thick slab limit in order to depict surface properties [13,17], and to describe screening effects for an in-plane perturbation in thin slabs [16].

The optical absorption spectrum is the imaginary part of the macroscopic average of the dielectric tensor. In Ref. [18], we solved the Dyson equation to obtain the response function to the external potential $\chi$. For 3D crystals, it is directly related to the macroscopic average of the inverse dielectric matrix $\epsilon_{00}^{-1}(\mathbf{q}; \omega) = 1 + \frac{4\pi}{|\mathbf{q}|^2} \chi_{00}(\mathbf{q}; \omega)$, which is expected to give the plasmon, and the use of the Adler and Wiser formula [19–21] $\epsilon_M(\omega) = \lim_{\mathbf{q}\to0} 1/\epsilon_{00}^{-1}(\mathbf{q}; \omega)$ allows for the calculation of the absorption spectrum in the optical limit. In our previous paper, we showed that the Adler and Wiser formula

<footnote>*christine.giorgetti@polytechnique.edu</footnote>
could not be applied for a quasi-2D object. For the in-plane component, we demonstrated that the macroscopic dielectric function is directly proportional to the average of the density response function (to the external potential). For the out-of-plane component, where the Adler and Wiser formula should in principle be valid, we showed that the difficulty in defining the thickness of matter and calculating a reliable value for the average of the induced electric field leads to a continuous shift of the absorption resonance. This behavior is similar to the effective-medium theory between matter and vacuum [11,13], but it occurs in the region that seems to define the interfaces, and the Adler-Wiser formula clearly gives wrong spectra [18,22]. Nevertheless, we also showed that the different in-plane and out-of-plane macroscopic dielectric functions calculated with the same value for the thickness of matter gave the same reflectance and transmittance spectra [18,22]. This means that the response function calculated within TD-DFT, which is a longitudinal formalism, still contains the electronic properties at the origin of the optical response, which is in essence a transverse phenomenon.

The purpose of this paper is to extract these optical properties from a calculation of the macroscopic average of the dielectric tensor that avoids the inversion of the quantity $\epsilon_{00}^{-1}(\mathbf{q}; \omega) = 1 + \frac{4\pi}{|\mathbf{q}|^2} \chi_{00}(\mathbf{q}; \omega)$, where $\chi_{00}$ is ambiguously normalized. We thus propose calculating the response function to the macroscopic part of the total potential, called $\bar{\chi}$, which for an infinite crystal leads directly to the absorption spectrum through $\epsilon_M(\mathbf{q}; \omega) = 1 - \frac{4\pi}{|\mathbf{q}|^2} \bar{\chi}_{00}(\mathbf{q}; \omega)$, where $\mathbf{q}$ is a vanishing reciprocal space vector allowing for the definition of the excitation direction. The direction of this momentum transfer mimics the direction of the electric field of the electromagnetic wave. The amplitude of $\epsilon_M(\mathbf{q}; \omega)$ does not depend on the choice of $|\mathbf{q}|$ in the optical limit. With such an expression, the ambiguity of the normalization of the response function would have only a small effect on the amplitude of the spectrum, but at least we expect the spectral shape to be correct. The difference between the macroscopic average of $\chi$ and $\bar{\chi}$ comes from the long-range component of the Coulomb potential, which is set to zero when solving the Dyson equation in reciprocal space for $\bar{\chi}$ as compared to $\chi$. This difference means that only microscopic components of the density response function exist when the excitation is due to a photon, contrary to what occurs when the perturbation is an electron. Nevertheless, beyond the scaling effect, the role of the interfaces will be shown to be crucial, and since they are characterized by an abrupt change of the electronic density, we will work within the random phase approximation (RPA) to focus on the LFEs. We will use the Selected-G formalism with the slab potential [13,16], as well as a mixed-space approach [13,18], where the in-plane directions are treated in reciprocal space, and the direction perpendicular to the slab is described in real space, in order to tackle the question of the definition of the interfaces. This leads us to provide a criterion to define the thickness of a quasi-2D object, and to question the relation between the longitudinal-longitudinal and transverse-transverse contractions of the dielectric tensor for the out-of-plane direction of a 2D object. We will still use a silicon slab as a model system, since the absorption and plasmon resonance are separated of more that 10 eV, which allows for a clear distinction [18,22].

But, we will show that the results are transferable to traditional 2D materials such as the stacking of a few layers of graphene.

This paper is organized as follows: In Sec. II, we calculate the macroscopic dielectric function of an isolated slab using the reciprocal space framework with the Selected-G method and the slab potential [13,16]. We evidence that a very small variation of the definition of the thickness of the slab leads to a change in the nature of the response function for the out-of-plane component. To clarify this effect, we adopt the mixed space framework in Sec. III, where we first establish how to suppress the long-range component of the Coulomb potential to obtain the response function to the total macroscopic potential, and we propose how to define the thickness of the matter. We also question the relation between the longitudinal-longitudinal and transverse-transverse contractions of the dielectric tensor for 2D objects. For this purpose, in Sec. IV, we analyze the optical response of a thin slab obtained by Airy's formula for reflectance and transmittance, and we show the connection with the macroscopic dielectric function calculated within TD-DFT. Finally, in Sec. V, we show that these findings are also valid for traditional 2D materials.

## II. A QUASI-2D OBJECT IN RECIPROCAL SPACE

To depict a surface or a 2D object within periodic boundary codes, one defines a supercell that contains the slab and vacuum in the perpendicular direction. It has been shown that, when the response function includes the local field effects, the spectra behave like an effective-medium theory with vacuum [13], which does not provide a correct physical picture. The Selected-G formalism [13] allows one to overcome the problem of the presence of the spurious vacuum for surfaces [13] as well as for 2D objects [16]. In this framework, one solves the Dyson equation on a set of reciprocal space vectors $\{\tilde{\mathbf{G}}\}$ defined according to the thickness of the matter $L_z^{\text{mat}}$, i.e., $\tilde{\mathbf{G}} = \mathbf{G}_{\parallel} + \tilde{\mathbf{G}}_z$ with $\tilde{\mathbf{G}}_z = n_z * 2\pi/L_z^{\text{mat}}$, $n_z \in \mathbb{Z}$ [13]. The remaining delicate point is to relate this density response function to the macroscopic dielectric function to obtain the absorption spectrum [18].

It can be shown that in reciprocal space, for a thin slab, the longitudinal-longitudinal contraction of the macroscopic dielectric tensor in the optical limit $\tilde{\epsilon}_M(\omega)$ can be calculated using the expression [23]

$$ \epsilon_M^{\text{LL}}(\mathbf{q}; \omega) = 1 - \frac{4\pi}{|\mathbf{q}|^2} \bar{\chi}_{00}(\mathbf{q}; \omega), \quad (1) $$

where $\bar{\chi}$ is the solution of a modified slab Dyson equation:

$$ \bar{\chi}_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}'}(\mathbf{q}; \omega) = \chi_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}'}^0(\mathbf{q}; \omega) + \sum_{\tilde{\mathbf{G}}_1 \tilde{\mathbf{G}}_2} \chi_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}_1}^0(\mathbf{q}; \omega) \bar{V}_{\tilde{\mathbf{G}}_1 \tilde{\mathbf{G}}_2}(\mathbf{q}) \bar{\chi}_{\tilde{\mathbf{G}}_2 \tilde{\mathbf{G}}'}(\mathbf{q}; \omega) \quad (2) $$

in which $\bar{V}_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}'}$ is derived from the slab potential $\tilde{V}_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}'}$ according to

$$ \bar{V}_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}'} = \begin{pmatrix} 0 & 0 \\ \tilde{V}_{\tilde{\mathbf{G}}0} & \tilde{V}_{\tilde{\mathbf{G}}\tilde{\mathbf{G}}'} \end{pmatrix} $$
with {$\tilde{G}$} the reciprocal space vectors defined according to the Selected-G formalism [13,16]. $\chi^0$ is the Kohn-Sham independent-particle response function. The slab potential diverges for $G_\parallel = 0$, i.e., in the optical limit. For this reason, a vanishing $q_\parallel$ vector must always be present in the calculation, prohibiting us from obtaining the $\varepsilon_M^{LL}(q_z; \omega)$ from a direct calculation. It must result instead from a linear combination of spectra calculated with ($q_\parallel$) and ($q_\parallel \pm q_z$) (see Appendix A).

To study the influence of the thickness of matter, we built a slab of 16 bilayers of silicon atoms corresponding to a stacking of four conventional cells of silicon ($a_{cell} = 10.263$ bohrs) leading to an “atomic” thickness of $L_z^{atom} = 41.052$ bohrs [24]. Then we embedded it in a supercell of height $L_z = 205.26$ bohrs, corresponding to five times the atomic thickness. To vary the $L_z^{mat}$ size, we solve the Dyson equation [Eq. (2)] within the Selected-G formalism, with different values of ratio $L_z/L_z^{mat}$ equal to 5, 4, 3, 2, and 1 defining, respectively, $L_z^{mat} = 41.052$ bohrs ($=L_z^{atom}$), 51.315 bohrs, 68.42 bohrs, 102.630 bohrs, and 205.26 bohrs. These different values for $L_z^{mat}$ depict a system with a constant thickness of matter defined according to the extension of atomic positions ($L_z^{atom} = 41.052$ bohrs, which we will called the “zero vacuum” case), plus the addition of a slice of vacuum of increasing thickness.

The imaginary part of the macroscopic dielectric function [Eq. (1)] for these different thicknesses is shown in Fig. 1. We neglected the in-plane LFE by setting $G_\parallel = 0$ in Eq. (2). The in-plane components are plotted as continuous and dashed (blue) lines. They are all composed of one peak located at around 4 eV, and their amplitude is scaled by a factor $1/L_z^{mat}$. The out-of-plane components are drawn as

<figure><img src="image_4.png" /><figcaption>FIG. 1. Im[$\varepsilon_M^{LL}(q_z; \omega)$] [Eq. (1)] for in-plane ($q_x$) (lines without symbols, blue) and out-of-plane excitations ($q_z$) (lines with symbols, brown and orange) for slabs of different thicknesses [$L_z^{mat} = 41.052$ bohrs ($=L_z^{atom}$) (continuous lines/diamond), 51.315 bohrs (small dashed lines/circle), 68.42 bohrs (large dashed lines/square), 102.630 bohrs (dotted-small dashed lines/up-triangle), and 205.26 bohrs (dotted-large dashed lines/down-triangle)], made of always the same quantity of silicon, and an increasing slice of vacuum. $G_\parallel = 0$ in [Eq. (2)]. Inset: Im[$\varepsilon_M^{LL}(q_z; \omega)$] for $L_z^{mat} = 41.052$ bohrs (diamond, brown), 44 bohrs (left-triangle, red), 46 bohrs (cross, turquoise), 48 bohrs (right-triangle, dark green), and 51.315 bohrs (circle, orange).</figcaption></figure>

lines with symbols: diamond (brown) for $L_z^{mat} = L_z^{atom}$, and circle, square, up-, and down-triangles (orange) for the other values. The spectra exhibit a completely different behavior: for $L_z^{atom}$ (Fig. 1, brown diamond), the peak of $\varepsilon_M^{LL}(q_z)$ appears at around 4 eV, while for all other values, it is located at 17 eV, with an amplitude scaled by a factor $1/L_z^{mat}$.

The peak at 4 eV for the in-plane component (Fig. 1, lines without symbols, blue) corresponds to the absorption resonance of the bulk silicon, and it is in agreement with the expected absorption spectrum as it results from the Lorentz model for thin slabs [18,22]. For the out-of-plane component (Fig. 1, lines with symbols, brown and orange), the result is more puzzling. It is not an effect of the effective-medium theory with vacuum, where the spectrum Im[$\varepsilon_M^{LL}(q_z; \omega)$] shifts regularly toward the plasmon peak when the size of the supercell increases [13]. Here, there is quite an abrupt jump from 4 to 17 eV. Actually, as can be seen in the inset of Fig. 1, the spectral weight is transferred from ~4 to ~17 eV on a lengthscale of around 10 bohrs. The peak at 4 eV for $L_z^{mat} = L_z^{atom}$ (Fig. 1, diamond, brown) means that for such a thickness, one performs an average on a bulklike unit cell. It also confirms that the Dyson equation given by Eq. (1) cuts the long-range component of the Coulomb potential regardless of the $q$ direction of the perturbation. For all the other values of $L_z^{mat}$, the peak is at 17 eV, corresponding to the plasmon frequency $\omega_p$ of the bulk silicon. This peak cannot be attributed to a surface plasmon, which is expected at $\omega_p/\sqrt{2}$ [25].

As could be expected from Fig. 1, there is a change in the nature of the spectra, which originates in the totally different shape of the ($q_\parallel \pm q_z$) spectra calculated to extract the $q_z$ component of $\varepsilon_M^{LL}(\omega)$ by linear combination [Eq. (A1)]. This puzzling effect is illustrated in Fig. 2 for the two values of $L_z^{mat} = L_z^{atom}$, corresponding to the situation in which the peak of the $q_z$ component appears at 4 eV (“zero vacuum”) and

<figure><img src="image_5.png" /><figcaption>FIG. 2. Im[$\varepsilon_M^{LL}(q_z; \omega)$] [Eq. (1)] for in-plane ($q_x$) [continuous (left) and dashed (right) blue curves], coupled ($q_x \pm q_z$) [up-triangle/continuous line (red) and down-triangle/dashed line (green), respectively], and out-of-plane excitations ($q_z$) [diamond/continuous line (brown) (left) and circle/dashed line (orange) (right)] extracted using Eq. (A1). Two thicknesses for the slab are considered: left: $L_z^{mat} = 41.052$ bohrs (zero vacuum) and right: $L_z^{mat} = 51.315$ bohrs (minimum vacuum).</figcaption></figure>
$L_z^{\text{mat}} = 51.315$ bohrs, which is the smallest value for which the peak for the $\mathbf{q}_z$ component appears at 17 eV, called for this reason “minimum vacuum.”

For $L_z^{\text{mat}} = L_z^{\text{atom}}$ (Fig. 2, left), the peak of the $(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ spectra (red up-triangle and green down-triangle) is slightly reduced in amplitude as compared to the $(\mathbf{q}_{\parallel})$ one, but it is still located at 4 eV and no structure at 17 eV is present: $\text{Im}[\varepsilon_M^{\text{LL}}(\mathbf{q}_z; \omega)]$ turns out to be ~4 eV (brown diamond). On the contrary, as soon as a slice of “vacuum” (~10 bohrs) is introduced in the slab thickness for the Selected-G procedure [13] (Fig. 2, right), the $(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ spectra show a peak at 4 eV reduced by a factor 2 as compared to the $(\mathbf{q}_{\parallel})$ one, and a peak at 17 eV arises. Consequently, $\text{Im}[\varepsilon_M^{\text{LL}}(\mathbf{q}_z; \omega)]$ exhibits a unique peak at 17 eV (orange circle). Actually, it can be seen that the induced density does not go to zero abruptly: this slice of ~10 bohrs corresponds to the two frontiers of the slab where the induced density, or the induced electric field (see Ref. [18], Fig. 2) or the response function (cf. Sec. III C), decrease to zero. This “minimum vacuum” case actually includes the interfaces.

All these spectra were calculated within the Selected-G formalism using the slab potential [13,16], which allows one to remove the artificial vacuum and calculate the response of an isolated object [16] of size defined according to $L_z^{\text{mat}}$. The remarkable difference of the spectra between the zero vacuum and minimum vacuum cases does not come from a spurious effect of vacuum, but from the fundamental change of electronic properties of the object we are describing. The zero vacuum case depicts a system that recovers the local electronic environment of the bulk, leading to $\bar{\chi}_{00}(\mathbf{q})$ having a resonance at the absorption frequency of the bulk material, whatever $\mathbf{q}$ is. On the contrary, the “minimum vacuum” case, due to the inclusion of the interfaces, depicts a quasi-2D object whose response function $\bar{\chi}_{00}(\mathbf{q})$ is highly anisotropic: it contains a peak at 4 eV, coming from the in-plane perturbation and corresponding to the absorption resonance, and a peak at 17 eV arising from the out-of-plane perturbation, and located at the plasmon resonance. One notes that, since this latter component can only be calculated from a linear combination of $(\mathbf{q}_{\parallel})$ and $(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ components, it turns out that the response function calculated for $(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ is a linear combination of what would be the absorption and loss spectrum of bulk silicon. This extremely puzzling spectrum directly results from our framework (Selected-G and Slab potential), which allows us to isolate the slab, remove the vacuum up to the bulk limit, and calculate the coupled $(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ components.

### III. A QUASI-2D OBJECT IN MIXED SPACE

#### A. In-plane and out-of-plane components of the absorption spectrum

To have a better understanding of the properties of a quasi-2D object, and in particular the role of the extra slice beyond $L_z^{\text{atom}}$, we solved the Dyson equation in a mixed-space approach [13,18]. In this framework, the in-plane components are described in reciprocal space $(q_x, q_y)$, since the system is still infinite and periodic, and the direction perpendicular to the slab, where the matter has been cut, is treated in real space $(z)$. The Dyson equation reads (we omit the frequency for the

sake of clarity)

$$ \chi_{G_{\parallel}G'_{\parallel}}(\mathbf{q}_{\parallel}, z, z') = \chi_{G_{\parallel}G'_{\parallel}}^{0}(\mathbf{q}_{\parallel}, z, z') + \sum_{G_{1\parallel}} \int dz_1 \int dz_2 \chi_{G_{\parallel}G_{1\parallel}}^{0}(\mathbf{q}_{\parallel}, z, z_1) \times v_{G_{\parallel}}(\mathbf{q}_{\parallel}, z_1, z_2) \chi_{G_{\parallel}G'_{\parallel}}(\mathbf{q}_{\parallel}, z_2, z'), $$

where $v_{G_{\parallel}}(\mathbf{q}_{\parallel}, z_1, z_2)$ is the two-dimensional Fourier transform of the Coulomb potential:

$$ v_{G_{\parallel}}(\mathbf{q}_{\parallel}, z, z') = \frac{2\pi}{|G_{\parallel} + \mathbf{q}_{\parallel}|} e^{-|G_{\parallel} + \mathbf{q}_{\parallel}|z - z'}, \quad (3) $$

and $\chi^0$ is the Kohn-Sham response function. Moreover, since we expect the in-plane local field effects to be small [26], we hypothesize further that only the $G_{\parallel} = G'_{\parallel} = 0$ case is sufficient to give reliable results (from now on, we omit the in-plane $G_{\parallel} = 0_{\parallel}$), leading to

$$ \chi(\mathbf{q}_{\parallel}, z, z') = \chi^0(\mathbf{q}_{\parallel}, z, z') + \int dz_1 \int dz_2 \chi^0(\mathbf{q}_{\parallel}, z, z_1) \times v(\mathbf{q}_{\parallel}, z_1, z_2) \chi(\mathbf{q}_{\parallel}, z_2, z'). \quad (4) $$

In Ref. [13], since our goal was to study the optical properties of a surface, depicted as a thick slab, the wavelength of the perturbation was rescaled to ~1 in order to fulfill the condition $|\mathbf{q}_{\parallel}|L/2 \gg 1$ [27]. In Ref. [18], the purpose was to study the dielectric properties of an ultrathin slab, so we kept a vanishing value for $|\mathbf{q}|$. Moreover, to obtain independently the dielectric function for the component perpendicular to the slab, we did the further approximation that the slab was homogeneous in-plane, and we replaced Eq. (3) with $v(\mathbf{q}_{\parallel} \equiv 0, z, z') = -2\pi|z - z'|$. In this paper, we pursue the study of a quasi-2D object ($|\mathbf{q}_{\parallel}|L/2 \ll 1$), but using the usual 2D Coulomb potential given by Eq. (3) to compare with the reciprocal space with the slab potential framework. The Kohn-Sham response function $\chi^0$ is first evaluated in reciprocal space with the DP code [28] for a set of reciprocal-lattice vectors $(G = (0_{\parallel}, G_z))$, from the ground-state energies and eigenfunctions calculated with ABINIT [29,30] for the supercell. The direction of the perturbation is set through a vanishing $\mathbf{q}$ vector, used within $\mathbf{k} \cdot \mathbf{p}$ theory for the head and the wings (first line and column of the matrix). We then perform the inverse Fourier transform (IFT). Technical details and tests of consistency are given in Appendix B. The IFT of $\chi_{G_z, G'_z}^0(\mathbf{q}_{\parallel})$ reads

$$ \chi^0(\mathbf{q}_{\parallel}, z_n, z_{n'}) = \frac{1}{L_z} \sum_{m,m'} e^{iG_m z_n} \chi_{G_m G_{m'}}^0(\mathbf{q}_{\parallel}) e^{-iG_{m'} z_{n'}}. \quad (5) $$

As was the case for the slab potential in reciprocal space, the expression for the Coulomb potential [Eq. (3)] for $G_{\parallel} = 0$ diverges if $\mathbf{q}_{\parallel} \equiv 0$. Thus the $\mathbf{q}_z$ component must be extracted from a linear combination of different components (Appendix A). The macroscopic average of the response function in reciprocal space for the coupled $(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ perturbation is recovered by multiplying $\chi^0(\mathbf{q}_{\parallel}, z_n, z_{n'})$ by $e^{-iq_z(z_n - z_{n'})}$ leading to

$$ \chi_M^0(\mathbf{q}_{\parallel}, q_z) = \frac{\Delta z_z^2}{L_z} \sum_{n,n'} e^{-iq_z z_n} \chi^0(\mathbf{q}_{\parallel}, z_n, z_{n'}) e^{iq_z z_{n'}}. \quad (6) $$
<figure><img src="image_6.png" /><figcaption>FIG. 3. Imaginary part of the macroscopic average [Eq. (6)] of $\chi^0(\mathbf{q}_x, z, z')$ multiplied by $-4\pi/|\mathbf{q}|^2$, calculated for $\mathbf{q}_x$ ($\mathbf{q}_z \equiv 0$) (violet square) or $\mathbf{q}_x \pm \mathbf{q}_z$ [down-triangle, light green (+) and right-triangle, dark green (-)] ($q_z \neq 0$). The $\mathbf{q}_z$ contribution (magenta diamond) is obtained from linear combination [Eq. (A1)]. Calculations for $||\mathbf{q}_x|| = ||\mathbf{q}_z|| = 1 \times 10^{-4}$ a.u., but since, in the optical limit, $\chi_{00}^0(\mathbf{q}) \equiv \chi_M^0(\mathbf{q}_\|, q_z)$ is proportional to $|\mathbf{q}|^2$, the results are independent of the value of $|\mathbf{q}|$.</figcaption></figure>

In Eq. (6), the sum over $z_n$ is done over the full size of the simulation box defined by $L_z$. $\chi_M^0(\mathbf{q}_\|, q_z)$ can also be obtained as the FT of $\chi^0(\mathbf{q}_\|, z_n, z_{n'})$ [Eq. (B1) for $G_m = G_{m'} = 0$].

In the optical limit ($\mathbf{q} \to 0$), $\chi_{00}^0(\mathbf{q}; \omega)$ is related to the absorption spectrum when the local field effects are negligible, according to $\epsilon_M^{\text{NLF}}(\mathbf{q}; \omega) = 1 - 4\pi/|\mathbf{q}|^2\chi_{00}^0(\mathbf{q}; \omega)$, where NLF stands for nonlocal fields. Then, one can extract $\epsilon_M^{\text{NLF}}(\mathbf{q}_z; \omega)$ using the properties of the dielectric tensor $\overleftrightarrow{\epsilon}_M(\omega)$ (see Appendix A). The result is presented in Fig. 3.

Applying the linear combination [Eq. (A1)] to spectra calculated for ($\mathbf{q}_x$) (square, violet), ($\mathbf{q}_x + \mathbf{q}_z$) (down-triangle, light green), and ($\mathbf{q}_x - \mathbf{q}_z$) (right-triangle, dark green), we obtain the spectrum $\text{Im}[\epsilon_M^{\text{NLF}}(\mathbf{q}_z; \omega)]$ (diamond, magenta). We recover the spectrum calculated in reciprocal space directly for the $\mathbf{q}_z$ component (Fig. 12, right).

This result validates the procedure of calculating the out-of-plane component of the macroscopic dielectric function within this mixed-space approach from alternatively using the macroscopic average [Eq. (6)] or the Fourier transform [Eq. (B1) for $G_m = G_{m'} = 0$], followed by the linear combination [Eq. (A1)].

One notes that the absorption spectra calculated neglecting the local field effects contain one peak at 4 eV regardless of the direction of the perturbation $\mathbf{q}$ (Fig. 3).

## B. Susceptibility without the long-range component of the Coulomb potential in mixed-space

For three-dimensional crystals, the absorption spectrum can be calculated in reciprocal space as the optical limit ($\mathbf{q} \to 0$) of the imaginary part of

$$ \epsilon_M(\mathbf{q}; \omega) = 1 - \frac{4\pi}{|\mathbf{q}|^2} \bar{\chi}_{00}(\mathbf{q}; \omega), \quad (7) $$

where $\bar{\chi}$ is the response function to the macroscopic part of the total potential. It is calculated by solving a modified Dyson equation (we omit the frequency $\omega$ for clarity):

$$ \bar{\chi}_{GG'}(\mathbf{q}) = \chi_{GG'}^0(\mathbf{q}) + \sum_{\mathbf{G}_1} \chi_{GG_1}^0(\mathbf{q}) \bar{v}_{G_1}(\mathbf{q}) \bar{\chi}_{G_1 G'}(\mathbf{q}), \quad (8) $$

$$ \begin{gathered} \text{where } \bar{v}_{G}(\mathbf{q}) = 0 \text{ for } \mathbf{G} = \mathbf{0} \text{ and } \\ \bar{v}_{G}(\mathbf{q}) \equiv v_{G}(\mathbf{q}) = \frac{4\pi}{||\mathbf{G} + \mathbf{q}||^2} \text{ for } \mathbf{G} \neq \mathbf{0}. \end{gathered} \quad (9) $$

In Eq. (8), contrary to the regular Dyson equation [with full $v_G$, Eq. (C2)], we suppress the long-range component of the Coulomb potential $v_0$, which is at the origin of the difference between the plasmon and the absorption resonances in infinite systems. For OD systems, this difference naturally vanished.

In this paper, since we deal with quasi-2D systems which are still infinite in the $(x, y)$ directions, the question of the influence of a long-range component occurs. The approach resulting from Eqs. (8) and (9) is well suited for reciprocal space, but it is less obvious in real space. To achieve this goal, we developed a procedure (see Appendix C) that consists in (i) the definition, in reciprocal space, of a new $\bar{\chi}^0$ identical to $\chi^0$, except for the first column which is set to zero: $\bar{\chi}_{G0}^0 \equiv 0, \forall \mathbf{G}$; (ii) the calculation of the inverse Fourier transform of both $\chi^0$ and $\bar{\chi}^0$ to get $\chi^0(\mathbf{q}_\|, z, z')$ and $\bar{\chi}^0(\mathbf{q}_\|, z, z')$; and (iii) the resolution of

$$ \begin{aligned} \bar{\chi}(\mathbf{q}_\|, z, z') = \chi^0(\mathbf{q}_\|, z, z') &+ \int dz_1 \int dz_2 \bar{\chi}^0(\mathbf{q}_\|, z, z_1) \\ & \times v(\mathbf{q}_\|, z_1, z_2) \bar{\chi}(\mathbf{q}_\|, z_2, z') \end{aligned} \quad (10) $$

with $v(\mathbf{q}_\|, z, z') = \frac{2\pi}{|\mathbf{q}_\|} e^{-|\mathbf{q}_\||z-z'|}$. The difference between $\bar{\chi}_{00}(\mathbf{q})$ and $\chi_{00}(\mathbf{q})$ resulting from the solution of the Dyson equation in real space (Fig. 13) is detailed in Appendix C for the structure at 4 eV.

While $\chi_{00}^0(\mathbf{q}_\|)$ and $\chi_{00}^0(\mathbf{q}_\| \pm \mathbf{q}_z)$ present only one peak at 4 eV (Fig. 3), the introduction of the local field effects through the solution of the Dyson equation for $\bar{\chi}$ leads to the appearance of a peak at 17 eV (Fig. 4, top) for ($\mathbf{q}_x \pm \mathbf{q}_z$) (magenta left-triangle, orange right-triangle).

Moreover, one can see that the peak at 17 eV in $\bar{\chi}_{00}(\mathbf{q}_x \pm \mathbf{q}_z)$ (magenta left-triangle, orange right-triangle) is strictly the same as the one in $\chi_{00}(\mathbf{q}_x \pm \mathbf{q}_z)$ (red square, green diamond). The small shift observed for the peak at 4 eV (Fig. 13) does not exist. Applying the linear combination to $\bar{\chi}_{00}$ [Eq. (A3)] as well as to $\chi_{00}$, we can extract exactly the same peak for the $\mathbf{q}_z$ component: $\bar{\chi}_{00}(\mathbf{q}_z) \equiv \chi_{00}(\mathbf{q}_z)$. It is not obvious that one can apply the linear combination [Eq. (A3)] to the $\chi_{00}$ response function. In three-dimensional materials, this quantity is related to the loss function $\epsilon_{00}^{-1}(\mathbf{q}) = 1 + \frac{4\pi}{|\mathbf{q}|^2} \chi_{00}(\mathbf{q})$, which is not expected to behave as linear combinations of spectra for the different components. It is $\epsilon_M = \lim_{q \to 0} [1/\epsilon_{00}^{-1}(\mathbf{q})]$, which verifies this property. For the case of the thin slab, the similarity between $\bar{\chi}$ and $\chi$ seems to indicate that the linear combinations of spectra [Eq. (A3)] can be applied to both response functions. The small dependence in $|\mathbf{q}|$ of the resonance around $\omega_0$ does not affect the extracted $\mathbf{q}_z$ component. Moreover, it can be verified that the peak at 17 eV is exactly the same as the one obtained in the direct calculation of
<figure><img src="image_7.png" /><figcaption>FIG. 4. Imaginary part of FT [Eq. (B1) for $G_m = G_{m'} = 0$] of $\bar{\chi}(\mathbf{q}_x, z, z')$ and $\chi(\mathbf{q}_x, z, z')$, respectively, after resolution of the Dyson equations [Eqs. (10) and (4)]. Spectra for $\mathbf{q}_x$ (blue up-triangle and black circle, respectively) and $(\mathbf{q}_z \pm \mathbf{q}_x)$ (magenta left-triangle/orange right-triangle and red square/green diamond, respectively). Bottom: Linear combination [Eq. (A3)] for $\bar{\chi}_{00}(\mathbf{q})$ and $\chi_{00}(\mathbf{q})$ (calculated in mixed-space), to extract the $\mathbf{q}_z$ component (violet down-triangle and turquoise cross, respectively) compared to the $\mathbf{q}_z$ component coming from the direct calculation of Ref. [18] (brown diamond). Calculation for the slab where the atomic positions extend on 41.052 bohrs, introduced in a supercell of height 82.104 bohrs.</figcaption></figure>

$\chi_{00}(\mathbf{q}_z)$ in the case of a planar homogeneous slab (diamond brown) (see Ref. [18]). These results show that for an isolated 2D object, $\bar{\chi}_{00}(\mathbf{q}_{\perp}; \omega) \equiv \chi_{00}(\mathbf{q}_{\perp}; \omega)$: there is no effect of the long-range component of the Coulomb potential for the out-of-plane perturbation in the optical limit. The same result could have been obtained in reciprocal space with the slab potential. Such a finding is finally less surprising than for the case of the in-plane component, which probes the infinite directions, since for $\mathbf{q}_{\perp}$, the perturbation is in the direction where the matter was cut. The surprising point is that the quantity $\bar{\chi}_{00}(\mathbf{q}_{\perp}; \omega)$, which should be related to the absorption spectrum, exhibits a resonance at the plasma frequency. In this calculation, done in the simulation box defined by $L_z$ much

larger than $L_z^{\text{mat}}$, where the interfaces are naturally present, we confirm the puzzling result obtained in reciprocal space for the minimum vacuum case (see Sec. II).

### C. How do we define the thickness of the slab?

Up to now, the spectra (Fig. 3) have been obtained performing a macroscopic average [Eq. (6)] of quantities calculated in mixed-space in the full simulation box, resulting from the inverse Fourier transform (IFT) of the $\chi^0$ calculated in reciprocal space in a supercell. From the IFT process, the size of the box is the height of the supercell $L_z$ (Fig. 5, top, left). As a consequence, the amplitude of spectra $\epsilon_M^{\text{NLF}}$ (Fig. 5, bottom, left) suffers from the same problem as in reciprocal space, i.e., it is scaled by the factor $1/L_z$, which appears in Eq. (6). Indeed, $\chi^0(\mathbf{q}_{\parallel}, z_n, z_{n'})$ [Eq. (5)] is independent of the height of the supercell used in reciprocal space, and $\bar{\chi}(\mathbf{q}_{\parallel}; z_n, z_{n'})$ and $\chi(\mathbf{q}_{\parallel}; z_n, z_{n'})$ as well. Thus the absorption spectrum is smaller by a factor $L_z^{\text{mat}}/L_z$ as compared to the one for bulk silicon, where $L_z^{\text{mat}}$ is the thickness of the slab. Beyond the technical point, we know that it is a manifestation of the effective-medium theory with vacuum [13]. Moreover, the question of the normalization of the Kohn-Sham response function, when it is introduced in the Dyson equation, is also an issue.

In Ref. [13], to depict a surface, $L_z^{\text{mat}}$ was defined according to the atomic positions, and $\chi^0$ in reciprocal space was scaled with the ratio $L_z/L_z^{\text{atom}}$. Nevertheless, it has been shown in Refs. [18,22] that for a thin slab, this quantity cannot be defined so easily due to the fact that the induced density and electric field extend in $(z, z')$ much larger than the extremes of the atomic positions, or in other words, that the thickness and the role of the interfaces are not negligible. To define $L_z^{\text{mat}}$, we propose to adjust the size of the simulation box at the limits where the shape of the spectra is recovered. This value must result directly from the inverse Fourier transform procedure. Since the size of the box is defined by the length of the smallest $G_z$ vector involved in the IFT, we propose to apply the Selected-G procedure [13] to the IFT, and to replace Eq. (5) with

$$ \chi^0(\mathbf{q}_{\parallel}, z_n, z_{n'}) = \frac{1}{L_z^{\text{mat}}} \sum_{m,m'} e^{i\tilde{G}_m z_n} \chi_{\tilde{G}_m \tilde{G}_{m'}}^0(\mathbf{q}_{\parallel}) e^{-i\tilde{G}_{m'} z_{n'}}, \quad (11) $$

where $\tilde{G}_m = \frac{2\pi}{L_z^{\text{mat}}}$, and no longer $\frac{2\pi}{L_z}$.

The result is shown in Fig. 5 for two values of $L_z^{\text{mat}}$: $L_z^{\text{mat}} = L_z^{\text{atom}} = 41.052$ bohrs, corresponding to the extension of the atomic positions (center), and $L_z^{\text{mat}} = 51.315$ bohrs, corresponding to the extension of the response functions (right). In both cases, the ratio between $L_z$ and $L_z^{\text{mat}}$ was chosen to be 2, which means that the original $\chi_{\tilde{G}_x \tilde{G}_x}^0(\mathbf{q}_{\parallel})$ were evaluated in a supercell of height $L_z = 84.104$ and 102.63 bohrs, respectively. In Fig. 5 (top, center, and right), the response functions have not been cut by the plot, but they are defined as they are plotted: the length of the simulation box has been resized by the Selected-G procedure to the value that we will identify as $L_z^{\text{mat}}$.

The value of $L_z$ for the initial supercells was chosen in order to contain the slab of silicon plus the vacuum needed to allow the apparition of the surface states, and the isolation of the replicas at the level of the DFT calculations. This value is
<figure><img src="image_8.png" /><figcaption>FIG. 5. Top: Im[$\chi^0(\mathbf{q}_\|, z, z')$] resulting from the IFT of $\chi^0_{G_z; G'_z}(\mathbf{q}_x)$ for a box of size $L_z = 82.104$ bohrs, coming from the supercell height (left), a box of size corresponding to $L_z^{\text{mat}} = L_z^{\text{atom}} = 41.052$ bohrs (center), and to $L_z^{\text{mat}} = 51.315$ bohrs (right). The last two cases result from the Selected-G procedure (ratio 2) applied when performing the inverse Fourier transform of $\chi^0_{G_z; G'_z}(\mathbf{q}_x)$ [Eq. (11)] calculated with the original supercells of size $L_z = 82.104$ and 102.63 bohrs, respectively. The different symbols correspond to different values of $z'$. Bottom: Imaginary part of $\epsilon_M^{\text{NLF}}(\mathbf{q}; \omega) = 1 - 4\pi/|\mathbf{q}|^2\chi^0_{00}(\mathbf{q}; \omega)$, where $\chi^0_{00}(\mathbf{q}; \omega)$ is the macroscopic average of $\chi^0(\mathbf{q}_x, z, z'; \omega)$ [Eq. (6)]. For the three cases, the summation over $z_n$ extends up to the full simulation box, which has been resized to the dedicated value of $L_z^{\text{mat}}$.</figcaption></figure>

overestimated, but it must be at least twice the thickness of the silicon slab to apply the Selected-G procedure, which requires an integer ratio between $L_z$ and $L_z^{\text{mat}}$. A choice of $L_z = 51.315$ bohrs, which would have avoided the Selected-G procedure, is, on the contrary, not large enough to provide noninteracting replicas in the calculation of the Kohn-Sham structure. The left panels correspond to the IFT and FT without Selected-G, plotted for reference.

We can immediately see that for $L_z^{\text{mat}} = 41.052$ bohrs (Fig. 5, top, center), $\chi^0(\mathbf{q}_\|, z, z')$ is folded in the borders of the box. This is particularly visible for $z' = \pm 20.50$ bohrs (open black circle and plain cyan circle, respectively) where the curves seem to add to each other, instead of slowly going back to zero (Fig. 5, top, left). As a consequence, the calculation of the macroscopic average [Eq. (6), for which the summation over $z_n$ extends up to the limit of the response function, which corresponds to the new size of the simulation box] does not allow one to recover the original spectrum $\epsilon_M^{\text{NLF}}(\mathbf{q}_z)$ (dashed blue line in Fig. 5, bottom, center), as can be seen from a comparison with the reference spectrum (Fig. 5, bottom, left). Indeed, the $\mathbf{q}_x \pm \mathbf{q}_z$ spectra (red up-triangle and green down-triangle in Fig. 5, bottom, center) are wrong, since, due to the folding, the $z$ and $z'$ values appearing in the phase factor $e^{-iq_z(z-z')}$ are not associated with the correct corresponding values of $\chi^0(\mathbf{q}_\|, z, z')$. On the contrary, for $L_z^{\text{mat}} = 51.315$ bohrs (Fig. 5, top, right), $\chi^0(\mathbf{q}_\|, z, z')$ exhibits the correct development, no folding occurs, all the $z'$ cuts reach zero in the border of the simulation box, and the $\mathbf{q}_z$ spectrum is

recovered (dashed blue line in Fig. 5, bottom, right). Moreover, the amplitude of the spectra $\epsilon_M^{\text{NLF}}(\mathbf{q})$ (Fig. 5, bottom, right) is now of the order of magnitude of the one of the dielectric function of silicon. The exact amplitude of the dielectric function of silicon bulk would have been obtained by performing the summation over any box large enough to allow the full development of $\chi^0(\mathbf{q}_\|; z, z')$ [since the numerator of Eq. (6) is independent of the size of the box], and normalizing this quantity by $L_z^{\text{atom}} = 41.052$ bohrs, corresponding to the “zero vacuum” case. Nevertheless, if one defines the thickness from the simulation box with the minimum size, which allows one to recover the original $\chi^0_{00}(\mathbf{q})$ (see Fig. 5), it cannot be chosen smaller than the extension where the Kohn-Sham response function goes to zero. In our calculation, it is achieved for $L_z^{\text{mat}} = 51.315$ bohrs, where $\chi^0_{00}(\mathbf{q}_\|)$ is of the order of $5 \times 10^{-5}$ (Fig. 5, right panels).

This justifies the term “minimum vacuum” given to this case. This value actually accounts for the presence of the interfaces, and we consider that it is the correct way to define the thickness of the quasi-2D object. We emphasize that it is different from the usual practice [32,33], which indeed corresponds to taking for the thickness the value $L_z^{\text{atom}}$, in the spirit of a van der Waals atomic radius, which is also quite close to the value given by the ground-state density [16], [34]. From a practical point of view, it should correspond to the value where the $\chi^0(\mathbf{q}_\|; z, z')$ reaches zero. We evaluate this extra slice beyond the extreme atomic positions to be of the order of magnitude of 5–6 bohrs per interface. It has been
established in the case of a silicon slab, but the scheme also affects traditional 2D materials. Indeed, Fig. 1 of Ref. [31] shows that the induced density in one layer of MoS₂ extends beyond the atoms, and has the same features as Fig. 2 of Ref. [18]. This will be evidenced in Sec. V on a stacking of four graphene layers.

### D. Macroscopic dielectric function versus longitudinal-longitudinal contraction of $\vec{\epsilon}_M(\omega)$ for a quasi-2D object

Now that we have defined a thickness that accounts correctly for the electronic properties of the slab, we solve the Dyson equation [Eq. (10)] with the corresponding $\chi^0(\mathbf{q}_{\|}; z, z')$ (Fig. 5, top, right) [and $\bar{\chi}^0(\mathbf{q}_{\|}; z, z')$], obtained by applying the Selected-G procedure when performing the inverse Fourier transform of $\chi^0_{G_z G'_z}(\mathbf{q}_x)$ [and $\bar{\chi}^0_{G_z G'_z}(\mathbf{q}_x)$] [Eq. (11)], in order to adjust the size of the simulation box at the “minimum vacuum” case ($L_z^{\text{mat}} = 51.315$ bohrs). The result is presented in Fig. 6, which depicts the imaginary part of $[1 - \frac{4\pi}{|\mathbf{q}|^2} \bar{\chi}_{00}(\mathbf{q}; \omega)]$ labeled $\epsilon_M^{\text{LL}}(\mathbf{q}; \omega)$, where $\bar{\chi}_{00}(\mathbf{q}; \omega)$ is the macroscopic average [Eq. (6)] of $\bar{\chi}(\mathbf{q}_{\|}; z, z')$. The summation is done over the simulation box resized to the value $L_z^{\text{mat}} = 51.315$ bohrs.

The first observation is that the response function $\bar{\chi}_{00}(\mathbf{q}; \omega)$ has exactly the same spectral shape as if they were calculated on a box resulting from the original supercell height (Fig. 4, top), showing that the Selected-G process at the level of IFT [Eq. (11)] does not introduce any spurious effect. The advantage is that, with this process, the amplitude of the spectra is associated with the thickness of the slab, and no longer to the size of the supercell. Contrary to $\chi_{00}^0(\mathbf{q})$ (Fig. 3), $\bar{\chi}_{00}(\mathbf{q})$ exhibits a peak at 4 eV for the in-plane excitation (Fig. 6, continuous black line), and two peaks at 4 and 17 eV for the coupled $(\mathbf{q}_x \pm \mathbf{q}_z)$ excitation (Fig. 6, red up-triangle and green down-triangle), leading to a contribution at 17 eV for

<figure><img src="image_9.png" /><figcaption>FIG. 6. Im[$\epsilon_M^{\text{LL}}(\mathbf{q}; \omega)$] [Eq. (7)] with $\bar{\chi}_{00}(\mathbf{q})$ calculated according to Eqs. (10) and (6), for $\mathbf{q}_x$ (continuous black line), $(\mathbf{q}_x + \mathbf{q}_z)$ (red up-triangle), $(\mathbf{q}_x - \mathbf{q}_z)$ (green down-triangle), and $\mathbf{q}_z$ (blue dashed line) [Eq. (A1)]. Calculation for $L_z^{\text{mat}} = 51.315$ bohrs with Selected-G procedure applied when performing the inverse Fourier transform of $\chi^0_{G_z G'_z}(\mathbf{q}_x)$ [Eq. (11)].</figcaption></figure>

the extracted $\mathbf{q}_z$ spectrum [Eq. (A1)] (Fig. 6, dashed blue line). This spectrum results from the addition of the very thin slice of vacuum to $L_z^{\text{atom}}$, allowing us to rule out that it is due to a mechanism of the effective-medium theory with vacuum [11,13,18].

We recover the results obtained within the reciprocal space calculations (Sec. II, Fig. 2, right). The presence of the peak at 17 eV, in the absence of the long-range component of the Coulomb potential, when the perturbation is perpendicular to the slab, is actually a consequence of the inclusion of the interfaces, and is a characteristic feature of longitudinal-longitudinal contraction of the dielectric tensor of the quasi-2D object. We conclude from these results that Ehrenreich’s demonstration [21], established for 3D (cubic) systems, which shows that in the optical limit the longitudinal-longitudinal contraction of the dielectric tensor is equal to the transverse-transverse contraction (allowing the calculation of the absorption spectrum from the longitudinal calculation), is no longer valid for 2D objects. Indeed, the resonance of Im[$\epsilon_M^{\text{LL}}(\mathbf{q}_{\|}; \omega)$] occurs at the plasma frequency, while the resonance of the dielectric function is still expected at the optical absorption frequency (4 eV for silicon). This failure comes from the finiteness of the object, and it should also occur for 1D and 0D objects.

### E. Angular dependence of the longitudinal macroscopic dielectric function for a quasi-2D object

Up to now, the presented calculations were done for $|\mathbf{q}_x| = |\mathbf{q}_z|$, corresponding to an angle between the direction of the perturbation and the interfaces of the slab of $\theta = 45^\circ$. Based on the direct calculation of $\epsilon_M^{\text{LL}}(\mathbf{q}_{\|})$ (Fig. 6, continuous black

<figure><img src="image_10.png" /><figcaption>FIG. 7. Imaginary part of [$\epsilon_M^{\text{LL}}(\mathbf{q}_x + \mathbf{q}_z) + \epsilon_M^{\text{LL}}(\mathbf{q}_x - \mathbf{q}_z)]/2$ for different angles $\theta$ of $\mathbf{q}$ with the surface plane direction (x): light green circle ($\theta = 30^\circ$), blue square ($\theta = 45^\circ$), and magenta diamond ($\theta = 60^\circ$) compared with the linear combination Im[$\epsilon_M^{\text{LL}}(\mathbf{q}_x)]\cos^2(\theta) + \text{Im}[\epsilon_M^{\text{LL}}(\mathbf{q}_z)]\sin^2(\theta)$ with Im[$\epsilon_M^{\text{LL}}(\mathbf{q}_x)$] (Fig. 6, continuous black) and Im[$\epsilon_M^{\text{LL}}(\mathbf{q}_z)$] (Fig. 6, dashed blue) for $\theta = 30^\circ$ (dark green cross), $\theta = 45^\circ$ (cyan up-triangle), and $\theta = 60^\circ$ (orange down-triangle). Calculations are done in real space for $L_z^{\text{mat}} = 51.315$ bohrs (“minimum vacuum” case) with Selected-G procedure applied when performing the inverse Fourier transform of $\chi^0_{G_z G'_z}(\mathbf{q}_x)$ [Eq. (11)].</figcaption></figure>
line), and using the linear combination [Eq. (A1)], we extracted what appears to be $\epsilon_M^{LL}(\mathbf{q}_z)$ (Fig. 6, dashed blue line). It can be checked that the spectrum $\epsilon_M^{LL}(\mathbf{q}_z)$ is independent of the angle chosen. Indeed, Eq. (A1) can be rewritten as

$$ \frac{\epsilon_M^{LL}(\hat{\mathbf{q}}_x + \hat{\mathbf{q}}_z) + \epsilon_M^{LL}(\hat{\mathbf{q}}_x - \hat{\mathbf{q}}_z)}{2} = \epsilon_M^{LL}(\hat{\mathbf{q}}_x) \cos^2 \theta + \epsilon_M^{LL}(\hat{\mathbf{q}}_z) \sin^2 \theta. $$

In Fig. 7, we compare the sum of the spectra calculated for $(\mathbf{q}_x \pm \mathbf{q}_z)$ for different values of $|\mathbf{q}_x|$ and $|\mathbf{q}_z|$ corresponding to $\theta = 30^{\circ}$ (light green circle), $\theta = 45^{\circ}$ (blue square), and $\theta = 60^{\circ}$ (magenta diamond) with $\text{Im}[\epsilon_M^{LL}(\mathbf{q}_x)] \cos^2(\theta) + \text{Im}[\epsilon_M^{LL}(\mathbf{q}_z)] \sin^2(\theta)$ (dark green cross for $\theta = 30^{\circ}$, cyan up-triangle for $\theta = 45^{\circ}$, and orange down-triangle for $\theta = 60^{\circ}$), where $\text{Im}[\epsilon_M^{LL}(\mathbf{q}_x)]$ and $\text{Im}[\epsilon_M^{LL}(\mathbf{q}_z)]$ correspond to the spectra plotted as a continuous black line and a dashed blue line in Fig. 6. Such a result would have been considered trivial for a bulk material, but in the case of the slab, the quantity that we label $\epsilon_M^{LL}(\mathbf{q}_z)$ corresponds to the plasmon of the bulk counterpart, and it cannot describe the absorption of an optical photon in silicon.

## IV. OPTICAL PROPERTIES OF A QUASI-2D OBJECT: AIRY'S FORMULAS

To understand this confusing result, namely the linear combination of the absorption peak and the plasmon peak in a quantity that we expected to represent the optical response of our quasi-2D object, we investigated the problem by means of classical electromagnetism. Indeed, for lack of experimental results, at least to our knowledge, we considered that Airy's formula should be representative of the measured quantity.

Let us consider first a surface separating the vacuum and a semi-infinite medium having a dielectric function $\epsilon(\omega)$. The reflection coefficient of an electromagnetic field of wave vector $\mathbf{k}_i$ impinging the surface with an angle $\theta$ (between $\mathbf{k}_i$ and the normal to the surface) is the so-called Fresnel coefficient [35–37], which for the p-polarization of light reads

$$ r_p \equiv \frac{\epsilon \cos \theta - \sqrt{\epsilon - \sin^2 \theta}}{\epsilon \cos \theta + \sqrt{\epsilon - \sin^2 \theta}}. \qquad (12) $$

We omit the frequency dependence for clarity. We focus on the p-polarization since it depicts the case in which the electric field $\mathbf{E}$ lies in the incident plane, leading to a component perpendicular to the surface ($\mathbf{E}_z$), in addition to the parallel one ($\mathbf{E}_\parallel$) for $\theta$ different from zero.

If the medium now has a finite thickness $d$, the electromagnetic field will be subjected to multiple reflections and transmissions on the two interfaces, and the coefficients for the reflection and transmission are given by Airy's formulas [36,38–40]:

$$ r^{\text{slab}} = r \frac{1 - e^{2i\beta}}{1 - r^2 e^{2i\beta}}, \qquad (13) $$

$$ t^{\text{slab}} = \frac{(1 - r^2) e^{i\beta}}{1 - r^2 e^{2i\beta}}, \qquad (14) $$

where $\beta = k_z d$, with $k_z$ the component perpendicular to the interface of the wave vector transmitted in the medium:

$k_z = k_0 \sqrt{\epsilon - \sin^2 \theta}$, with $k_0 = \omega/c$ the wave vector in vacuum [36].

To evaluate the reflectance $R = |r^{\text{slab}}|^2$, the transmittance $T = |t^{\text{slab}}|^2$, and the absorbance $A = 1 - R - T$, we used for $\epsilon(\omega)$ a model Lorentz function:

$$ \epsilon(\omega) = 1 - \frac{\omega_p^2}{\omega^2 - \omega_0^2 + i\omega\gamma} $$

with $\omega_0 = 4$ eV, $\omega_p = 17$ eV, and $\gamma = 1.9$ eV, in order to reproduce the dielectric function of bulk silicon.

<figure><img src="image_11.png" /><figcaption>FIG. 8. Spectra of reflectance $R = |r_p^{\text{slab}}|^2$ with $r_p^{\text{slab}}$ given by Eq. (13) (continuous black line), transmittance $T = |t_p^{\text{slab}}|^2$ with $t_p^{\text{slab}}$ given by Eq. (14) (red up-triangle), and absorbance $A = 1 - R - T$ (green down-triangle) calculated with Airy's formula, compared with the reflection on the surface $|r_p|^2$ (Fresnel, blue dashed line) for $\theta = 30^{\circ}$ (top), $\theta = 45^{\circ}$ (center), and $\theta = 60^{\circ}$ (bottom). Thickness: $d = 20 \AA$.</figcaption></figure>
The results are presented in Fig. 8. While the reflection on a surface [Eq. (12)] is a plateau between $\omega_0$ and $\omega_p$ (dashed blue lines), the reflectance and transmittance calculated with Airy's formula are completely different. The reflectance (black continuous line), transmittance (red up-triangle), and absorbance (green down-triangle) are composed of two peaks, one located at $\omega_0 = 4$ eV, and one located at $\omega_p = 17$ eV. The respective spectral weight of these two features evolves with the incident angle $\theta$.

It can be shown that, for a thin slab, the reflection and transmission coefficients can be expended at first order in $\beta$, which gives (see Appendix D)

$$r_p^{\text{slab}} \approx \frac{-2ik_0 d [(\varepsilon - 1) \cos^2 \theta + (1/\varepsilon - 1) \sin^2 \theta]}{4 \cos \theta - 2ik_0 d [(\varepsilon - 1) \cos^2 \theta - (1/\varepsilon - 1) \sin^2 \theta]}, \quad (15)$$

$$i_p^{\text{slab}} \approx \frac{4 \cos \theta}{4 \cos \theta - 2ik_0 d [(\varepsilon - 1) \cos^2 \theta - (1/\varepsilon - 1) \sin^2 \theta]} \quad (16)$$

These expressions are the key to understanding the structure of the reflectance, transmittance, and absorbance spectra (Fig. 8). The peak at $\omega_0 = 4$ eV arises from the $(\varepsilon - 1)$ contribution, and the peak at $\omega_p = 17$ eV from the $(1/\varepsilon - 1)$ term. The angular dependence is also perfectly explained. The comparison of $R = |r_p^{\text{slab}}|^2$ with $r_p^{\text{slab}}$ given by Eq. (13) (Fig. 8, continuous black line) with the reflectance on the semi-infinite media [Eq. (12)] (Fig. 8, dashed blue line) shows that, for the thin slab, the part of the intensity of the electric field that is reflected is strongly reduced. The introduction of the second interface, allowing reflections and transmissions, leads to some destructive interferences. Moreover, the sum of the reflectance and the transmittance at 17 eV is not equal to 1, which would have given zero for the absorbance at this energy. The remarkable conclusion is that the resulting spectra contain both the absorption and plasmon resonances of the bulk material, as is evidenced by the expressions Eqs. (15) and (16).

The presence of a peak at the frequency of the bulk plasmon in an electromagnetic radiation was already pointed out in the context of the so-called transition radiation [41], measured when a fast electron beam impinges a thin slab of metal perpendicular to the interface [42,43]. In a debate to interpret this effect with surface plasma oscillations [44,45], Economou [25] suggested on the contrary that the peak in the transition radiation was due to a switch from conditions of total reflection at $\omega = \omega_p$, the bulk plasma frequency, to conditions of large transmission at the neighboring points, and he supported his explanation by noting that the reflection coefficient of Airy [Eq. (13)] had a maximum at the bulk plasmon frequency. Since the experiment was done on metal, no peak at $\omega = \omega_0$ was observed. In our case, $R = |r_p^{\text{slab}}|^2$ (Fig. 8, continuous black line) is far from the total reflection. A non-negligible part of the beam is transmitted and absorbed.

It can be seen by considering the case of a biaxial material [46], where the dielectric functions for the parallel and perpendicular excitations are explicitly taken into account in the reflection and transmission coefficients [Eqs. (D10) and (D11)], that the resonance at $\omega_0$ is associated with the in-plane

dielectric function $(\varepsilon_{\parallel} - 1)$ and the structure at $\omega_p$ to the out-of-plane term $(1/\varepsilon_{\perp} - 1)$. To explain the peak at $\omega_p$ arising from the multiple reflection and transmission processes of the electromagnetic field on the interfaces of a thin slab, we suggest that the perpendicular component of the electric field could create a collective oscillation of the electronic density, which resonates at the plasma frequency.

The major result of our work is that the longitudinal-longitudinal contraction of the macroscopic dielectric tensor $\varepsilon_M^{\text{LL}}(\mathbf{q} \to \mathbf{0}; \omega)$ for a thin slab contains two structures that arise from the absorption and the plasmon peaks of the bulk counterpart, depending of the direction of the perturbation (Fig. 7). These two structures are the resonances of the response function of the system to the macroscopic part of the total potential, calculated within TD-DFT, by solving the Dyson equation [Eq. (2) or (10)]. These bulk features enter in the reflectance and transmittance coefficients of a thin slab. As a consequence, the TD-DFT, even if it is a longitudinal formalism, allows one to calculate the quantities involved in the optical response, which is by nature the response to a transverse field. Such a result was already pointed out in Ref. [47].

The structure coming from the plasmon peak arises in the so-called minimum vacuum configuration, in reciprocal space (Fig. 2, left), or in mixed-space (Fig. 6), which allows one to correctly include the interfaces (Fig. 5, right). This effect is similar to Airy's formula. In the classical approach, the dielectric layer is embedded between two media with a dielectric function equal to 1 (vacuum), and in the ab initio framework, the fact that we let the response function reach zero accounts for the vacuum environment around the slab. This effect, arising from the inclusion of the local field effects, reveals capable to mimic the multiple reflections and transmissions experienced by the electromagnetic field.

The silicon slab turns out to be a well-suited model system since the resonances for the absorption process and collective excitations are largely separated. However, the materials of interest for technology involving quasi-2D objects are the so-called traditional 2D materials, namely piles of layers bonded through van der Waals interactions.

## V. FOUR-GRAPHENE-LAYER STACKING

We consider an AB stacking of four layers of graphene [16]. The bulk counterpart is the AB graphite, for which the lattice parameter perpendicular to the hexagonal planes was taken as $c = 12.588$ bohrs. The imaginary parts of the macroscopic dielectric function ($\text{Im}[\varepsilon_M(\mathbf{q})]$), giving the absorption for vanishing $\mathbf{q}$) and of the macroscopic inverse dielectric function ($\text{Im}[\varepsilon_{00}^{-1}(\mathbf{q})]$, leading to the loss function) for the in-plane and out-of-plane components, calculated within TD-DFT, are plotted as a reminder in Fig. 9. $\text{Im}[\varepsilon_M(\mathbf{q}_{\parallel})]$ (continuous green line) presents two large peaks between 0 and 5 eV coming from the $\pi$ transitions and a wide resonance between 10 and 25 eV corresponding to the $\pi + \sigma$ ones. Due to the anisotropy of graphite, $\text{Im}[\varepsilon_M(q_z)]$ (dashed red line) differs: the $\pi$ contribution almost vanishes and the $\pi + \sigma$ one is composed of a first structure peaked at 11–12 eV and a second one at 16 eV. The in-plane loss function (double-dashed-dotted orange line, see the inset for an extended energy
<figure><img src="image_12.png" /><figcaption>FIG. 9. Imaginary parts of $\varepsilon_M(\mathbf{q})$ and $\varepsilon_{00}^{-1}(\mathbf{q})$ for the in- and out-of-plane excitations of bulk graphite: absorption spectrum in-plane (green continuous line) and out-of-plane (red dashed line); loss function in-plane (orange double-dashed-dotted line) (large energy range in the inset) and out-of-plane (blue dashed-dotted line).</figcaption></figure>

range) exhibits the well-known $\pi$ and $\pi + \sigma$ plasmons at 7 and 28 eV, respectively. The out-of-plane spectrum (dashed-dotted blue line) shows a vanishing $\pi$ contribution like the absorption, and it presents a resonance at 18 eV.

Based on the usual definition of the thickness of one graphene layer [32,33], we get for the thickness defined by the atomic positions $L_z^{\text{atom}} = 4 \times 6.294 = 25.176$ bohrs. Including the interfaces, which we evaluated to ~12 bohrs, we define $L_z^{\text{mat}} = 6 \times 6.294 = 37.764$ bohrs. We checked for one graphene layer that such an extension of the interfaces allows us to reproduce correctly the spectra for the different directions of the perturbation by applying the procedure of Sec. III C.

We applied to this quasi-2D object the reciprocal space formalism using the slab potential with the Selected-G scheme (see Sec. II). Im[$\varepsilon_M^{\text{LL}}(\mathbf{q}_y + \mathbf{q}_z)$] is shown in Fig. 10 as a plain-circle brown line for the object of thickness $L_z^{\text{mat}} = L_z^{\text{atom}} = 25.176$ bohrs, corresponding to the so-called zero vacuum case, or describing an object without interfaces, and as a continuous cyan line for the thickness $L_z^{\text{mat}} = 37.764$ bohrs, describing the object including the interfaces (minimum vacuum case). In Fig. 10, we also added {Im[$\varepsilon_M(\mathbf{q}_\|)$] + Im[$\varepsilon_M(\mathbf{q}_z)$]}/2 (dashed red line) and {Im[$\varepsilon_M(\mathbf{q}_\|)$] - Im[$\varepsilon_{00}^{-1}(\mathbf{q}_z)$]}/2 (blue dashed-dotted line), where the spectra for bulk graphite are the ones of Fig. 9.

The $\pi$ resonance is not well-suited to determine if the slab spectrum comes from the linear combination of the in-plane and out-of-plane macroscopic dielectric functions of the bulk graphite (dashed red line), or the linear combination of the in-plane macroscopic dielectric function and the out-of-plane macroscopic inverse dielectric function of the bulk graphite (dashed-dotted blue line), since the out-of-plane components are negligible compared to the in-plane ones. The difference of the amplitudes only comes from the “thickness” scaling factor. On the contrary, the $\pi + \sigma$ resonance, due to the presence of the $\sigma$ bonds, which contribute to the width of

<figure><img src="image_13.png" /><figcaption>FIG. 10. Comparison of the longitudinal-longitudinal contraction of the dielectric tensor of the four-layer-graphene object for the coupled excitation with the linear combinations of the different resonances of the bulk graphite (Fig. 9). Without interfaces ($L_z^{\text{mat}} = L_z^{\text{atom}} = 25.176$ bohrs): plain-circle brown line. With interfaces ($L_z^{\text{mat}} = 37.764$ bohrs): continuous cyan line. {Im[$\varepsilon_M(\mathbf{q}_\|$)] + Im[$\varepsilon_M(\mathbf{q}_z)$]}/2: dashed red line. {Im[$\varepsilon_M(\mathbf{q}_\|$)] - Im[$\varepsilon_{00}^{-1}(\mathbf{q}_z)$]}/2: dashed-dotted blue line.</figcaption></figure>

the interfaces, is more appropriate to get insight on the definition of the thickness of the quasi-2D object. The spectra of Fig. 9 behave exactly as those for the silicon slab (see Fig. 2). Without interfaces, the coupled component of the slab spectrum $\varepsilon_M^{\text{LL}}(\mathbf{q}_y + \mathbf{q}_z)$ is very well reproduced by the sum of the in- and out-of-plane absorption spectra of the bulk graphite (plain-circle brown versus dashed red line), exhibiting a bulk behavior. Including the interfaces, the coupled component of the quasi-2D spectrum $\varepsilon_M^{\text{LL}}(\mathbf{q}_y + \mathbf{q}_z)$ is in very good agreement with the sum of the in-plane absorption and out-of-plane plasmon of the bulk counterpart (continuous cyan versus dashed-dotted blue lines). One notices nevertheless a small discrepancy at the level of the peak at 18 eV arising from Im[$\varepsilon_{00}^{-1}(\mathbf{q}_z)$], which appears at 17 eV for the slab spectrum. Due to the $sp^2$ bonding in the plane of the layers, the $\pi$ and $\pi + \sigma$ resonances of the absorption and the loss function coexist in the same energy region (0–10 and 10–25 eV, respectively) (Fig. 9), so the identification of the respective contributions is less intuitive than for the silicon slab, but the same mixing occurs.

Finally, in order to verify that the longitudinal-longitudinal contraction of the dielectric tensor of the four-layer-graphene system could be related to the reflectance or transmittance features of an electromagnetic field impinging on this quasi-2D object, we calculated these quantities using Eqs. (D8) and (D9) derived from the transfer-matrix formalism [46] for the p-polarization in the biaxial case [18]. We used for the anisotropic dielectric functions the ones of graphite shown in Fig. 9. The results are plotted in Fig. 11.

The spectrum Im[$\varepsilon_M^{\text{LL}}(\mathbf{q}_y + \mathbf{q}_z)$] for the object including the interfaces (continuous cyan line) and the sum of the in-plane absorption and the out-of-plane plasmon of the bulk graphite (dashed blue line) have been multiplied by the energy ($\omega$) to
<figure><img src="image_14.png" /><figcaption>FIG. 11. Spectra of reflectance $R = |r_{pp}|^2$ with $r_{pp}$ given by Eq. (D8) (up-triangle violet line), $(1-T)$ where the transmittance $T = |t_{pp}|^2$ with $t_{pp}$ given by Eq. (D9) (down-triangle magenta line), compared with $\text{Im}[\varepsilon_M^{LL}(\mathbf{q}_y + \mathbf{q}_z)]$ for the slab with interfaces (continuous cyan line) and $\{\text{Im}[\varepsilon_M(\mathbf{q}_y)] - \text{Im}[\varepsilon_{00}^{-1}(\mathbf{q}_z)]\}/2$ (blue dashed-dotted line). The bulk spectra are the ones of Fig. 9. The last two spectra have been scaled by the energy. $R$ and $(1-T)$ have been scaled to fit the amplitude.</figcaption></figure>

account for the $k_0 = \omega/c$ dependence of the reflectance and transmittance. The reflectance and transmittance spectra have been scaled to fit the amplitude. It turns out that the transmittance spectrum (down-triangle magenta line) reproduces very well the spectral shape of the sum of the in-plane absorption and the out-of-plane plasmon of the bulk graphite (dashed-dotted blue line), as could be expected from the expression of the transmission coefficient [Eq. (D11)]. The coupled component $\varepsilon_M^{LL}(\mathbf{q}_y + \mathbf{q}_z)$ calculated for the slab including interfaces (continuous cyan line) is in very good agreement with the reflectance spectrum (up-triangle violet line), even if the depression at 17 eV is not present. The shift in the structure at 17-18 eV seems to result from the more complex expression of the reflection coefficient with the bulk dielectric functions [Eq. (D10)]. Further investigations are needed to establish the exact relationship between the TD-DFT response function and the reflectance, transmittance, or absorbance spectra.

Nevertheless, these results show that the scheme evidenced for the silicon slab also applies to traditional 2D materials. The interfaces must be included in the definition of the thickness of the slab to correctly account for the optical properties of the quasi-2D object. The LL contraction of the dielectric tensor, as it results from the TD-DFT, is a mixture of the in-plane absorption and the out-of-plane plasmon of the bulk counterpart, in a way similar to the transverse response.

## VI. CONCLUSION

In this paper, we studied the optical properties of a quasi-2D object within the *ab initio* TD-DFT formalism, including the local field effects. We proposed calculating the macroscopic average of the dielectric tensor of a quasi-2D object directly from a response function that contains only micro-

scopic excitations in order to avoid the use of the Adler and Wiser formula.

Comparing the result of the Dyson equation without and with the long-range contribution of the Coulomb potential, we evidenced that they are almost similar. We identified a small shift of the resonance at $\omega_0$ for the in-plane component, depending on the amplitude of the in-plane wave vector of the perturbation. For out-of-plane excitation they are identical, indicating that the long-range part of the Coulomb potential is always cut by the slab. The most spectacular result is that this resonance arises at the plasma frequency. From the reciprocal space calculations, using the Selected-G method to remove the spurious vacuum effect and the slab potential to isolate the slab, we evidenced that the appearance of the peak at the plasma frequency for the out-of-plane component of the response function to the macroscopic part of the total potential $[\tilde{\chi}_{00}(q_z)]$ occurs only when the interfaces are included in the 2D object. If the thickness of the slab is defined according to the atomic positions, called the zero vacuum case, and which corresponds to the bulklike definition of the unit cell, the Selected-G method allows one to build an object exhibiting a resonance of $\varepsilon_M^{LL}(\mathbf{q}_z)$ at $\omega_0$, as is expected in the bulk material. On the contrary, extending the thickness of the slab slightly in order to incorporate the interfaces, we obtained a resonance of $\varepsilon_M^{LL}(\mathbf{q}_z)$ at $\omega_p$, evidencing a complete change in the nature of the response for the $\mathbf{q}_z$ excitation.

To understand the origin of this unexpected excitation, we investigated the optical response of a thin slab by means of the classical electromagnetism formalism (Airy's formula, with a model Lorentz dielectric function). We showed that for the p-polarization of the light for a nonzero incident angle, the reflectance, transmittance, and absorbance spectra exhibit two peaks at $\omega_0$ and $\omega_p$. The expansion of the biaxial version of Airy's formula in the limit of the thin slab allows us to evidence that the two resonances appear in conjunction with $(\varepsilon_{\parallel} - 1)$ and $(1/\varepsilon_{\perp} - 1)$, which confirms the link of the out-of-plane resonance with the bulk plasmon. This optical response of a 2D object is by nature the response to a transverse field. The structure at $\omega_p$ cannot result from an absorption process, but it arises from the reflection and transmission processes of the interfaces. We suggest that the back and forth reflections of the component of the electric field perpendicular to the slab create an oscillation of the electronic density, which naturally resonates at the plasma frequency. This leads us to the conclusion that the thickness of the quasi-2D object must include the interfaces.

The remarkable conclusion that for an isolated quasi-2D object the in-plane response function presents a resonance at the frequency of the absorption, and the out-of-plane response function presents a peak governed by the plasma frequency, has the consequence that the Ehrenreich result, demonstrated for infinite cubic crystals, cannot apply anymore for 2D objects. For these highly anisotropic systems, the longitudinal-longitudinal contraction of the macroscopic dielectric tensor is equal to the transverse-transverse one only for the in-plane direction, while for the out-of-plane direction it is closely related to the plasmon of the bulk counterpart. For this reason, we need to label with "LL" the quantities calculated within TD-DFT, which are in essence the longitudinal-longitudinal contraction of the macroscopic
dielectric tensor. Nevertheless, our work shows that the (transverse) optical response of a 2D object can still be calculated within the longitudinal formalism of the TD-DFT, as was the case for bulk material, even if the phenomena involved are different.

## ACKNOWLEDGMENTS

For this work, access was granted to the HPC/AI resources of IDRIS and TGCC under the allocation (Grant No. 0900544 during the years 2018–2021) made by GENCI.

## APPENDIX A: LINEAR COMBINATIONS OF $\bar{\chi}_{00}(\mathbf{q}; \omega)$

The linear combinations use the properties of the dielectric tensor $\overleftrightarrow{\epsilon}_M(\omega)$ (in the following, we omit $\omega$ for clarity),

$$ \overleftrightarrow{\epsilon}_M = \begin{pmatrix} \epsilon_M^{xx} & \epsilon_M^{xy} & \epsilon_M^{xz} \\ \epsilon_M^{yx} & \epsilon_M^{yy} & \epsilon_M^{yz} \\ \epsilon_M^{zx} & \epsilon_M^{zy} & \epsilon_M^{zz} \end{pmatrix} $$

and

$$ \epsilon_M^{12} = \hat{\mathbf{q}}_1 \epsilon_M \hat{\mathbf{q}}_2 \quad \text{with} \quad |\hat{\mathbf{q}}_i| = 1. $$

Let $\hat{\mathbf{q}}$ be a normalized vector,

$$ \hat{\mathbf{q}} = (\hat{q}_x, \hat{q}_y, \hat{q}_z) = \left( \frac{q_x}{\sqrt{q_x^2 + q_y^2 + q_z^2}}, \frac{q_y}{\sqrt{q_x^2 + q_y^2 + q_z^2}}, \frac{q_z}{\sqrt{q_x^2 + q_y^2 + q_z^2}} \right). $$

We then obtain (taking $q_y = 0$ for simplification)

$$ \begin{aligned} \epsilon_M(\hat{\mathbf{q}}_x + \hat{\mathbf{q}}_z) &= (\hat{q}_x \quad 0 \quad \hat{q}_z) \begin{pmatrix} \epsilon_M^{xx} & \epsilon_M^{xy} & \epsilon_M^{xz} \\ \epsilon_M^{yx} & \epsilon_M^{yy} & \epsilon_M^{yz} \\ \epsilon_M^{zx} & \epsilon_M^{zy} & \epsilon_M^{zz} \end{pmatrix} \begin{pmatrix} \hat{q}_x \\ 0 \\ \hat{q}_z \end{pmatrix} \\ &= \hat{q}_x^2 \epsilon_M^{xx} + \hat{q}_x \hat{q}_z \epsilon_M^{xz} + \hat{q}_z \hat{q}_x \epsilon_M^{zx} + \hat{q}_z^2 \epsilon_M^{zz}. \end{aligned} $$

In an equivalent manner, this yields

$$ \epsilon_M(\hat{\mathbf{q}}_x - \hat{\mathbf{q}}_z) = \hat{q}_x^2 \epsilon_M^{xx} - \hat{q}_x \hat{q}_z \epsilon_M^{xz} - \hat{q}_z \hat{q}_x \epsilon_M^{zx} + \hat{q}_z^2 \epsilon_M^{zz}, $$

which allows us to extract $\epsilon_M^{zz}$:

$$ \epsilon_M^{zz} = \frac{q_x^2 + q_z^2}{q_z^2} \frac{\epsilon_M(\hat{\mathbf{q}}_x + \hat{\mathbf{q}}_z) + \epsilon_M(\hat{\mathbf{q}}_x - \hat{\mathbf{q}}_z)}{2} - \frac{q_x^2}{q_z^2} \epsilon_M^{xx}. \quad (A1) $$

We obtain immediately that $\epsilon_M^{yx} \equiv \epsilon_M(\hat{\mathbf{q}}_x)$ and $\epsilon_M^{zz} \equiv \epsilon_M(\hat{\mathbf{q}}_z)$. Using now the definition of $\epsilon_M(\hat{\mathbf{q}})$:

$$ \epsilon_M(\hat{\mathbf{q}}) = \lim_{\hat{\mathbf{q}} \to 0} \left[ 1 - \frac{4\pi}{|\hat{\mathbf{q}}|^2} \bar{\chi}_{00}(\hat{\mathbf{q}}) \right], \quad (A2) $$

and introducing Eq. (A2) into Eq. (A1), one finally obtains

$$ \bar{\chi}_{00}(\mathbf{q}_z) = \frac{\bar{\chi}_{00}(\mathbf{q}_x + \mathbf{q}_z) + \bar{\chi}_{00}(\mathbf{q}_x - \mathbf{q}_z)}{2} - \bar{\chi}_{00}(\mathbf{q}_x). \quad (A3) $$

One can either use Eq. (A1) for $\epsilon_M$ spectra, or Eq. (A3) for $\bar{\chi}$. To use the linear combination (A3), the $(\mathbf{q}_x)$ used to calculate $\bar{\chi}_{00}(\mathbf{q}_x)$ must be the same as the one in $(\mathbf{q}_x \pm \mathbf{q}_z)$, and the $(\mathbf{q}_z)$ for which we have calculated $\bar{\chi}_{00}(\mathbf{q}_z)$ is the one in $(\mathbf{q}_x \pm \mathbf{q}_z)$.

## APPENDIX B: TECHNICAL DETAILS OF THE FOURIER TRANSFORM AND TESTS OF CONSISTENCY IN MIXED-SPACE

The inverse Fourier transform (IFT) of $\chi_{G_z G'_z}^0(\mathbf{q}_\|)$ is performed using a 2D FFTW3 algorithm, where even $N$ integer values of $G_z = G_m$ (in units of $2\pi/L_z$, where $L_z$ is the size of the supercell) are taken into account. It reads [Eq. (5)]

$$ \chi^0(\mathbf{q}_\|, z_n, z_{n'}) = \frac{1}{L_z} \sum_{m, m'} e^{iG_m z_n} \chi_{G_m G_{m'}}^0(\mathbf{q}_\|) e^{-iG_{m'} z_{n'}}. $$

The result is an $N \times N$ matrix for $z_n = \frac{n}{2N}L_z$ with $n \in [-N, N-1]$, thus $\chi^0(\mathbf{q}_\|, z_n, z_{n'})$ extends in a box of size $L_z \times L_z$ and it can be checked that with such a normalization, it is independent of the height of the supercell used in reciprocal space.

To recover in reciprocal space the response function for the coupled $(\mathbf{q}_\| \pm \mathbf{q}_z)$ perturbation, we introduce the $\mathbf{q}_z$ direction of the perturbation by multiplying $\chi^0(\mathbf{q}_\|, z_n, z_{n'})$ by $e^{-iq_z(z_n - z_{n'})}$ in the Fourier transform (FT) process, leading to

$$ \chi_{G_m G_{m'}}^0(\mathbf{q}_\|, q_z) = \frac{L_z}{N^2} \sum_{n, n'} e^{-i(G_m + q_z)z_n} \chi^0(\mathbf{q}_\|, z_n, z_{n'}) e^{i(G_{m'} + q_z)z_{n'}}. \quad (B1) $$

The macroscopic average [21] is defined with the same methodology, yielding

$$ \begin{aligned} \chi_M^0(\mathbf{q}_\|, q_z) &= \frac{\Delta z_n^2}{L_z} \sum_{n, n'} e^{-iq_z z_n} \chi^0(\mathbf{q}_\|, z_n, z_{n'}) e^{iq_z z_{n'}} \\ &= \frac{1}{L_z N^2} \sum_{n, n'} e^{-iq_z z_n} \chi^0(\mathbf{q}_\|, z_n, z_{n'}) e^{iq_z z_{n'}} \\ &= \chi_{00}^0(\mathbf{q}_\|, q_z). \end{aligned} $$

We recover the FT [Eq. (B1)] for $G_m = G_{m'} = 0$. The sum over $z_n$ is done over the full size of the simulation box defined by $L_z$.

To describe a quasi-2D object, the wave vector used to define the direction of the excitation takes a very small value ($|\mathbf{q}_\|$ = $1 \times 10^{-4}$ a.u. and $|\mathbf{q}_z| = 1 \times 10^{-4}$ a.u. here, but we checked that $1 \times 10^{-5}$ a.u. or $1 \times 10^{-3}$ a.u. gave the same results). Since the head and the wings of $\chi_{G_z G'_z}^0(\mathbf{q}_\|)$ are proportional, respectively, to $|\mathbf{q}_\|$ and $|\mathbf{q}_z|$, they are quantities much smaller than the elements of the body. To ensure that we did not lose any information in the inverse Fourier transform leading to $\chi^0(\mathbf{q}_\|, z_n, z_{n'})$, we did some tests of consistency. The results are presented in Fig. 12.

We first checked that the Fourier transform algorithm [Eq. (B1)] allows us to recover the original $\chi^0$ calculated in reciprocal space. The left panel of Fig. 12 shows the comparison of the imaginary part of the FT of $\chi^0(\mathbf{q}_x, z_n, z_{n'})$ [Eq. (B1) for $G_m = G_{m'} = 0$] for $\mathbf{q}_x$ (violet square) and for $(\mathbf{q}_x + \mathbf{q}_z)$ (green down-triangle) with the $\chi_{00}^0(\mathbf{q})$ calculated in reciprocal space ($\chi_{\text{orig}}^0(\mathbf{q})$ for $\mathbf{q}_x$ (cyan circle) and $(\mathbf{q}_x + \mathbf{q}_z)$ (red up-triangle): both in-plane and coupled components $(\mathbf{q}_x + \mathbf{q}_z)$ are recovered.

To rule out the possibility that numerical errors could be compensated for by the reversibility of the FFT procedure,
<figure><img src="image_15.png" /><figcaption>FIG. 12. Left: Imaginary part of FT of $\chi^0(\mathbf{q}_x, z_n, z_{n'})$ [Eq. (B1) for $G_m = G_{m'} = 0$], calculated for $q_z = 0$ (violet square) or $q_z > 0$ (green down-triangle), compared to the original (reciprocal space) $\chi_{00}^0(\mathbf{q})$ (labeled $\chi_{\text{orig}}^0$) calculated for the polarizations $\mathbf{q}_x$ (cyan circle) and $(\mathbf{q}_x + \mathbf{q}_z)$ (red up-triangle). Center: Imaginary part of the macroscopic average [Eq. (6)] of $\chi^0(\mathbf{q}_x, z_n, z_{n'})$ normalized by $-4\pi/|\mathbf{q}|^2$, compared to NLF spectra in reciprocal space for $\mathbf{q}_x$ (violet square vs cyan circle) or $(\mathbf{q}_x + \mathbf{q}_z)$ (green down-triangle vs red up-triangle). Right: $\mathbf{q}_z$ component extracted from the linear combination of spectra calculated for $\mathbf{q}_x$, $\mathbf{q}_x + \mathbf{q}_z$ and $\mathbf{q}_x - \mathbf{q}_z$, in mixed-space (magenta diamond), compared with the $\mathbf{q}_z$ spectrum calculated directly in reciprocal space (turquoise square). The calculations are done for $||\mathbf{q}_x|| = ||\mathbf{q}_z|| = 1 \times 10^{-4}$ a.u., but the results are independent of the value of $|\mathbf{q}|$ in the optical limit, except for $\chi_{00}^0(\mathbf{q})$, which is proportional to $|\mathbf{q}|^2$. In Eqs. (B1) and (6), the summation is done on a box of size $L_z = 82.104$ bohrs.</figcaption></figure>

we also compared the spectra resulting from the macroscopic average [Eq. (6)] to the absorption spectra calculated without local fields (NLF) in reciprocal space, which are defined according to $\text{Im}[\varepsilon_M^{\text{NLF}}(\mathbf{q}; \omega)] = \text{Im}[1 - 4\pi/|\mathbf{q}|^2 \chi_{00}^0(\mathbf{q}; \omega)]$ in the optical limit $\mathbf{q} \to \mathbf{0}$. The center panel shows the imaginary part of the macroscopic average normalized by $-4\pi/|\mathbf{q}|^2$ for $\mathbf{q}_x$ (violet square) and $(\mathbf{q}_x + \mathbf{q}_z)$ (green down-triangle) compared to the absorption spectra $\text{Im}[\varepsilon_M^{\text{NLF}}(\mathbf{q}; \omega)]$ for $\mathbf{q}_x$ (cyan circle) and $(\mathbf{q}_x + \mathbf{q}_z)$ (red up-triangle). All the corresponding spectra are on top of each other. The difference in the amplitudes of the corresponding spectra between the left and center panels comes from the $1/|\mathbf{q}|^2$ factor, since in our calculations $|\mathbf{q}_x| = 1 \times 10^{-4}$ a.u. and $|\mathbf{q}_x + \mathbf{q}_z| = \sqrt{2} \times 10^{-4}$ a.u.

These results show that $\chi^0(\mathbf{q}_\parallel, z_n, z_{n'})$ does not suffer from the relative amplitude of the head and wings as compared to the body of the initial reciprocal space, $\chi_{G_z G'_z}^0(\mathbf{q}_\parallel)$. They also

show that the $\mathbf{q}_z$ dependence of $\chi_{00}^0(\mathbf{q})$ arose correctly from the multiplication of the $(z_n, z_{n'})$ element with the phase factor $e^{-iq_z(z_n-z_{n'})}$.

One can extract $\varepsilon_M^{\text{NLF}}(\mathbf{q}_z; \omega)$ using the properties of the dielectric tensor $\overleftrightarrow{\varepsilon}_M(\omega)$ (see Appendix A). Using the linear combination given by Eq. (A1), we obtain the spectrum $\text{Im}[\varepsilon_M^{\text{NLF}}(\mathbf{q}_z; \omega)]$ (magenta diamond) in Fig. 12, right: it is identical to the spectrum calculated in reciprocal space directly for the $\mathbf{q}_z$ component (turquoise square).

## APPENDIX C: $\tilde{\chi}$ IN MIXED-SPACE

For three-dimensional crystals, the response function associated with the absorption spectrum is the response function to the total macroscopic potential ($\tilde{\chi}$) solution of the reciprocal
space Dyson equation [Eq. (8)]:

$$ \bar{\chi}_{GG'}(\mathbf{q}) = \chi_{GG'}^0(\mathbf{q}) + \sum_{\mathbf{G}_1} \chi_{GG_1}^0(\mathbf{q}) \bar{v}_{\mathbf{G}_1}(\mathbf{q}) \bar{\chi}_{\mathbf{G}_1 G'}(\mathbf{q}), $$

where $\bar{v}_G(\mathbf{q}) = 0$ for $G = 0$ and

$$ \bar{v}_G(\mathbf{q}) \equiv v_G(\mathbf{q}) = \frac{4\pi}{||\mathbf{G} + \mathbf{q}||^2} \quad \text{for } \mathbf{G} \neq \mathbf{0}. $$

This approach, which aims to suppress the long-range component of the Coulomb potential $v_0$, is well suited for reciprocal space, since the $G$ elements can be addressed in a distinguished manner. The equivalent scheme in real space is more tricky. To go further, we separate the different equations for the head, wings, and body of the $\bar{\chi}_{GG'}$, and we isolate the term that involved the long-range component of the Coulomb potential:

$$ \bar{\chi}_{00} = \chi_{00}^0 + \chi_{00}^0 \bar{v}_0 \bar{\chi}_{00} + \sum_{G_1 \neq 0} \chi_{0G_1}^0 \bar{v}_{G_1} \bar{\chi}_{G_1 0}, $$

$$ \bar{\chi}_{0G'} = \chi_{0G'}^0 + \chi_{00}^0 \bar{v}_0 \bar{\chi}_{0G'} + \sum_{G_1 \neq 0} \chi_{0G_1}^0 \bar{v}_{G_1} \bar{\chi}_{G_1 G'}, $$

$$ \bar{\chi}_{G0} = \chi_{G0}^0 + \chi_{G0}^0 \bar{v}_0 \bar{\chi}_{00} + \sum_{G_1 \neq 0} \chi_{GG_1}^0 \bar{v}_{G_1} \bar{\chi}_{G_1 0}, $$

$$ \bar{\chi}_{GG'} = \chi_{GG'}^0 + \chi_{G0}^0 \bar{v}_0 \bar{\chi}_{0G'} + \sum_{G_1 \neq 0} \chi_{GG_1}^0 \bar{v}_{G_1} \bar{\chi}_{G_1 G'}, $$

where $G \neq 0$ and $G' \neq 0$.

By definition, $\bar{v}_0 = 0$, so the second term in each right member is suppressed. One sees that to obtain the Dyson equation for $\bar{\chi}$, it is equivalent to suppressing $v_0$ in the equation for $\chi$, or to suppressing the first column of the matrix $\chi^0$: $\chi_{G0}^0 \equiv 0 \forall G$. Indeed, these elements do not intervene in the other terms of the equations. It is thus equivalent to define a new matrix $\bar{\chi}^0$ identical to $\chi^0$, except for the first column, which is set to zero, and to solve

$$ \bar{\chi}_{GG'}(\mathbf{q}) = \chi_{GG'}^0(\mathbf{q}) + \sum_{\mathbf{G}_1} \bar{\chi}_{GG_1}^0(\mathbf{q}) v_{G_1}(\mathbf{q}) \bar{\chi}_{G_1 G'}(\mathbf{q}), \quad (C1) $$

$$ \text{where } v_G(\mathbf{q}) = \frac{4\pi}{||\mathbf{G} + \mathbf{q}||^2} \quad \forall \mathbf{G}. \quad (C2) $$

In reciprocal space, this is equivalent to solving Eq. (8) or Eq. (C1). In real space, it is not possible to suppress the long-range part of the Coulomb potential, but since our $\chi^0(z, z')$ is obtained by a FT of a reciprocal space matrix, one can define $\bar{\chi}_{GG'}^0(\mathbf{q})$, perform the inverse Fourier transform, and solve Eq. (C1) in mixed-space. The resulting equation is

$$ \bar{\chi}(\mathbf{q}_{\parallel}, z, z') = \chi^0(\mathbf{q}_{\parallel}, z, z') + \int dz_1 \int dz_2 \bar{\chi}^0(\mathbf{q}_{\parallel}, z, z_1) \times v(\mathbf{q}_{\parallel}, z_1, z_2) \bar{\chi}(\mathbf{q}_{\parallel}, z_2, z') $$

with $v(\mathbf{q}_{\parallel}, z, z') = \frac{2\pi}{|\mathbf{q}_{\parallel}|} e^{-|\mathbf{q}_{\parallel}|z-z'}$.

<figure><img src="image_16.png" /><figcaption>FIG. 13. Imaginary part of FT [Eq. (B1) for $G_m = G_{m'} = 0$] of $\chi(\mathbf{q}_x, z, z')$ and $\bar{\chi}(\mathbf{q}_x, z, z')$, respectively, after resolution of the Dyson equations [Eqs. (4) and (10)]. Zoom on the peak at 4 eV: spectra for $\mathbf{q}_x$ (black circle and blue up-triangle, respectively) and $(\mathbf{q}_z \pm \mathbf{q}_z')$ (red square/green diamond and magenta left-triangle/orange right-triangle, respectively). Calculation is for the slab where the atomic positions extend to 41.052 bohrs, introduced in a supercell of height 82.104 bohrs.</figcaption></figure>

The spectra obtained by solving the Dyson equation for $\bar{\chi}$ and $\chi$ in real space (with the FT [Eq. (B1) for $G_m = G_{m'} = 0$]), are presented in Fig. 13 for $\mathbf{q}_x$ and $(\mathbf{q}_x \pm \mathbf{q}_z)$. We zoomed on the peak at 4 eV. The three spectra $\bar{\chi}_{00}(\mathbf{q}_x)$ (blue up-triangle) and $\bar{\chi}_{00}(\mathbf{q}_x \pm \mathbf{q}_z)$ (magenta left-triangle and orange right-triangle, respectively) are identical. The same occurs for the three components of $\chi_{00}$ (black circle, red square, and green diamond, respectively). There is a tiny difference between $\bar{\chi}_{00}$ and $\chi_{00}$. It comes from the fact that the resonance of $\chi(\mathbf{q}_{\parallel})$ occurs at $\sqrt{\omega_0^2 + |\mathbf{q}_{\parallel}|L/2} \omega_p$, as we evidenced in Ref. [18], where $\omega_0$ is the absorption resonance and $\omega_p$ is the plasmon frequency (4 and 17 eV, respectively, for bulk silicon), and $L$ is the thickness of the slab, while the resonance of $\bar{\chi}_{00}(\mathbf{q}_{\parallel})$ is expected exactly at $\omega_0$. It can be checked numerically that $\chi_{00}(\mathbf{q})$ reaches $\bar{\chi}_{00}(\mathbf{q})$ as $|\mathbf{q}_{\parallel}| \to 0$. The difference between the two response functions illustrates the effect of the long-range component of the Coulomb potential for the perturbation parallel to the infinite directions: it leads to a small shift of the resonance, but in the limit of the ultrathin slab ($|\mathbf{q}_{\parallel}|L/2 \ll 1$), it does not have the capacity to move the resonance of $\chi_{00}(\mathbf{q}_{\parallel})$ up to the plasmon frequency. We remind the reader that it is a consequence of the huge reduction of the screening in the 2D object for the in-plane perturbation [18,22]. This result also confirms that the procedure we developed to remove the long-range component of the Coulomb potential in real space is adequate.

More puzzling is the fact that the amplitudes of the spectra for $\mathbf{q}_x$ and $(\mathbf{q}_x \pm \mathbf{q}_z)$ are the same. This was not the case for $\chi^0$, where the amplitude of $\chi_{00}^0(\mathbf{q}_{\parallel} \pm \mathbf{q}_z)$ was twice that of $\chi_{00}^0(\mathbf{q}_{\parallel})$ for the considered case $||\mathbf{q}_x|| = ||\mathbf{q}_z||$ (Fig. 12, left).
## APPENDIX D: ANALYSIS OF AIRY’S FORMULAS FOR A SLAB IN VACUUM

### 1. p-polarization

The Fresnel coefficient of reflection for a p-polarized beam with an angle of incidence $\theta$ is [Eq. (12)]

$$r_p = \frac{\varepsilon \cos \theta - \sqrt{\varepsilon - \sin^2 \theta}}{\varepsilon \cos \theta + \sqrt{\varepsilon - \sin^2 \theta}} \equiv \frac{A}{B} \equiv \frac{c-d}{c+d} \quad (D1)$$

The reflection coefficient for a slab in vacuum is [Eq. (13)]

$$r_p^{\text{slab}} = r_p \frac{1 - e^{2i\beta}}{1 - r_p^2 e^{2i\beta}} \\ \text{with } \beta = k_0 d \sqrt{\varepsilon - \sin^2 \theta} \quad \text{and} \quad k_0 = \frac{\omega}{c}.$$

$d$ is the thickness of the slab, $c$ is the speed of light in vacuum, and $\omega$ is the energy of the photon.

Considering the situation in which the thickness of the slab is much smaller than the wavelength, $\beta \ll 1$, one can expand $r_p^{\text{slab}}$ at first order in $\beta$:

$$r_p^{\text{slab}} = \frac{A}{B} \frac{-2i\beta}{1 - (\frac{A}{B})^2(1 + 2i\beta)} = \frac{-2i\beta(c^2 - d^2)}{4cd - 2i\beta(c^2 + d^2 - 2cd)},$$

which gives (see Fig. 14)

$$r_p^{\text{slab}} = \frac{-2ik_0 d \sqrt{\varepsilon - \sin^2 \theta} [\varepsilon^2 \cos^2 \theta - (\varepsilon - \sin^2 \theta)]}{4 \varepsilon \cos \theta \sqrt{\varepsilon - \sin^2 \theta} - 2ik_0 d \sqrt{\varepsilon - \sin^2 \theta} [\varepsilon^2 \cos^2 \theta + \varepsilon - \sin^2 \theta - 2 \varepsilon \cos \theta \sqrt{\varepsilon - \sin^2 \theta}]} \quad (D2)$$

$$= \frac{-2ik_0 d [\varepsilon \cos^2 \theta - 1 + 1/\varepsilon \sin^2 \theta]}{4 \cos \theta - 2ik_0 d [\varepsilon \cos^2 \theta + 1 - 1/\varepsilon \sin^2 \theta - 2 \cos \theta \sqrt{\varepsilon - \sin^2 \theta}]}$$

$$= \frac{-2ik_0 d [(\varepsilon - 1) \cos^2 \theta + (1/\varepsilon - 1) \sin^2 \theta]}{4 \cos \theta (1 + ik_0 d \sqrt{\varepsilon - \sin^2 \theta}) - 2ik_0 d [(\varepsilon - 1) \cos^2 \theta - (1/\varepsilon - 1) \sin^2 \theta] - 4ik_0 d \cos^2 \theta}$$

$$\approx \frac{-2ik_0 d [(\varepsilon - 1) \cos^2 \theta + (1/\varepsilon - 1) \sin^2 \theta]}{4 \cos \theta - 4ik_0 d \cos^2 \theta - 2ik_0 d [(\varepsilon - 1) \cos^2 \theta - (1/\varepsilon - 1) \sin^2 \theta]} \quad (D3)$$

$$\approx \frac{-2ik_0 d [(\varepsilon - 1) \cos^2 \theta + (1/\varepsilon - 1) \sin^2 \theta]}{4 \cos \theta - 2ik_0 d [(\varepsilon - 1) \cos^2 \theta - (1/\varepsilon - 1) \sin^2 \theta]} \quad (D4)$$

The transmission coefficient of p-polarized light for a slab in vacuum is [Eq. (14)]

$$t_p^{\text{slab}} = \frac{(1 - r_p^2) e^{i\beta}}{1 - r_p^2 e^{2i\beta}}.$$

<figure><img src="image_17.png" /><figcaption>FIG. 14. Results of different approximations of Airy’s formulas ($\theta = 45^\circ$). Reflectance = $|r_p^{\text{slab}}|^2$ (left) [transmittance = $|t_p^{\text{slab}}|^2$ (right)]: (black) Eq. (D2) [(D5)], (magenta) Eq. (D3) [(D6)], (orange) Eq. (D4) [(D7)].</figcaption></figure>
Using the quantities defined according to Eq. (D1), we obtain

$$t_p^{\text{slab}} = \frac{[1 - (\frac{A}{B})^2](1 + i\beta)}{1 - (\frac{A}{B})^2(1 + 2i\beta)} = \frac{4 c d (1 + i\beta)}{4 c d - 2i\beta (c^2 + d^2 - 2 c d)},$$

$$t_p^{\text{slab}} = \frac{4 \epsilon \cos \theta \sqrt{\epsilon - \sin^2 \theta} (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta})}{4 \epsilon \cos \theta \sqrt{\epsilon - \sin^2 \theta} - 2ik_0 d \sqrt{\epsilon - \sin^2 \theta} [\epsilon^2 \cos^2 \theta + \epsilon - \sin^2 \theta - 2 \epsilon \cos \theta \sqrt{\epsilon - \sin^2 \theta}]}, \quad (D5)$$

$$= \frac{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta})}{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta}) - 2ik_0 d [\epsilon \cos^2 \theta + 1 - 1/\epsilon \sin^2 \theta]}$$

$$= \frac{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta})}{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta}) - 2ik_0 d [(\epsilon - 1) \cos^2 \theta - (1/\epsilon - 1) \sin^2 \theta] - 4 ik_0 d \cos^2 \theta}$$

$$\approx \frac{4 \cos \theta}{4 \cos \theta - 4 ik_0 d \cos^2 \theta - 2ik_0 d [(\epsilon - 1) \cos^2 \theta - (1/\epsilon - 1) \sin^2 \theta]} \quad (D6)$$

$$\approx \frac{4 \cos \theta}{4 \cos \theta - 2ik_0 d [(\epsilon - 1) \cos^2 \theta - (1/\epsilon - 1) \sin^2 \theta]} \quad (D7)$$

This decomposition of $r_p^{\text{slab}}$ and $t_p^{\text{slab}}$ allows one to understand the reason why the reflectance and the transmittance will exhibit two peaks (see Fig. 14), one coming from $(\epsilon - 1)$ corresponding to the resonance of the absorption frequency, and one coming from $(1/\epsilon - 1)$, associated with the plasmon frequency.

## 2. s-polarization

The Fresnel coefficient of reflection for an s-polarized beam with an angle of incidence $\theta$ is

$$r_s = \frac{\cos \theta - \sqrt{\epsilon - \sin^2 \theta}}{\cos \theta + \sqrt{\epsilon - \sin^2 \theta}}.$$

The same expansion as for the p case gives

$$r_s^{\text{slab}} = \frac{2ik_0 d (\epsilon - 1)}{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta}) - 4ik_0 d \cos^2 \theta - 2ik_0 d (\epsilon - 1)} \\ \approx \frac{2ik_0 d (\epsilon - 1)}{4 \cos \theta - 4ik_0 d \cos^2 \theta - 2ik_0 d (\epsilon - 1)} \approx \frac{2ik_0 d (\epsilon - 1)}{4 \cos \theta - 2ik_0 d (\epsilon - 1)},$$

$$t_s^{\text{slab}} = \frac{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta})}{4 \cos \theta (1 + ik_0 d \sqrt{\epsilon - \sin^2 \theta}) - 4ik_0 d \cos^2 \theta - 2ik_0 d (\epsilon - 1)} \\ \approx \frac{4 \cos \theta}{4 \cos \theta - 4ik_0 d \cos^2 \theta - 2ik_0 d (\epsilon - 1)} \approx \frac{4 \cos \theta}{4 \cos \theta - 2ik_0 d (\epsilon - 1)}.$$

It is obvious that for the s-polarized light, no peak at the plasmon frequency should appear.

## 3. p-polarization for a biaxial slab in vacuum

In Ref. [18], based on the work of Ref. [46], we have established that for a slab of biaxial material in vacuum, the reflection and transmission coefficients are expressed as

$$r_{pp} = \frac{i \sin(\kappa_p) [1 - 1/\epsilon_{\perp} \sin^2(\theta) - \epsilon_{\parallel} \cos^2(\theta)]}{2 \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} \cos(\kappa_p) \cos(\theta) - i \sin(\kappa_p) [1 - 1/\epsilon_{\perp} \sin^2(\theta) + \epsilon_{\parallel} \cos^2(\theta)]}, \quad (D8)$$

$$t_{pp} = \frac{2 \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} \cos(\theta)}{2 \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} \cos(\kappa_p) \cos(\theta) - i \sin(\kappa_p) [[1 - 1/\epsilon_{\perp} \sin^2(\theta)] + \epsilon_{\parallel} \cos^2(\theta)]}, \quad (D9)$$

with $\kappa_p = k_0 d \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]}$, and $\epsilon_{\parallel}$ and $\epsilon_{\perp}$ are the in-plane and out-of-plane dielectric functions of the biaxial material. Considering the case in which the thickness of the slab is much smaller than the wavelength of the electromagnetic field, one has

$$\sin(\kappa_p) \approx \kappa_p = k_0 d \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]}, \\ \cos(\kappa_p) \approx 1,$$
and the reflection and transmission coefficients of a biaxial slab in vacuum become

$$
\begin{align*}
r_{pp} &= \frac{-2 i k_0 d \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} [\epsilon_{\parallel} \cos^2(\theta) - [1 - 1/\epsilon_{\perp} \sin^2(\theta)]]}{4 \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} \cos(\theta) - 2 i k_0 d \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} [[1 - 1/\epsilon_{\perp} \sin^2(\theta)] + \epsilon_{\parallel} \cos^2(\theta)]} \\
&= \frac{-2 i k_0 d [[(\epsilon_{\parallel} - 1) \cos^2(\theta) + (1/\epsilon_{\perp} - 1) \sin^2(\theta)]]}{4 \cos(\theta) - 4 i k_0 d \cos^2(\theta) - 2 i k_0 d [(\epsilon_{\parallel} - 1) \cos^2(\theta) - (1/\epsilon_{\perp} - 1) \sin^2(\theta)]} \\
&\approx \frac{-2 i k_0 d [[(\epsilon_{\parallel} - 1) \cos^2(\theta) + (1/\epsilon_{\perp} - 1) \sin^2(\theta)]]}{4 \cos(\theta) - 2 i k_0 d [(\epsilon_{\parallel} - 1) \cos^2(\theta) - (1/\epsilon_{\perp} - 1) \sin^2(\theta)]} \tag{D10}
\end{align*}
$$

$$
\begin{align*}
t_{pp} &= \frac{2 \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} \cos(\theta)}{2 \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} \cos(\theta) - i k_0 d \sqrt{\epsilon_{\parallel}} \sqrt{1 - [1/\epsilon_{\perp} \sin^2(\theta)]} [[1 - 1/\epsilon_{\perp} \sin^2(\theta)] + \epsilon_{\parallel} \cos^2(\theta)]} \\
&= \frac{4 \cos(\theta)}{4 \cos(\theta) - 4 i k_0 d \cos^2(\theta) - 2 i k_0 d [(\epsilon_{\parallel} - 1) \cos^2(\theta) - (1/\epsilon_{\perp} - 1) \sin^2(\theta)]} \\
&\approx \frac{4 \cos(\theta)}{4 \cos(\theta) - 2 i k_0 d [(\epsilon_{\parallel} - 1) \cos^2(\theta) - (1/\epsilon_{\perp} - 1) \sin^2(\theta)]} \tag{D11}
\end{align*}
$$

This decomposition of $t_{pp}$ and $r_{pp}$ gives expressions similar to those obtained from Airy's formula [Eqs. (D2) and (D5)] where we treated the isotropic case, but it shows that, for the biaxial case, the resonance at the absorption energy comes from the in-plane dielectric function, while the resonance at the plasmon frequency comes from the inverse of the out-of-plane dielectric function. These expressions explain why the spectra of Ref. [18] (Figs. 3 and 5), considered to be anisotropic dielectric functions of effective slabs of different thicknesses $d$, give the same reflectance and absorbance spectra (Ref. [18], Fig. 9). Indeed, since these anisotropic dielectric functions result from a mechanism that leads to a scaling with $1/d$ of $(\epsilon_{\parallel} - 1)$ and $(1/\epsilon_{\perp} - 1)$, the expressions [Eqs. (D10) and (D11)] show why an exact compensation occurs: the dielectric functions are involved through a linear combination of $d$ $(\epsilon_{\parallel} - 1)$ and $d$ $(1/\epsilon_{\perp} - 1)$, leading to a cancellation of the scaling effect [22].

[1] E. C. Dreaden, A. M. Alkilany, X. Huang, C. J. Murphy, and M. A. El-Sayed, *Chem. Soc. Rev.* **41**, 2740 (2012).

[2] S. Öğüt, J. R. Chelikowsky, and S. G. Louie, *Phys. Rev. Lett.* **79**, 1770 (1997).

[3] C. Delerue, G. Allan, and M. Lannoo, *Phys. Rev. Lett.* **90**, 076803 (2003).

[4] L. Wirtz, A. Marini, and A. Rubio, *Phys. Rev. Lett.* **96**, 126104 (2006).

[5] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, *Phys. Rev. Lett.* **105**, 136805 (2010).

[6] G. Onida, L. Reining, and A. Rubio, *Rev. Mod. Phys.* **74**, 601 (2002).

[7] E. Runge and E. K. U. Gross, *Phys. Rev. Lett.* **52**, 997 (1984).

[8] E. K. U. Gross and W. Kohn, *Phys. Rev. Lett.* **55**, 2850 (1985).

[9] R. van Leeuwen, *Phys. Rev. Lett.* **82**, 3863 (1999).

[10] L. Reining, V. Olevano, A. Rubio, and G. Onida, *Phys. Rev. Lett.* **88**, 066404 (2002).

[11] D. Aspes, *Thin Solid Films* **89**, 249 (1982).

[12] F. Sottile, F. Bruneval, A. G. Marinopoulos, L. K. Dash, S. Botti, V. Olevano, N. Vast, A. Rubio, and L. Reining, *Int. J. Quantum Chem.* **102**, 684 (2005).

[13] N. Tancogne-Dejean, C. Giorgetti, and V. Véniard, *Phys. Rev. B* **92**, 245308 (2015).

[14] C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, *Phys. Rev. B* **73**, 205119 (2006).

[15] S. Ismail-Beigi, *Phys. Rev. B* **73**, 233103 (2006).

[16] C. Giorgetti, I. Iagupov, and V. Véniard, *Phys. Rev. B* **101**, 035431 (2020).

[17] N. Tancogne-Dejean, C. Giorgetti, and V. Véniard, *Phys. Rev. B* **94**, 125301 (2016).

[18] S. Mazzei and C. Giorgetti, *Phys. Rev. B* **106**, 035431 (2022).

[19] S. L. Adler, *Phys. Rev.* **126**, 413 (1962).

[20] N. Wiser, *Phys. Rev.* **129**, 62 (1963).

[21] H. Ehrenreich, *The Optical Properties of Solids* (Academic, New York, 1965), p. 106.

[22] S. Mazzei, Second Harmonic Generation from silicon surfaces functionalized with DNA nucleobases: An *ab initio* description, Ph.D. thesis, Ecole polytechnique, 2021, https://etsf.polytechnique.fr/system/files/105736_MAZZEI_2021_archivage-2b_ETSF.pdf.

[23] C. Giorgetti and V. Véniard (unpublished).

[24] It is cut in the (001) direction, with a surface reconstruction (2 × 1) [see M. Palummo, G. Onida, R. Del Sole, and B. S. Mendoza, *Phys. Rev. B* **60**, 2522 (1999)]. Actually, in the present paper, the (x, y) axes have been rotated from 90°, with the consequence that the x component shown here is the direction perpendicular to the dimers and corresponds to the usual y.

[25] E. N. Economou, *Phys. Rev.* **182**, 539 (1969).

[26] V. M. Silkin, E. V. Chulkov, and P. M. Echenique, *Phys. Rev. Lett.* **93**, 176801 (2004).

[27] N. Tancogne-Dejean, *Ab initio* description of second-harmonic generation from crystal surfaces, Ph.D. thesis, Ecole Polytechnique, 2015.

[28] The dp code, Theoretical Spectroscopy Group-LSI-ETSF-Palaiseau; available at https://etsf.polytechnique.fr/Software/Ab_Initio.

[29] The abinit code, The ABINIT code is a common project of the Université Catholique de Louvain, Corning Incorporated, and other contributors, http://www.abinit.org.

[30] X. Gonze, F. Jollet, F. A. Araujo, D. Adams, B. Amadon, T. Applencourt, C. Audouze, J.-M. Beuken, J. Bieder, A.
Bokhanchuk, E. Bousquet, F. Bruneval, D. Caliste, M. Côté, F. Dahm, F. D. Pieve, M. Delaveau, M. D. Gennaro, B. Dorado, C. Espejo et al., *Comput. Phys. Commun.* **205**, 106 (2016).

[31] T. Tian, D. Scullion, D. Hughes, L. H. Li, C.-J. Shih, J. Coleman, M. Chhowalla, and E. J. G. Santos, *Nano Lett.* **20**, 841 (2020).

[32] Z. H. Ni, H. M. Wang, J. Kasim, H. M. Fan, T. Yu, Y. H. Wu, Y. P. Feng, and Z. X. Shen, *Nano Lett.* **7**, 2758 (2007).

[33] G. E. P. Dresselhaus and M. Dresselhaus, *Science of Fullerenes and Carbon Nanotubes Their Properties and Applications* (Elsevier, Amsterdam, 1996).

[34] In Ref. [16], $L_{\text{mat}}^*$ was chosen as an integer multiple (corresponding to the number of layers) of half of the distance between the atomic planes in the bulk counterpart material, as is usually done in the community. Nevertheless, the results published in Ref. [16] are still valid up to a normalization factor, since they only concern the in-plane component.

[35] J. D. Jackson, *Classical Electrodynamics*, 3rd ed. (John Wiley & Sons, New York, 2012).

[36] M. Born, E. Wolf, A. B. Bhatia, P. C. Clemmow, D. Gabor, A. R. Stokes, A. M. Taylor, P. A. Wayman, and W. L. Wilcock, *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light*, 7th ed. (Cambridge University Press, Cambridge, 1999).

[37] A.-J. Fresnel, *Memoire sur la loi des modifications que la reflexion imprime à la lumière polarisee* ("Memoir on the law of the modifications that reflection impresses on polarized light"), read 7 January 1823; reprinted in Fresnel, 1866, pp. 767-799 (full text, published 1831), pp. 753-762 (extract, published 1823). See especially, pp. 773 (sine law), 757 (tangent law), 760-761, and 792-796 (angles of total internal reflection for given phase differences), Académie des sciences (France); Mémoires de l'Académie des sciences (Firmin Didot Frères, Paris, 1823).

[38] G. Airy, *London, Edinburgh Dublin Philos. Mag. J. Sci.* **2**, 20 (1833).

[39] F. Abelès, *Ann. Phys.* **12**, 596 (1950).

[40] P. Yeh, *Optical Waves in Layered Media*, Wiley Series in Pure and Applied Optics (Wiley-Interscience, Hoboken, NJ, 2005).

[41] V. L. Ginzburg and I. M. Frank, *J. Expt. Theor. Phys.* (URSS) **16**, 15 (1946).

[42] W. Steinmann, *Phys. Rev. Lett.* **5**, 470 (1960).

[43] R. W. Brown, P. Wessel, and E. P. Trounson, *Phys. Rev. Lett.* **5**, 472 (1960).

[44] E. A. Stern, *Phys. Rev. Lett.* **8**, 7 (1962).

[45] R. A. Ferrell, *Phys. Rev.* **111**, 1214 (1958).

[46] M. Schubert, *Phys. Rev. B* **53**, 4265 (1996).

[47] R. Del Sole and E. Fiorino, *Phys. Rev. B* **29**, 4631 (1984).