Analysis of Resonant Soft X-ray Reflectivity of Anisotropic Layered Materials

Abstract

We present here a method for the quantitative prediction of the spectroscopic specular reflectivity line-shape in anisotropic layered media. The method is based on a 4 × 4 matrix formalism and on the simulation from the first principles (through density functional theory—DFT) of the anisotropic absorption cross-section. The approach was used to simulate the reflectivity at the oxygen K-edge of a 3,4,9,10-perylene-tetracarboxylic dianhydride (PTCDA) thin film on Au(111). The effect of film thickness, orientation of the molecules, and grazing incidence angle were considered. The simulation results were compared to the experiment, permitting us to derive information on the film geometry, thickness, and morphology, as well as the electronic structure.

1. Introduction

The study of surfaces and thin films often involves the use of visible, UV, and X-ray radiation as a probe. In this context, photon-in, photon-out techniques have the valuable advantage of being in general non-destructive and they can be applied either to conducting and non-conductive samples.

Non-resonant specular elastic scattering (reflectivity), in particular, is widely carried out in the hard X-ray range for the study of crystal structure and morphology, exploiting the electron density contrast among different atoms [1]. Reflectivity at resonance provides additional advantages in terms of atomic and depth-resolved investigation of the chemical, structural, and magnetic properties of a variety of systems, including polymeric and organic materials, liquid interfaces, spintronic systems, and inorganic layered materials [2,3]. Element selectivity together with high contrast between different structures is obtained. The sampling depth of reflectivity is not limited to the near-surface region, as for electron yield spectroscopies, but buried interfaces can be accessed as well (up to several tens of nanometers below the surface). Moreover, the technique is quantitative by nature and measurements are relatively simple, besides the fact that in the X-rays and at resonance they often require synchrotron radiation [4]. In contrast, reflectivity spectra can be extremely difficult to interpret without a suitable model. Several factors contribute to the energy-resolved reflectivity signal—these include the wavelength-dependent optical properties (optical constants) of the materials and the morphology of the system, including layer thickness and stacking, the quality of the interfaces, the geometrical arrangement of the constituents down to the atomic scale. These are indeed the unknowns that most likely the experiment would like to explore. The richness of information can definitely represent an obstacle for a rapid understanding of the reflectivity spectra.

We developed a protocol, as well as a related computer code, to obtain simultaneous quantitative information on the structure, interface morphology, chemical properties, magnetic properties, and optical properties with sub-nanometer depth resolution of layered materials that can present structural (lattice) and/or magnetic dichroism (in this work we restrict to “structural” dichroism). Our method can be applied in any photon energy range, both at resonant and non-resonant frequencies, and it assumes that the response is linear, that is the medium and its response are described in terms of a dielectric framework. Our method takes account materials with anisotropic optical properties, which can be determined by the specific orientation of intrinsic anisotropic crystalline layers, or by the anisotropic arrangement of microscopic constituents (molecules or molecular building blocks) inside layered structures, as for organic molecular films. The method is based on the quantitative prediction of the spectral line-shape also across specific elemental edges through (1) the numerical calculation from the first principles of the elements of the dielectric tensor of each material in a given layered system, (2) the simulation of the propagation of the electromagnetic field in the stack of layers and the computation of the optical response (i.e., specular reflectivity in the present case), and (3) the comparison (and fitting) of the simulation with the experiment.

In some previous works [5,6,7,8], the method has been successfully applied by us in the soft X-rays at resonance to study anisotropic thin and ultrathin organic films on metals and insulators to get quantitative information on chemical composition, layer-resolved molecular orientation, optical anisotropies, and overall thickness of the films, simultaneously. We addressed in particular some prototypical cases that are interesting for molecular electronics applications, focusing on reflectivity at the carbon K-edge: 1,4-benzenedimethantiol single layers on Au(111) [5], 3,4,9,10-perylene-tetracarboxylic dianhydride (PTCDA) on Au(111) [6], pentacene on SiO₂ [7], and tetracene single crystal [8].

In this work, we illustrate the details of our method to calculate the optical properties of anisotropic layered materials (OPAL) and to simulate the specular reflectivity. We include an application example, referred to the simulation of the reflectivity at the oxygen K-edge of a thin film of PTCDA on Au(111). In particular, we apply here our approach for the first time to the O 1s edge, since in our previous studies [6] we concentrated specifically on the edge of C 1s.

The method transfers the knowledge gained in the years in the field of visible and infrared spectroscopy of thin anisotropic films [9,10,11] to the soft-X range. Moreover, it offers a tool of quantitative investigation in those cases where electron spectroscopy cannot be applied because of the presence of nonconductive materials, poor vacuum, or materials buried at distances from the surfaces higher than the electron-free collision path.

The approach that we implemented to simulate the optical properties of anisotropic layers is of general validity. It can be applied in any wavelength range. It allows one to calculate the propagation of light through a layered medium composed of an unlimited number of anisotropic layers interposed between vacuum on one side and a semi-infinite substrate on the other side, as represented in Figure 1a. The procedure extends to anisotropic and absorbing media the method introduced in the fundamental works of Parratt [12], Bertrand et al. [13], Berreman [14], and Yeh [15,16].

First, the simulation method assumes that each layer in the stack is represented by its specific complex dielectric tensor. The entries of the dielectric tensor can be derived either by direct measurements on the anisotropic bulk materials (i.e., by means of absorption measurements) or through simulations. Parameterized tensors can be also used in order to fit the calculation to the experiment.

Second, each layer is also characterized by its thickness (a scalar d). This implicitly corresponds to assuming each layer to have sharp and perfectly flat surfaces. Note, however, that this hypothesis can be partially relaxed by introducing effective buffers between the constituent layers.

In the present case, we deal with organic molecular films. They can present macroscopic anisotropic optical properties if they are composed by molecules arranged in an ordered manner. This is indeed the case for organic molecular crystals, but even in those cases that a crystal is not formed, as for self-assembled layers (SAMs), a preferential orientation of the molecules is often obtained, due to the interaction of the molecules with the substrate and between themselves in each layer. An effective dielectric tensor can therefore be derived, on the basis of the specific molecular arrangement in each layer and the optical properties of the single isolated molecule, assuming that these are not affected substantially by molecule–molecule interactions.

2. Materials and Methods

2.1. Evaluation of the Organic Film Dielectric Tensor in Correspndence of Absorption Resonances

The first step for the construction of the tensor is the calculation of the absorption cross-section of a single isolated molecule. This can be achieved thanks to density functional theoretical calculations (DFT), by computing electric dipole transitions along 3 mutually perpendicular directions that correspond to the main molecular axes. To calculate the dipole transition moments and oscillator strengths of a PTCDA-free molecule (Figure 1b), we used the code StoBe [17], applying the transition potential method [18]. The evaluation of the cross-sections along the molecular axes at the O K-edge follows closely what has been described in [6] for the C K-edge. The calculated absorption cross-sections are related to the imaginary part of the molecular scattering factors f₂(ℏω) along the 3 axes, which are linked, well outside the edge, to the Henke tabulated values [19]. In doing this, the sum is considered of all contributions to f₂(ℏω) of all the atomic constituents of the molecule (for PTCDA 24 C atoms, 6 O atoms, 8 H atoms), up to ℏω = 30,000 eV. The real part of the molecular scattering factors f₁ is then evaluated through Kramers–Krönig (KK) integral transformations [6]. In the present case, this was achieved in practice by first constructing the odd functions f_{2_odd}(ℏω), where f_{2_odd}(−ℏω) = −f₂(ℏω) and f_{2_odd}(ℏω) = f₂(ℏω), and then by calculating the Hilbert transforms of f_{2_odd}(ℏω) on the domain from −30,000 to 30,000 eV, in steps of 0.1 eV. This procedure can be implemented straightforward in most software for data analysis (including Mathematica [20], Matlab [21], IgorPro [22]). The complex refractive index along the three principal axes is derived from f₁(ℏω) and f₂(ℏω):

$${\tilde{n}}_j=n_{1,j}+i n_{2,j}=1-\frac{r_e}{2\pi }{\lambda }^{2}N\left(f_{1,j}-i f_{2,j}\right)$$

(1)

where j = x_m,y_m,z_m denotes the molecular axis index, r_e is the classical electron radius, λ is the light wavelength, and N is the density of the molecules. For PTCDA, we used a value of N = 1.37 × 10²¹ molecules/cm³ [6].

The complex dielectric constants along the main molecular axes are finally obtained:

$${\tilde{\epsilon }}_j={\left({\tilde{n}}_j\right)}^2={\epsilon }_{1,j}+i {\epsilon }_{2,j}=\left(n_{1,j}^2-n_{2,j}^2\right)+i 2n_{1,j}{n}_{2,j}$$

(2)

The three dielectric constants are used to fabricate the dielectric tensor of the organic film, taking into account the possible geometrical organization of the molecules in the local environment around each molecule. In doing this, we neglect the interaction between the molecules, and we assume that the macroscopic optical properties are de facto determined by the optical properties of the single molecules and by their geometrical arrangement. This is a reasonable assumption in the soft X-ray regime, where absorption edges are mainly determined by the molecular properties and to a much less extent by solid-state effects. In the present case, it can be assumed that molecules organize in domains, in which molecules show parallel z_m axes and an alternate herringbone arrangement in the x_my_m plane, with adjacent molecules laying nearly perpendicular to each other [23,24,25,26]. This can be described by a diagonal dielectric tensor of a film in which the xx and yy elements are taken as the azimuthal angular average (in the x_my_m plane) of ${\tilde{\epsilon }}_{x_m}$ and ${\tilde{\epsilon }}_{y_m}$ and the zz element is given by ${\tilde{\epsilon }}_{z_m}$. This model is further supported by the fact that the beam footprint (higher than 100 μm²) in a reflectivity experiment is much larger than the domains and the experimental signal integrates over a large number of domains that are typically oriented following the threefold symmetry of the Au (111) substrate. In fact, on the (111) surface of Au, the PTCDA molecules are known to organize in equivalent threefold domains, oriented at 120° with respect to each other [23,24,25,27,28]. This results in an overall in-plane isotropy of the reflectivity, assuming that the scattered intensity from the different domains is added up incoherently. Therefore, the dielectric tensor of a PTCDA layer can be written as

$${\overleftrightarrow{\epsilon }}_{PTCDA}=\left(\begin{array}{ccc}{\tilde{\epsilon }}_{xx}& 0& 0\\ 0& {\tilde{\epsilon }}_{yy}& 0\\ 0& 0& {\tilde{\epsilon }}_{zz}\end{array}\right),\mathrm{with}{\tilde{\epsilon }}_{xx}={\tilde{\epsilon }}_{yy}$$

(3)

However, it is worth stressing that the method and the software we developed also permit one to handle dielectric tensors with non-vanishing off-diagonal elements (such as magnetic media, chiral, or distorted systems).

2.2. Calculation of Light Propagation in the Layered Medium

The system is treated as a sequence of layers labelled l, with l varying from 1 to n, each described by a dielectric tensor ${\overleftrightarrow{\epsilon }}_l$ and by a thickness d_l. Layers are separated by flat planar and parallel interfaces. To describe a possible reorientation of the axes of each layer with respect to the x-, y-, and z-axes of the sample frame of reference (Figure 1a), we applied coordinate rotation matrices R_l in terms of the Euler angles ϕ_l, θ_l, ψ_l (describing a sequence of rotations around the z-x-z-axes) to the dielectric tensor so that the dielectric tensor becomes ${\overleftrightarrow{\epsilon }}_{Rl}=R_l\left({\varphi }_l,{\theta }_l,{\psi }_l\right){\overleftrightarrow{\epsilon }}_{Rl}{R}_l^{-1}\left({\varphi }_l,{\theta }_l,{\psi }_l\right)$. Since R_l is orthogonal, the dielectric tensor in x-, y-, and z-coordinates must be symmetric (for non-magnetic systems). Eventually, if multiple domains with different orientations are present in each layer, an average over the possible orientation angles is also applied.

For each layer, the electric field is assumed to propagate as a plane wave ${\overrightarrow{E}}_l{e}^{i\left({\overrightarrow{k}}_l·\overrightarrow{r}-\omega t\right)}$, and the wave equation in the momentum-frequency space yields [3,29,30]

$${\overrightarrow{k}}_l\times \left({\overrightarrow{k}}_l\times {\overrightarrow{E}}_l\right)+{\omega }^2{\mu }_0{\epsilon }_0{\overleftrightarrow{\epsilon }}_{Rl}{\overrightarrow{E}}_l={\overrightarrow{k}}_l\left({\overrightarrow{k}}_l·{\overrightarrow{E}}_l\right)-k_l^2{\overrightarrow{E}}_l+{\omega }^2{\mu }_0{\epsilon }_0{\overleftrightarrow{\epsilon }}_{Rl}{\overrightarrow{E}}_l=0$$

(4)

In the present case, the impinging wavevector ${\overrightarrow{k}}_0$ (with $k_0^2={\omega }^2{\mu }_0{\epsilon }_0$) is taken in the xz plane, so that, in the absence of roughness, assuming homogeneous layers in the xy plane with step-like potentials along z at each interface, in each layer ${\overrightarrow{k}}_l=\left(k_{lx},0,k_{lz}\right)$. Moreover, $k_{lx}$ is conserved in the propagation through the layers and $k_{lx}=k_{0x}$. In matrix form, the wave equation then becomes

$$\left(\begin{array}{ccc}{k}_0^2{\epsilon }_{R{l}_{xx}}-k_{lz}^2& k_0^2{\epsilon }_{R{l}_{xy}}& k_0^2{\epsilon }_{R{l}_{xz}}+k_{0x}{k}_{lz}\\ k_0^2{\epsilon }_{R{l}_{yx}}& k_0^2{\epsilon }_{R{l}_{yy}}-k_{0x}^2-k_{lz}^2& k_0^2{\epsilon }_{R{l}_{yz}}\\ k_0^2{\epsilon }_{R{l}_{zx}}+k_{lz}{k}_{0x}& k_0^2{\epsilon }_{R{l}_{zy}}& k_0^2{\epsilon }_{R{l}_{zz}}-k_{0x}^2\end{array}\right)\left(\begin{array}{c}{E}_{lx}\\ E_{ly}\\ E_{lz}\end{array}\right)=0.$$

(5)

Non-zero solutions for ${\overrightarrow{E}}_l$ exist if the determinant of the coefficient matrix is zero. The equation

$$det\left(\begin{array}{ccc}{k}_0^2{\epsilon }_{R{l}_{xx}}-k_{lz}^2& k_0^2{\epsilon }_{R{l}_{xy}}& k_0^2{\epsilon }_{R{l}_{xz}}+k_{0x}{k}_{lz}\\ k_0^2{\epsilon }_{R{l}_{yx}}& k_0^2{\epsilon }_{R{l}_{yy}}-k_{0x}^2-k_{lz}^2& k_0^2{\epsilon }_{R{l}_{yz}}\\ k_0^2{\epsilon }_{R{l}_{zx}}+k_{lz}{k}_{0x}& k_0^2{\epsilon }_{R{l}_{zy}}& k_0^2{\epsilon }_{R{l}_{zz}}-k_{0x}^2\end{array}\right)=0$$

(6)

reduces to a 4th degree polynomial in $k_{lz}$ with complex coefficients. In practice, in our Matlab implementation, we directly computed the roots of this polynomial by the roots command, which calculates the roots of a polynomial by computing the eigenvalues of its companion matrix (see, e.g., https://it.mathworks.com/help/matlab/ref/roots.html). This leads to 4 complex solutions for the z-projection of the wavevector, $k_{l{z}_{\sigma }}$ ($\sigma =1,2,3,4)$. These are grouped into 2 pairs, for wavevectors directed downwards ($\sigma =1,2$) and upwards ($\sigma =3,4$).

The four complex solutions of $k_{l{z}_{\sigma }}$ are used, in turn, to solve relation (5), in order to obtain the 4 complex electric field unit vectors ${\widehat{u}}_{l_{\sigma }}$ associated to each of these wavevectors. This is because the amplitudes of the 4 fields are left undetermined by solving the homogeneous linear system represented by expression (5). Numerically, these electric field unit vectors, that is, the null-space of the coefficient matrix in (5), may be found applying singular value decomposition algorithms [20,21,22]. In this regard, we remark that we did not use the null command of Matlab, since roundoff errors may cause matrices to be considered nonsingular by the program (implying an empty null space). Instead, we explicitly employed the singular value decomposition by the svd command (which is based on the algorithms employed in LAPACK [31]) and appropriately selected the vector corresponding to the null space.

It should be noted that if the dielectric tensor is diagonal with equal components, which is the case for vacuum and isotropic media (i.e., the substrate), the $k_{l{z}_{\sigma }}$ solutions are coincident in pairs (for $\sigma =1,2$ and $\sigma =3,4$) and the null-space is bidimensional, corresponding to the plane perpendicular to the propagation direction ${\widehat{k}}_{l_{\sigma }}$. In this case, since any linear combination of the 2 degenerate propagation vectors is a solution, we are free to take ${\widehat{u}}_{l_1}$ and ${\widehat{u}}_{l_4}$ parallel to $\widehat{y}$ (perpendicular with the scattering plane) and ${\widehat{u}}_{l_2}$ and ${\widehat{u}}_{l_3}$ perpendicular to $\widehat{y}$ (parallel to the scattering plane). This is shown for vacuum in Figure 1a, where the 2 coincident incoming wavevectors relate to the 2 perpendicular ${\overrightarrow{E}}_{01}$ (s-polarization) and ${\overrightarrow{E}}_{02}$ (p-polarization) electric fields, while the 2 coincident outgoing wavevectors relate to ${\overrightarrow{E}}_{03}$ and ${\overrightarrow{E}}_{04}$, laying in the scattering plane and perpendicular to it, respectively.

The total electric field in layer l is then given by a linear combination of the 4 solutions ${\widehat{u}}_{l_{\sigma }}$, with

$${\overrightarrow{\mathcal{E}}}_l\left(\overrightarrow{r},t\right)=\left\{{\displaystyle \sum }_{\sigma =1}^4{C}_{l\sigma }{\widehat{u}}_{l\sigma }\exp \left[i{k}_{lz\sigma }\left(z-z_l\right)\right]\right\}\exp \left[i\left(k_{0x}x-\omega t\right)\right]$$

(7)

where $z_l$ represents the bottom level of layer l and $C_{l\sigma }$ are the field amplitudes at $z_l$ (with $z_{\left(l-1\right)}\ge z\ge z_l$). The $C_{l\sigma }$ amplitudes have to be determined by imposing the continuity conditions of the fields at the interfaces. Without considering the $\exp \left(i\left(k_{0x}x-\omega t\right)\right)$ factor, we can reduce expression (7) to

$${\overrightarrow{E}}_l={{\displaystyle \sum }}_{\sigma =1}^4{C}_{l\sigma }{\widehat{u}}_{l_{\sigma }}$$

(8)

at $z=z_l.$ This is valid for layers l = 0,…,n. The substrate is semi-infinite and therefore it is convenient to refer the amplitudes of the fields for the substrate to the interface with layer n, at $z=z_n$. In this case, the field in the substrate becomes

$${\overrightarrow{\mathcal{E}}}_{sub}\left(\overrightarrow{r},t\right)=\left\{{\displaystyle \sum }_{\sigma =1}^4{C}_{sub,\sigma }{\widehat{u}}_{sub,\sigma }\exp \left[i{k}_{sub,z\sigma }\left(z_n-z\right)\right]\right\}\exp \left[i\left(k_{0x}x-\omega t\right)\right]$$

(9)

with $z\le z_n$.

In turn, the magnetic field in layer l is given by

$${\overrightarrow{\mathcal{H}}}_l\left(\overrightarrow{r},t\right)=\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\frac{1}{\left|{\overrightarrow{k}}_0\right|}\left\{{\displaystyle \sum }_{\sigma =1}^4{\overrightarrow{k}}_{l\sigma }\times C_{l\sigma }{\widehat{u}}_{l\sigma }\exp \left[i{k}_{lz\sigma }\left(z-z_l\right)\right]\right\}\exp \left[i\left(k_{0x}x-\omega t\right)\right]$$

(10)

with ${\overrightarrow{k}}_{l_{\sigma }}=\left(k_{0x},0,k_{l{z}_{\sigma }}\right)$. At $z=z_l$ and without considering the factor $\exp \left(i\left(k_{0x}x-\omega t\right)\right)$, this is reduced to

$${\overrightarrow{H}}_l=\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\frac{1}{\left|{\overrightarrow{k}}_0\right|}{{\displaystyle \sum }}_{\sigma =1}^4\left({\overrightarrow{k}}_{l_{\sigma }}\times C_{l\sigma }{\widehat{u}}_{l_{\sigma }}\right)$$

(11)

The magnetic field in the substrate, following the same considerations leading to expression (9), is given by

$${\overrightarrow{\mathcal{H}}}_{sub}\left(\overrightarrow{r},t\right)=\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\frac{1}{\left|{\overrightarrow{k}}_0\right|}\left\{{\displaystyle \sum }_{\sigma =1}^4{\overrightarrow{k}}_{sub,\sigma }\times C_{sub,\sigma }{\widehat{u}}_{sub,\sigma }\exp \left[i{k}_{sub,z\sigma }\left(z_n-z\right)\right]\right\}\exp \left[i\left(k_{0x}x-\omega t\right)\right]$$

(12)

The continuity for the electric and magnetic fields ${\mathcal{E}}_{lx}, {\mathcal{E}}_{ly}, {\mathcal{H}}_{lx}, {\mathcal{H}}_{ly}$ between layers l and ($l-1$) yields

$$\begin{gathered} {{\displaystyle \sum }}_{\sigma =1}^4{C}_{\left(l-1\right)\sigma }{\widehat{u}}_{{\left(l-1\right)}_{\sigma }}·\widehat{x}={{\displaystyle \sum }}_{\sigma =1}^4{C}_{l\sigma }{\widehat{u}}_{l_{\sigma }}·\widehat{x}\exp \left({ik}_{{lz}_{\sigma }}{d}_l\right) \\ {{\displaystyle \sum }}_{\sigma =1}^4{C}_{\left(l-1\right)\sigma }{\widehat{u}}_{{\left(l-1\right)}_{\sigma }}·\widehat{y}={{\displaystyle \sum }}_{\sigma =1}^4{C}_{l\sigma }{\widehat{u}}_{l_{\sigma }}·\widehat{y}\exp \left({ik}_{{lz}_{\sigma }}{d}_l\right) \\ {{\displaystyle \sum }}_{\sigma =1}^4\left({\overrightarrow{k}}_{{\left(l-1\right)}_{\sigma }}\times C_{\left(l-1\right)\sigma }{\widehat{u}}_{{\left(l-1\right)}_{\sigma }}\right)·\widehat{x}={{\displaystyle \sum }}_{\sigma =1}^4\left({\overrightarrow{k}}_{l_{\sigma }}\times C_{l\sigma }{\widehat{u}}_{l_{\sigma }}\right)·\widehat{x}\exp \left({ik}_{{lz}_{\sigma }}{d}_l\right) \\ {{\displaystyle \sum }}_{\sigma =1}^4\left({\overrightarrow{k}}_{{\left(l-1\right)}_{\sigma }}\times C_{\left(l-1\right)\sigma }{\widehat{u}}_{{\left(l-1\right)}_{\sigma }}\right)·\widehat{y}={{\displaystyle \sum }}_{\sigma =1}^4\left({\overrightarrow{k}}_{l_{\sigma }}\times C_{l\sigma }{\widehat{u}}_{l_{\sigma }}\right)·\widehat{y}\exp \left({ik}_{{lz}_{\sigma }}{d}_l\right) \end{gathered}$$

(13)

Exponentials in expression (13) account for the propagation of the fields trough layer l, from top interface with layer ($l-1$) to the bottom interface of layer l. These relations can be compacted in a 4 × 4 projection matrix

$$P_l=\left(\begin{array}{cccc}{\widehat{u}}_{l_1}·\widehat{x}& {\widehat{u}}_{l_2}·\widehat{x}& {\widehat{u}}_{l_3}·\widehat{x}& {\widehat{u}}_{l_4}·\widehat{x}\\ {\widehat{u}}_{l_1}·\widehat{y}& {\widehat{u}}_{l_2}·\widehat{y}& {\widehat{u}}_{l_3}·\widehat{y}& {\widehat{u}}_{l_4}·\widehat{y}\\ \left({\overrightarrow{k}}_{l_1}\times {\widehat{u}}_{l_1}\right)·\widehat{x}& \left({\overrightarrow{k}}_{l_2}\times {\widehat{u}}_{l_2}\right)·\widehat{x}& \left({\overrightarrow{k}}_{l_3}\times {\widehat{u}}_{l_3}\right)·\widehat{x}& \left({\overrightarrow{k}}_{l_4}\times {\widehat{u}}_{l_4}\right)·\widehat{x}\\ \left({\overrightarrow{k}}_{l_1}\times {\widehat{u}}_{l_1}\right)·\widehat{y}& \left({\overrightarrow{k}}_{l_2}\times {\widehat{u}}_{l_2}\right)·\widehat{y}& \left({\overrightarrow{k}}_{l_3}\times {\widehat{u}}_{l_3}\right)·\widehat{y}& \left({\overrightarrow{k}}_{l_4}\times {\widehat{u}}_{l_4}\right)·\widehat{y}\end{array}\right)$$

(14)

and a 4 × 4 transport matrix

$$T_l=\left(\begin{array}{cccc}\exp \left(i{k}_{l{z}_1}{d}_l\right)& 0& 0& 0\\ 0& \exp \left(i{k}_{l{z}_2}{d}_l\right)& 0& 0\\ 0& 0& \exp \left(i{k}_{l{z}_3}{d}_l\right)& 0\\ 0& 0& 0& \exp \left(i{k}_{l{z}_4}{d}_l\right)\end{array}\right),$$

(15)

so that the electric field amplitudes can be expressed as four-component vectors ${\overrightarrow{C}}_l=\left(C_{l_1},C_{l_2},C_{l_3},C_{l_4}\right)$ and

$${\overrightarrow{C}}_{\left(l-1\right)}=P_{\left(l-1\right)}^{-1}{P}_l{T}_l{\overrightarrow{C}}_l$$

(16)

This approach can be iterated to relate the fields in successive layers. For a system composed of n layers followed by a semi-infinite substrate,

$$\begin{array}{c}{\overrightarrow{C}}_0=P_0^{-1}{P}_1{T}_1{P}_1^{-1}{P}_2{T}_2\dots P_{\left(n-1\right)}^{-1}{P}_n{T}_n{P}_n^{-1}{P}_{sub}{\overrightarrow{C}}_{sub}=M{\overrightarrow{C}}_{sub}\\ \mathrm{where} M=P_0^{-1}{P}_1{T}_1{P}_1^{-1}{P}_2{T}_2\dots P_{\left(n-1\right)}^{-1}{P}_n{T}_n{P}_n^{-1}{P}_{sub}\end{array}$$

(17)

It can be assumed that the electric fields coming upwards from the substrate are zero (Figure 1a), that is, $C_{su{b}_3}= C_{su{b}_4}=0$, and that the impinging fields at the vacuum side are known $(C_{0_1}$ and $C_{0_2}$, with ${\widehat{u}}_{0_1}$ and ${\widehat{u}}_{0_2}$ corresponding to two orthogonal unit vectors parallel to $\widehat{y}$ and perpendicular to it, respectively). Equation (17) results in a linear system of 4 equations in 4 unknowns, the reflected fields $C_{0_3}$ and $C_{0_4}$ and the transmitted fields $C_{su{b}_1}$ and $C_{su{b}_2}$. In this case the expression

$${\overrightarrow{C}}_0=M{\overrightarrow{C}}_{sub} \to \left(\begin{array}{c}{C}_{0_1}\\ C_{0_2}\\ C_{0_3}\\ C_{0_4}\end{array}\right)=\left(\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& M_{14}\\ M_{21}& M_{22}& M_{23}& M_{24}\\ M_{31}& M_{32}& M_{33}& M_{34}\\ M_{41}& M_{42}& M_{43}& M_{44}\end{array}\right) \left(\begin{array}{c}{C}_{su{b}_1}\\ C_{su{b}_2}\\ 0\\ 0\end{array}\right)$$

(18)

can be rewritten to obtain

$$\left(\begin{array}{c}{C}_{0_1}\\ C_{0_2}\\ 0\\ 0\end{array}\right)=\left(\begin{array}{cccc}{M}_{11}& M_{12}& 0& 0\\ M_{21}& M_{22}& 0& 0\\ M_{31}& M_{32}& -1& 0\\ M_{41}& M_{42}& 0& -1\end{array}\right) \left(\begin{array}{c}{C}_{su{b}_1}\\ C_{su{b}_2}\\ C_{0_3}\\ C_{0_4}\end{array}\right)$$

(19)

The linear system (19) can be solved numerically, for instance, by the linsolve command of Matlab, which, for square matrices, uses LU factorization with partial pivoting (see, e.g., https://it.mathworks.com/help/matlab/ref/linsolve.html). Once the 4 unknown field amplitudes are determined, reflectivity is obtained through

$$R=\frac{{\left|C_{0_3}\right|}^2+{\left|C_{0_4}\right|}^2}{{\left|C_{0_1}\right|}^2+{\left|C_{0_2}\right|}^2}$$

(20)

which yields for s-polarization and p-polarization reflectivity

$$R^s={\left(\frac{{\left|C_{0_3}\right|}^2+{\left|C_{0_4}\right|}^2}{{\left|C_{0_1}\right|}^2}\right)}_{C_{0_2}=0}{\mathrm{and}R}^p={\left(\frac{{\left|C_{0_3}\right|}^2+{\left|C_{0_4}\right|}^2}{{\left|C_{0_2}\right|}^2}\right)}_{C_{0_1}=0}$$

(21)

In analogy, in case that the substrate is vacuum, transmitted intensity from layer n yields

$$T=\frac{{\left|C_{{sub}_1}\right|}^2+{\left|C_{{sub}_2}\right|}^2}{{\left|C_{0_1}\right|}^2+{\left|C_{0_2}\right|}^2}$$

(22)

2.3. Experimental Methods

The X-ray absorption and reflectivity experiments were carried out at the IOM-CNR synchrotron BEAR beamline [32,33,34] at Elettra (Trieste, Italy). The Au(111) single crystal substrate (Mateck) was cleaned by repeated cycles of Ar+ sputtering (2 keV at grazing incidence for 10 min) and annealing (up to 450 °C for 20 min). Surface cleanliness was controlled by photoemission (by inspection of the Au valence band, C 1s and Au 4f core levels). The PTCDA molecules were evaporated in situ in ultra-high-vacuum (UHV), at a base pressure P = 1 × 10⁻⁹ mbar [6] from an home-made source, constituted by a tantalum pinhole source that was resistively heated at a temperature of 270 °C. An evaporation rate of 0.05–0.10 nm/min was used, which was monitored with a quartz microbalance. In particular, we focus here on a PTCDA film of 2 nm of nominal coverage, the same film studied in [6], to which the reader is addressed for further details. The specular reflectivity experiments were taken at a grazing incidence θ ≈ 8° in s- and p-polarization geometries. The photon energy was scanned across the O K-edge with a resolution of 0.2 eV, keeping fixed θ−2θ scattering conditions. The intensity of the direct (I₀) and reflected (I) beams was measured with a photodiode (IRD SXUV-100) in separate scans. Possible beam instabilities were taken into account by acquiring the photocurrent drained by the last optical element of the beamline at the same time during I₀ and I scans, just before the scattering chamber. The reflectivity was obtained by evaluating the ratio between I and I₀, after each of the two signals was corrected (normalized) for its respective beam instabilities. The degree of linear light polarization was selected of the order of 0.9, obtained by accepting photons only from the central portion of the beam from the bending magnet, with an angular acceptance of ± 0.1 mrad with respect to the plane of the synchrotron orbit. In simulations, we assumed completely linear polarized light, as expressed by Equation (21) above. A Si window filter (thickness 0.1 µm) was placed in the beam before incidence on the sample to reduce diffuse light contribution. X-ray absorption spectra (XAS) were acquired simultaneously, along with reflectivity, in total electron yield mode, by measuring the drain current from the sample. The XAS spectra were normalized to the beam flux, measuring the total electron yield from a carbon- and oxygen-free Au sample in the same energy region of the O K-edge.

3. Results and Discussion

Carbon K-edge of PTCDA on Au(111) was investigated by our group in [6], where reflectivity was simulated as a function of different parameters, as film thickness, molecular orientation, and scattering geometry, and simulations were compared to the experimental result. Here, we followed the same scheme, focusing instead on the O K-edge. This is the first time we performed the analysis at this specific edge. This is noteworthy, since the higher the photon energy, the lower the reflectivity signal in absolute units, and the comparison between simulation and experiment may become more critical.

In Figure 2a,b, the calculated diagonal elements of the molecular film dielectric tensor are shown. This refers to the idealized case of molecules perfectly aligned to the x-, y-, and z-axes of the substrate. In Figure 2c, the experimental near-edge X-ray absorption signal is also shown, as taken on a film of 2 nm of nominal coverage in s- and p-light scattering geometries. It can be seen that all the main features of the experimental spectra are reproduced by calculations of ${\epsilon }_2$.

The first sharp and asymmetric resonant structure, particularly pronounced in p-polarization in Figure 2c, was associated to 1s → π* transitions involving carbonyl oxygen atoms (C=O). The second feature at about 535 eV was also ascribed to 1s → π* transitions from the bridging O atoms. Features above 538 eV, mainly pronounced in s-scattering, were related to 1s → σ* transitions [35,36].

The experimental reflectivity spectra taken in s- and p-scattering are shown in Figure 3 for the same film of 2 nm of nominal coverage. The experimental curves are compared to the results of simulations obtained for PTCDA films of 2 nm of coverage, considering a flat lying arrangement, a herringbone configuration in the xy plane, and threefold-oriented domains, due to the (111) symmetry of the substrate. Slightly different angles of incidence of the radiation had to be considered for s- and p-scattering for the best matching between experiment and theory. These findings are compatible and consistent with the results presented previously by us at the C K-edge [6]. The effective thickness evaluated at the O K-edge in particular nicely matched the nominal thickness. It is noteworthy that in order to reproduce the experimental reflectivity, we did not have to consider any roughness, suggesting that the films were particularly flat.

It can be further noticed that the experimental spectrum in s-polarization presented a small but sizable feature at about 531 eV in Figure 3, while this was not present in the simulation. This is also consistent with s-polarization NEXAFS shown in Figure 2c. We associate this to some possible degree of disorder. In principle, this effect can also be partially ascribed to the non-complete linear polarization of the incoming light.

To show the sensitivity of reflectivity to the variation of the scattering and film parameters, as we did for the C K-edge [6], we report in Figure 4 the effect of film thickness (a), the effect of the grazing incidence angle (b), and the effect of the tilt angle of the molecules (c), from flat laying to almost perpendicular configuration. These results clearly show the high sensitivity of the technique to the different parameters. While the change of the incidence angle affected the overall intensity of reflectivity, both at the edge and outside, the variation of the film thickness principally influenced the intensity of the O K-edge features. The variation of the tilt angle of the molecules instead showed the progressive rising of π* features in s-polarization together with their progressive attenuation in p-polarization, while opposite behavior occurred for σ* structures. As was apparent also at the C K-edge [6], the sensitivity to the variation of the tilt angle was less pronounced with respect to the other parameters. If PTCDA molecules were not lying completely flat, tilts below 10° would be difficult to discriminate from the comparison with the experiment. In the present case in particular, we did not associate the oscillation of the experimental spectrum at 531 eV in s-polarization (Figure 3) to the effect of the tilt in the whole film, as this would also affect p-polarization with a drastic decrease of the reflectivity dip at this same energy, which was not observed. Instead, we tended to associate this to some level of disorder, with some number of molecules adopting a more perpendicular orientation, and/or to the ellipticity of the photon, with the incoming beam not being completely linear polarized. This effect was mostly evident in s-polarization, because for perfectly flat lying molecules and perfectly linear polarization, this spectral region would otherwise give a flat, featureless signal.

In the present study, the surface and interface roughness were not taken into account. Although we believe that this was not relevant in the present case, the surface and interlayer roughness can play and important role, especially for thicker films and systems composed of multiple layers. It may be expected that roughness can affect the overall intensity but not the overall line-shape and relative intensities. The effect of roughness will be considered in further evolutions of our code.

In perspective, the matrix elements of the different optical transitions that are used to construct the dielectric tensor could also be parameterized to improve the fitting with the experiment and to single out the spectral regions where the optical properties deviate from that of free, noninteracting molecules. However, it is remarkable that in spite of the simplifications introduced in the procedure and of the lack of adjustable parameters in the evaluation of the molecular cross-sections, the experimental features were nicely reproduced by the simulation, even in absolute magnitude.

4. Conclusions

Reflectivity is a powerful tool that can give rich information on layered media, once accompanied by suitable simulations. The 4 × 4 matrix method described here can be applied in any spectral range and it is specifically conceived for anisotropic media and in correspondence of spectral resonances. The method was used in this work to interpret the experimental line-shape of grazing incidence reflectivity experiments carried out at the O K-edge of an organic ultra-thin film of PTCDA on Au(111). The simulated line-shape is based on ab initio calculation of the optical absorption of the PTCDA molecules, and, therefore, is not based on entry parameters either taken form experiments or phenomenological. At this stage, the interaction between the molecules or with the substrate are not taken into account. The good agreement that can be reached with the experiment indicates that for molecular films solid state effects have a minor influence on the overall line-shape with respect to film morphology, molecular orientation, and film thickness. At a grazing incidence of the order of 8°, the reflectivity at the O K-edge was of the order of 10⁻², with oscillations in the 10⁻³ scale. In spite of this, high sub-nanometer sensitivity was obtained. The experimental spectra can be interpreted in terms of films of flat laying molecules, with the effective thickness obtained by simulations nicely corresponding with the nominal coverage. Small differences in the experimental spectra with respect to simulations were interpreted in terms of some level of disordered molecules, showing a more tilted alignment with respect to the rest of the film, and/or due to the effect of the non-complete linear polarization of the incoming photons.