Quantitative Predictions of Photoelectron Spectra in Amorphous Molecular Solids from Multiscale Quasiparticle Embedding

TL;DR

This paper introduces a multiscale first-principles simulation framework that accurately predicts ultraviolet photoelectron spectra of amorphous molecular solids, aiding in the precise determination of HOMO energies crucial for organic electronics.

Contribution

The novel multiscale approach combines deposition simulation, Green's Function Theory, and vibrational coupling to bridge experimental and theoretical spectral analysis of amorphous materials.

Findings

Excellent agreement with experimental spectra of MADN is achieved.

HOMO energies can be determined with an accuracy better than 0.1 eV.

The method provides molecular-level insights into spectral features.

Abstract

We present a first-principles-based multiscale simulation framework for quantitative predictions of the high-energy part of the Ultraviolet Photoelectron Spectroscopy (UPS) spectra of amorphous molecular solids. The approach combines a deposition simulation, many-body Green's Function Theory, polarizable film-embedding, and multimode electron-vibrational coupling and provides a molecular-level view on the interactions and processes giving rise to spectral features. This insight helps bridging the current gap between experimental UPS and theoretical models as accurate analyses are hampered by the energetic disorder, surface-sensitivity of the measurement and the complexity of excitation processes. In particular this is relevant for the unambiguous determination the highest occupied molecular orbital energy (HOMO) of organic semiconductors, a key quantity for tailoring and engineering new opto-electronic devices. We demonstrate the capabilities of the simulation approach studying the spectrum of two isomers of 2-methyl-9,10-bis(naphthalen-2-yl)anthracene (MADN) as archetypical materials showing a clearly separated HOMO peak in experiment. The agreement with experiment is excellent, suggesting that our approach provides a route for determining the HOMO energy with an accuracy better than 0.1eV.

Tables

Table 1. Table 1: Characteristics of the predicted UPS spectrum (see the caption of Fig.  6 ) for α 𝛼 \alpha -MADN and β 𝛽 \beta -MADN at the four different levels (0 − - 3) of the multiscale quasiparticle embedding procedure. The HOMO peak position, onset, and the full width at half maximum (FWHM) (all in eV) are compared to the respective experimental results. The table also gives the calculated bulk adiabatic ionization energy ε HOMO,bulk ad superscript subscript 𝜀 HOMO,bulk ad \varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}} .

	vacuum		$GW/\text{pMM}$
	KS	$GW$	Static	Polar.	Exp.
$𝜶$-MADN
${\epsilon }_{\text{HOMO,bulk}}^{\text{ad}}$				-5.91
UPS peak	-4.92	-6.78	-6.76	-6.10	-6.12
UPS onset	-4.71	-6.57	-6.44	-5.79	-5.79
UPS FWHM	0.25	0.25	0.39	0.39	0.37
$𝜷$-MADN
${\epsilon }_{\text{HOMO,bulk}}^{\text{ad}}$				-5.89
UPS peak	-4.80	-6.70	-6.67	-6.09	-6.14
UPS onset	-4.58	-6.48	-6.37	-5.75	-5.73
UPS FWHM	0.25	0.26	0.36	0.41	0.42

Equations

$$S_{\text{UPS}}(E)=\frac{1}{N_{\text{m}}}\sum_{j=1}^{N_{\textrm{m}}}S_{\text{el-vib}}(E;\varepsilon_{j})\exp\bigg{(}{-\frac{z_{0}(x,y)-z_{j}}{\Lambda}}\bigg{)},$$

$$[T_0+V_{\text{ext}}+V_{\text{Hartree}}+V_{\text{xc}}]|{\varphi }_i^{\text{KS}}⟩={ϵ}_i^{\text{KS}}|{\varphi }_i^{\text{KS}}⟩,$$

$$[H_0+\Sigma (E)]G_1(E)=E{G}_1(E),$$

$$\Sigma (\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\frac{i}{2\pi }\int d{\omega }^{\prime }{G}_1(\mathbf{r},{\mathbf{r}}^{\prime },\omega +{\omega }^{\prime })W(\mathbf{r},{\mathbf{r}}^{\prime },{\omega }^{\prime }).$$

$$H_{ij}^{GW}(E)={\epsilon }_i^{\text{KS}}{\delta }_{ij}+⟨{\varphi }_i^{\text{KS}}|\Sigma (E)-V_{\text{xc}}|{\varphi }_j^{\text{KS}}⟩.$$

$$\begin{array}{rl}{\epsilon }_i^{\text{QP},\mathrm{pert}}& ={\epsilon }_i^{\text{KS}}+\Delta {ϵ}_i^{GW}\\ & ={\epsilon }_i^{\text{KS}}+⟨{\varphi }_i^{\text{KS}}|\Sigma ({\epsilon }_i^{\text{QP},\mathrm{pert}})-V_{\text{xc}}|{\varphi }_i^{\text{KS}}⟩.\end{array}$$

$${\varphi }_i^{\text{KS}}(\mathbf{r})=\sum _{j=0}^M{X}_{ij}{\psi }_j(\mathbf{r}).$$

$$(ij|kl)=\iint d\mathbf{r}d{\mathbf{r}}^{\prime }{\psi }_i(\mathbf{r}){\psi }_j(\mathbf{r})v_c(\mathbf{r},{\mathbf{r}}^{\prime }){\psi }_k({\mathbf{r}}^{\prime }){\psi }_l({\mathbf{r}}^{\prime })$$

$$(ij|kl)\approx \sum _{\nu ,\mu }(ij|\nu )(\nu |\mu {)}^{-1}(\mu |kl).$$

$$\Delta Q_t^a=-\sum _{b\in B}\sum _{\begin{array}{c}{a}^{\prime }\in A\\ a^{\prime }\ne a\end{array}}{\alpha }_{t{t}^{\prime }}^{a{a}^{\prime }}{T}_{t^{\prime }u}^{a^{\prime }b}(Q_u^b+\Delta Q_u^b).$$

$$E_{\text{MM}}^{(s)}=\frac{1}{2}\sum _{a\in \mathcal{A}}\sum _{b\in \mathcal{B}}(Q_t^{a(s)}+\Delta Q_t^{a(s)})T_{tu}^{ab}{Q}_u^{b(s)}.$$

$$\text{SDOS}(E)=\frac{1}{N_{\text{m}}}\sum _j\delta (E-{\epsilon }_j(z_j))\exp \left(-\frac{z_0(x,y)-z_j}{\Lambda }\right).$$

$$S_{\text{el-vib}}^{\text{Marcus}}\big{(}\Delta E\big{)}=\frac{1}{\sqrt{4\pi\lambda_{\text{tot}}k_{B}T}}\exp\left(-\frac{(\Delta E-\lambda_{\text{tot}})^{2}}{4\lambda_{\text{tot}}k_{B}T}\right),$$

$$S_{\text{el-vib}}^{\text{FQ}}\big{(}\Delta E\big{)}=\frac{1}{2\pi\hbar}\int_{-\infty}^{+\infty}e^{i\frac{\Delta E}{\hbar}t}e^{-F(0)}e^{F(t)}\,dt,$$

Full text

I Introduction

Amorphous organic semiconductors are intensively applied in opto-electronic devices such as organic light-emitting diodes (OLEDs) Tang and VanSlyke (1987); Kido et al. (1995); Baldo et al. (1998); Adachi (2014), photovoltaic cells Halls et al. (1995); Yu et al. (1995); Zhao et al. (2017) and photodetectors Peumans et al. (2003); Jansen-van Vuuren et al. (2016); Kielar et al. (2016). Device properties can be tuned by varying the chemical building blocks or the material processing conditions, and by combining different molecular materials in complex blends or layer stacks. A key parameter determining the functioning of a material in a device is the ionization energy, often termed the highest occupied molecular orbital (HOMO) energy, $\varepsilon_{\textrm{HOMO}}$. Relative changes of $\varepsilon_{\textrm{HOMO}}$ on the order of 0.1 eV can already significantly alter the charge transport through host-guest materials or across internal interfaces between layers. However, even when using the perhaps most direct method for measuring $\varepsilon_{\textrm{HOMO}}$, Ultraviolet Photoelectron Spectroscopy (UPS), this level of accuracy has so far not been accomplished.

Excitation processes in organic semiconductors are complex because their localized nature gives rise to strong structural reorganization (polaron formation) and electron-vibration coupling, which leads to shifts, broadening and additional features in the UPS spectrum Salaneck et al. (1980); Malagoli et al. (2004); Kera et al. (2009); Kera and Ueno (2015). Combined with the energetic disorder originating from the amorphous structure and the surface-sensitivity of the measurement, this obstructs the unambiguous analysis of the spectra. As a result, the method used for deducing $\varepsilon_{\textrm{HOMO}}$ from the spectra (from the first peak energy or from an effective onset energy?) is a subject of debate Hill et al. (2000); D’Andrade et al. (2005); Krause et al. (2008). This uncertainty hampers the use of UPS studies for the rational design of new devices, and the combined use of high-resolution UPS, inverse UPS and photoluminescence spectra for obtaining accurate exciton binding energies Yoshida et al. (2015). Qualitative understanding of the spectra is often sought via gas-phase single-molecule calculations based on density-functional theory (DFT). However, the obtained energy levels need to be artificially shifted and broadened due to the well-known underestimation of the single-particle energy gap by DFT Perdew (1985); Mori-Sánchez et al. (2008) and due to the effects of intermolecular interactions Tadayyon et al. (2004). Those calculations lack an explicit link to the molecular morphology, cannot resolve surface and bulk contributions to the density of states (DOS), and do not account for the spectral consequences of the molecular ionization process resulting from the excitation of molecular vibrations. This lack of predictive power combined with the ambiguity in extracting the DOS from the experimental data is a big obstacle for the development of layer stacks for organic photovoltaics or next-generation OLEDs, for which the functioning and ultimate performance is already sensitive to energy level variations of only 100 meV.

In this work we present a first-principles-based multiscale simulation approach that bridges the current gap between experimental UPS and theoretical models by providing a quantitative prediction of the high-energy part of the UPS spectrum from which the ionization potential is derived. It consists of an accurate evaluation of (i) quasiparticle energy levels within the $GW$ approximation and (ii) thin-film embedding effects, using a hybrid quantum-mechanics/molecular-mechanics (QM/MM) approach that takes the molecular polarizabilities and the long-range interactions due to partially ordered static multipole moments into account, (iii) the inclusion of surface sensitivity via the electron attenuation length (EAL), $\Lambda$, and (iv) a full-quantum treatment of electron-vibration coupling. We focus here on only one type of initial state (frontier orbital), and simulate the UPS spectrum for perpendicular emission as a weighted sum of individual molecular environment-dependent densities of states according to

[TABLE]

with $\varepsilon_{j}$ the frontier orbital energy level of molecule $j$, $z_{j}$ the distance of the molecule’s center-of-mass (COM) to the corrugated surface at $z_{0}(x,y)$ not , and $S_{\text{el-vib}}(E;\varepsilon_{j})$ the energy-dependent spectral shape due to electron-vibration coupling. $N_{\textrm{m}}$ is the total number of molecules included in the summation, which is equal to the number of molecules for which from a vapour deposition simulation the atomistic morphology is obtained (see Sec. II.1). In view of the large optical absorption depth, we neglect optical matrix element effects.

As prototypical systems, we study the UPS spectrum for thin films of the $\alpha$ and $\beta$ isomers of 2-methyl-9,10-bis(naphthalen-2-yl)anthracene (MADN), whose chemical structures are shown as insets in Fig. 1. MADN is a morphologically stable amorphous wide-gap semiconductor Lee et al. (2004) that is used extensively as an ambipolar host material in OLEDs containing deep blue fluorescent emitter molecules Lee et al. (2005); Liao et al. (2005); Ho et al. (2006); Huang et al. (2012); Ho et al. (2010); Ahn et al. (2018). The methyl substituent disrupts the symmetry and stabilizes the material against crystallization. The type of coupling of the anthracene core and the naphtyl substituents ($\alpha$ or $\beta$) affects the planarity of the molecules, and thereby the frontier orbital energies and their distribution in a thin film. We regard MADN as particularly suitable for this study because it exhibits a HOMO peak that originates from a single non-degenerate state, located predominantly on the anthracene core. Experimental high resolution (low instrumental broadening) UPS measurements (see Fig. 1 and Sect. S1 of the Supplemental Material sup ), show that the peak full width at half maximum ($\sim 0.4\,\mathrm{eV}$) is significantly smaller than for many other often-used hole transporting and emitting materials in OLEDs. It is furthermore advantageous that the HOMO state is well-separated from the deeper levels. The selection of the two isomers enables us to study the effects of morphology differences and the related effects on energy level shifts due to the small molecular dipole moments.

The paper is organized as follows: In Sec. II, we provide an overview of the methodology and the computational details, including the specifics of the simulation of the atomistic thin-film morphologies, the calculation of quasiparticle energies in the $GW$ approximation, the QM/MM quasiparticle embedding schemes, and the carrier-vibration coupling. Section III contains the results obtained for the thin films of $\alpha$- and $\beta$-MADN, and a discussion of the effects of long-range interactions on the quasiparticle energies, the layer-resolved and surface density-of-states, and the final UPS simulations including carrier-vibration coupling. A brief summary (Section IV) concludes the paper.

II Methodology

II.1 Atomistic thin-film morphologies

As a first step, realistic thin-film morphologies are obtained using the Metropolis Monte Carlo-based simulated annealing protocol DEPOSIT Neumann et al. (2013). It mimics the vapor deposition (PVD) process and provides molecular morphologies that exhibit commonly observed PVD characteristic features Friederich et al. (2017, 2018). Each molecule was deposited using 32 simulated annealing cycles with 120000 Monte-Carlo steps each, with annealing temperatures decreasing from 4000 K to 300 K. Periodic boundary conditions were applied in the directions perpendicular to the growth direction with a side length equal to 10 nm. The deposition substrate is represented by a fixed dense layer of MADN. The energy at each simulation step was computed using customized force-fields generated using the Parametrizer module of the DEPOSIT code. These force-fields comprise Coulomb electrostatics based on partial charges obtained from an electrostatic potential fit Besler et al. (1990), Lennard-Jones potentials to account for the van-der-Waals interaction and the Pauli repulsion, and compound-specific dihedral force-fields with quantum chemistry accuracy generated by the Dihedral Parametrizer of the DEPOSIT code. The final deposited morphologies contain 1000 molecules and are about 10 nm thick. For the following analysis, we remove the bottom 2 nm of the film to avoid spurious effects from the artificial substrate.

II.2 Quasiparticle Energies in the GW Approximation

Next, the internal contributions to the HOMO energy of all individual molecules are calculated including quasiparticle corrections within the $GW$ approximation of many-body Green’s Functions theory Sham and Rice (1966); Hedin and Lundqvist (1970). At this level, the calculations already include the effects of molecular deformations, obtained from the morphology simulations, but not yet the effects of embedding in the polarizable thin-film environment. These will be discussed in Section II.3. Throughout this subsection, atomic (Hartree) units are used ($\hbar=1$, $m_{e}=1$ and $e^{2}/(4\pi\epsilon_{0})=1$, with $m_{e}$ the electron mass, $e$ the elementary electron charge, and $\epsilon_{0}$ the vacuum permittivity).

The properties of a closed shell system of $N$ electrons with spin singlet ground state can be calculated using DFT by solving the Kohn–Sham (KS) equations Onida et al. (2002)

[TABLE]

with $T_{0}$ the kinetic energy operator, $V_{\text{ext}}$ an external potential, $V_{\text{Hartree}}$ the Hartree potential, $V_{\text{xc}}$ the exchange-correlation potential, and $\epsilon^{\text{KS}}$ (${|\phi^{\text{KS}}\rangle}$) the KS energies (wave functions), respectively. Particle-like excitations, known as quasiparticles (QP), in which one electron is added to or removed from the $N$-electron ground state, are described by the one-body Green’s function, $G_{1}$ Sham and Rice (1966); Hedin and Lundqvist (1970) that obeys a Dyson-type equation of motion. This, in spectral representation, is given by

[TABLE]

with $H_{0}=T_{0}+V_{\text{ext}}+V_{\text{Hartree}}$ and with $\Sigma(E)$ an electron self-energy operator that describes the exchange-correlation effects. It can be shown that Eq. (3) is part of a closed set of coupled equations, known as Hedin equations Hedin (1965); Strinati (1988). KS-DFT is an approximate solution for the excited-electrons problem, in which $\Sigma\sim V_{\text{xc}}$. A beyond-KS-DFT solution of Eq. (3) is given by the GW approximation, in which the self energy is expressed as

[TABLE]

Here, $W({\mathbf{r}},{\mathbf{r}}^{\prime},\omega^{\prime})=\int\epsilon^{-1}({\mathbf{r}},{\mathbf{r}}^{\prime\prime},\omega^{\prime})\,v_{c}({\mathbf{r}}^{\prime\prime},{\mathbf{r}}^{\prime})\,d{\mathbf{r}}^{\prime\prime}$ is the screened Coulomb interaction between unit charges at position $\mathbf{r}$ and $\mathbf{r}$’, with $v_{c}({\mathbf{r}},{\mathbf{r}}^{\prime})=1/|{\mathbf{r}}-{\mathbf{r}}^{\prime}|$ being the bare Coulomb interaction and $\epsilon^{-1}({\mathbf{r}},{\mathbf{r}}^{\prime\prime},\omega)$ the inverse of the frequency-dependent dielectric function. We calculate the latter in the random-phase approximation (RPA) Hybertsen and Louie (1985). Employing the $GW$ approximation, we express the quasiparticle wave functions as linear combinations of KS states and obtain their energies $\varepsilon_{i}^{\text{QP}}$ after diagonalization of the energy-dependent $GW$ Hamiltonian Aulbur et al. (2000); Rohlfing (2000)

[TABLE]

In a perturbative treatment of Eq. (5), in which it is assumed that ${|\phi^{\text{QP}}_{i}\rangle}\approx{|\phi^{\text{KS}}_{i}\rangle}$, the quasiparticle energies are determined as diagonal elements of $H_{ij}^{GW}$. This leads to

[TABLE]

A comparison between these two approaches is given in Sect. S2 of the Supplemental Material sup for the two MADN derivatives.

As the self-energy is energy-dependent, and thus depends on $\varepsilon_{i}^{\text{QP}}$, the solution to Eq. (5) or Eq. (6) need to be found self-consistently. Both the correction term $\Delta\epsilon_{i}^{GW}$ and the non-local, energy-dependent microscopic dielectric function calculated within the RPA depend on $\varepsilon_{i}^{\text{QP}}$ (Hybertsen and Louie, 1985; Rohlfing et al., 1995). Within the ${GW}$ method employed in this work, we iteratively solve the Hamiltonian by updating as well the energy-dependent non-local dielectric function, until self-consistency of the eigenvalues (ev) is obtained. In literature this is often referred to as $\text{ev}GW$ Wehner et al. (2018).

We have performed the KS-DFT and the $GW$ steps using the ORCA Neese (2012) and VOTCA-XTP Wehner et al. (2018) software packages, respectively. VOTCA-XTP can read information from standard packages, as such as ORCA, using Gaussian-type orbitals as basis functions $\left\{\psi_{i}({\mathbf{r}})\right\}$ to express

[TABLE]

Matrix elements ${\langle\phi_{i}^{\text{KS}}|}V_{\text{xc}}{|\phi_{j}^{\text{KS}}\rangle}$ needed in Eq. (6) are numerically integrated using spherical Lebedev and radial Euler-Maclaurin grids as used in NWChem Valiev et al. (2010), with XC functionals provided by the LibXC library Marques et al. (2012).

In the evaluation of the self energy, four-center Coulomb integrals of the form

[TABLE]

need to be calculated. VOTCA-XTP makes use of the resolution-of-identity approximation (namely the RI-V approximation) to reduce the scaling from $N_{\text{b}}^{4}$ to $N_{\text{b}}^{3}$, with $N_{\text{b}}$ is the number of basis function. An auxiliary basis set $\{\chi_{\nu}({\mathbf{r}})\}$ is introduced so that the integrals in Eq. (8) are rewritten in the form

[TABLE]

Here, $(\nu|\mu)^{-1}$ are elements of the inverse of the two-center Coulomb matrix $(\nu|\mu)=\iint d{\mathbf{r}}_{1}\,d{\mathbf{r}}_{2}\,\chi_{\nu}({\mathbf{r}}_{1})v_{c}({\mathbf{r}}_{1},{\mathbf{r}}_{2})\chi_{\mu}({\mathbf{r}}_{2})$ and $(ij|\nu)=\iint d{\mathbf{r}}_{1}\,d{\mathbf{r}}_{2}\psi_{i}({\mathbf{r}}_{1})\psi_{j}({\mathbf{r}}_{2})v_{c}({\mathbf{r}}_{1},{\mathbf{r}}_{2})\chi_{\mu}({\mathbf{r}}_{2})$ are elements the three-center Coulomb matrix, respectively.

The KS-DFT eigenvalues, and thus the quasiparticle energies, may depend strongly on the exchange-correlation functional used. However, for both MADN isomers, the final $GW$ results show a negligible starting-point dependence. A comparison using the PBE functional Perdew et al. (1996) and the hybrid PBEh Adamo and Barone (1999); Ernzerhof and Scuseria (1999) is given in Sect. S2 of the Supplemental Material sup . All results reported in this paper, are obtained using the PBE functional and the cc-pVTZ basis Kendall et al. (1992) with its optimized auxiliary basis set Weigend et al. (2002) for resolution-of-identity techniques.

The numerical accuracy of the calculations depends on the convergence limit used, the number of levels included and the method for carrying out the frequency integration in Eq. (4). The convergence limit for the self-consistent $GW$-cycles in the $\text{ev}GW$ scheme was set to $10^{-5}$ Hartree (0.27 meV). The number of occupied and unoccupied levels taken into account for the QP calculations is 327, while for the calculation of the RPA dielectric function the full spectrum of the KS states (1385 levels) is used. The frequency integration in Eq. (4) can be performed in VOTCA-XTP using the Fully Analytical Approach Golze et al. (2019) (FAA) or a generalized plasmon-pole model (PPM) Rohlfing et al. (1995). The FAA expresses the frequency dependence of the self-energy in the eigenbasis of the full RPA Hamiltonian, which is in turn evaluated in the basis of KS product states. This approach is in principle exact. However, disadvantageously, the $N_{\text{b}}^{6}$ scaling of the FAA Golze et al. (2019) makes its application to molecules of the size of MADN computationally extremely demanding. As an alternative, the PPM allows for a fast evaluation of the self-energy. For inorganic semiconductors, the quasiparticle energy obtained using the PPM can show deviations of several tenths of eV from the exact result Golze et al. (2019); Larson et al. (2013). However, for geometry-optimized $\alpha$-MADN and $\beta$-MADN, we find a difference of only 0.02 eV between the HOMO energies obtained using the FAA and PPM approaches. All things considered, the results reported in this paper have been obtained using the PPM.

II.3 QM/MM quasiparticle embedding schemes

Intermolecular interactions in the thin-film environment give rise to additional non-uniform modifications of the quasiparticle energies. We determine the respective corrections to the intramolecular $GW$ energies in a coupled quantum-classical (QM/MM) procedure (“quasiparticle embedding”) May et al. (2012); Baumeier et al. (2014); Schwabe et al. (2012); Li et al. (2016, 2018); Wehner et al. (2018). Within the MM model, we employ a classical representation of the molecular electrostatic potential based on static and induced multipole moments, located on each of the atoms in the system. A region treated on QM (here $GW$) level is coupled to suitably defined MM regions with a special scheme, which properly includes the long-range character of electrostatic interactions among the excited and neutral MADN molecules. Below, the details of the approach are discussed.

The classical environmental contribution to the potential (the MM part) follows from the static atomic multipole moments $Q^{a}_{t}$ Stone (2013), where $t$ indicates the multipole rank and $a$ the associated atom in the molecule $A$, and from the induced moments $\Delta Q_{t}^{a}$ due to the field generated by moment $t^{\prime}$ of atom $a^{\prime}\neq a$ in molecule $A$ and the one generated by the moment $u$ of atom $b$ in molecule $B$:

[TABLE]

with $\alpha_{tt^{\prime}}^{aa^{\prime}}$ the atomic polarizability on each site and $T_{t^{\prime}u}^{a^{\prime}b}$ the tensor describing the interactions between the multipoles moments $Q^{a^{\prime}}_{t^{\prime}}$ and $Q^{b}_{u}$. A repeated-index summation convention is used for the multipole indices $t,t^{\prime},u$. The classical total energy due to the interaction between molecules with atomic indices $a$ and $b$ in regions $\mathcal{A}$ and $\mathcal{B}$ (case $\mathcal{A}=\mathcal{B}$ included) is given for the state $(s)$ (i.e neutral ($s=\text{n}$) or charged ($s=\text{qp}$, via quasiparticle excitation) by Stone (2013)

[TABLE]

Static atomic partial charges from a CHELPG fit to the neutral molecule’s electrostatic potential Breneman and Wiberg (1990) are used for the classical representation of molecules, and atomic polarizabilities, optimized to reproduce the polarizable volume of the molecule obtained from DFT, account for polarization effects via the induction of atomic dipoles (Thole model Stone (2013)).

Various schemes for coupling a QM-treated inner region to a MM-treated outer region have been described in the literature. Within an additive scheme Cao and Ryde (2018) (here termed $GW/\text{aMM}$), the potential of the MM environment is explicitly included in the $GW$ calculation as an additional external potential to the Hamiltonian. The QM region is directly polarized by the multipole distribution (and vice versa) and coupled solutions are found self-consistently Wehner et al. (2018). Within a subtractive scheme ($GW/\text{sMM}$), the QM region is replaced by a MM representation of the different states and a purely classical energy correction $E_{\text{MM}}^{(\text{n})}-E_{\text{MM}}^{(\text{qp})}$ is added to the $GW$ vacuum energies. When applying this scheme, we approximate the state of the MADN molecules after the creation of an hole by that of the cation. In both cases, the MM environment includes all molecules inside a region within a cutoff distance $r_{\text{c}}$ around the QM molecule.

However, both cutoff-based techniques rely on the assumption that only short-ranged local interactions affect the energies of the QM region. In the thin-films of MADN, this is not the case as the deposition simulations reveal a weak net ordering of the small molecular dipole moment (0.59 D and 0.56 D for the $\alpha$ and $\beta$ isomers). For the simulated morphologies considered, the cumulative calculated electric dipole moment parallel to the surface normal ($z$-direction) of in total 65.6 D ($\alpha$-MADN) and 35.7 D ($\beta$-MADN). A detailed analysis of the respective distributions resolved on molecular scale is available in Sect. S3 of the Supplemental Material sup . The two aforementioned cutoff-based approaches cannot account for long-range electrostatic effects Kratz et al. (2016) that result as a combination of the thin-film geometry and cumulative electrostatics. Our final calculations are therefore based on a third scheme, here termed $GW\!/\text{pMM}$. This is an extension of the $GW/\text{sMM}$, in which the long-range electrostatic interaction effects are included via an infinite periodic embedding based on the traditional classical Ewald summation method Arnold and Holm (2005).

We show in Section III that this scheme allows us to include electrostatic interactions up to an arbitrarily large cutoff distance. When only the static point charges are considered in both regions, we call this the ”static $GW\!/\text{pMM}$ scheme”. The final results are obtained by also including polarizable (polar) interactions up to a cut-off distance $r_{c}^{p}$ of 3 nm (”polarizable $GW\!/\text{pMM}$ scheme”). Outside that radius, $\alpha_{tt}^{aa^{\prime}}=0$. Including these polarizable interactions reveals that the molecular dipole moments are slightly screened, so that the total accumulated dipole moment in the film is reduced by 3.5 D ($\alpha$-MADN) and 2.8 D ($\beta$-MADN).

II.4 Surface Density of States and Carrier-Vibration Coupling

As a final step, we calculate from quasiparticle energies the HOMO contribution to the UPS spectrum, using Eq. (1), by taking the surface-sensitivity of the experiment and the carrier-vibration coupling into account. In the absence of the carrier-vibration coupling, the spectrum would be proportional to the surface density-of-states (SDOS). This is obtained by weighting the $z$-dependent frontier orbital energy with the exponential function from Eq. (1), i.e.

[TABLE]

We will adopt a value of $\Lambda=1\,\mathrm{nm}$ as suggested from experiment. It is compatible with estimates of the inelastic mean free path of the electrons within the random-phase approximation using $GW$ energies, which we consider an upper limit to the electron attenuation length (see Sect. S4 of the Supplemental Material sup and Ref. chr ). The SDOS does not include the effect of the intramolecular reorganization process upon charge removal and the associated shift and lineshape broadening via carrier-vibration coupling. Conventionally, the effect of the intramolecular reorganization process upon charge removal and the associated shift and lineshape broadening via carrier-vibration coupling is described using semi-classical Marcus theory Gerischer and Willig (1976); Pope and Swenberg (1999). The spectral shape due to the coupling of the photoelectrons with vibrational modes is then given by

[TABLE]

with $\lambda_{\text{tot}}$ total reorganization energy, $k_{\text{B}}$ the Boltzmann constant and $T$ the temperature. The energy difference is defined relative to the adiabatic excitation energy, $\varepsilon_{\text{HOMO}}^{\text{ad}}=\varepsilon_{\text{HOMO}}^{\text{vert}}+\lambda_{\text{tot}}$ (see the schematic representation in Fig. 2). However, significant coupling with vibrational modes with energies well above $k_{\mathrm{B}}T$, such as the C-C stretch vibrations of the phenyl rings in the 0.1 – 0.2 eV range, makes the semiclassical approach for most organic semiconductor materials invalid. This has been demonstrated for the related problems of the rates of electron or hole hopping and exciton transfer de Vries et al. (2018); de Vries et al. (2019). Analogous to the full-quantum (FQ) approach for inter-molecular charge transfer de Vries et al. (2018), which approximates the potential energy surface of the excited molecule in the independent mode displaced harmonic oscillator model Petrenko and Neese (2007), we find that the spectral shape due to coupling of the photoelectrons with vibrational modes $k$ with an energy $\hbar\omega_{k}$ and a coupling energy $\lambda_{k}$ is given by

[TABLE]

where $F(t)=\sum_{k}\frac{\lambda_{k}}{\hbar\omega_{k}}\left(\coth\left(\frac{\hbar\omega_{k}}{2k_{\text{B}}T}\right)\cos(\omega_{k}t)+i\sin(\omega_{k}t)\right)$. Evaluating the needed parameters (vibrational modes $\omega_{k}$ and coupling energy $\lambda_{k}$) can be simplified under the assumptions that (i) the ground- and excited-state potential energy surfaces are harmonic and (ii) no vibrational frequency alteration or normal mode rotation occurs in the excited state. With these two assumptions, equivalent to the premise of a linear electron-phonon coupling, only the ground-state vibrational modes frequencies $\omega_{k}$ and the gradient of the total energy for the charged (excited) system in the ground state geometry with respect to the phonon mode coordinates $q_{k}$, $\frac{\partial E}{\partial q_{k}}$, need to be evaluated as indicated in Fig. 2. The mode-specific coupling energies are then determined as $\lambda_{k}=\frac{1}{2}\frac{\partial E}{\partial q_{k}}$.

III Results and discussion

III.1 Quasiparticle embedding and long-range interactions

To assess the influence of long-range electrostatic interactions on the quasiparticle energies of the two MADN thin films, we analyze the results obtained with the different embedding schemes introduced in Sec. II.3. First, we consider the effects of the dependence of the cutoff radius in the additive $GW/\text{aMM}$ scheme. All results given in this subsection are obtained using only a static embedding approach, within which the polarizability of the molecules is switched off. We show in Fig. 3 (blue curve) the calculated $\epsilon_{\text{HOMO}}$ as a function of $r_{\text{c}}$ for a molecule at the surface and in the bulk-like region of the $\beta$-MADN film, respectively. The slow decrease of the energy with size of the embedding region indicates, in particular for the surface molecule, that even at a cutoff of 20 nm no converged result is obtained. Repeating the same analysis for the computationally less demanding subtractive $GW$/sMM scheme, in which we now allow the classical substitute of the QM molecule to be polarizable as it automatically is in $GW$/aMM, yields a cutoff dependence given by the red lines in Fig. 3. Comparison to the $GW/\text{aMM}$ data reveals deviations smaller than 0.02 eV. Based on this good agreement, we then use this parametrization and embed the classically represented QM molecule in a periodically repeated background in the $GW$/pMM setup. The results for both the surface and the bulk molecule are shown as data points for $r_{\text{c}}=\infty$ in  Fig. 3. For the surface molecule, the periodically embedded HOMO energy is 0.07 eV lower than as obtained from the $GW$/sMM calculation with the largest cutoff considered. In the bulk, the difference is slightly smaller, viz. 0.03 eV. We adopt the $GW\!/\text{pMM}$ scheme for the following analysis for thin-film energy levels, the SDOS and the UPS spectrum.

III.2 Layer-resolved energy levels and DOS

Figure 4 shows the laterally-averaged depth dependence of the HOMO energies as resulting from various levels of refinement, labeled ”0” to ”3”. It shows that ${GW}$ corrections (”level 1”) to the vacuum KS levels (”level 0”) lower the energies by up to 1.9  eV, nearly uniformly for both isomers. The gas-phase simulations include the molecular deformations in the thin-film morphology, but these cause only a small broadening of the DOS. It is mainly due to disorder in the anthracene-naphthalene torsion angle, which is largest for $\beta$-MADN. When long-range electrostatic interactions are included (static $GW\!/\text{pMM}$, level “2”, blue squares in Fig. 4), we find for both isomers a nearly linear $z$-dependence of the mean HOMO energy, which is symmetric with respect to the mean vacuum ${GW}$ energy $\varepsilon^{GW}$. This is due to accumulating net dipole moment contributions parallel to the surface normal during film growth. Even though the dipole moments of the individual molecules are small (0.59 D and 0.56 D for the $\alpha$ and $\beta$ isomers) and their average net components parallel to the growth axis are only 0.065 D and 0.040 D, respectively, the resulting energy gradients are a few hundredths of an eV/nm. Adding polarization effects (polar $GW\!/\text{pMM}$, level “3”, red circles in Fig. 4) leads to a shift of the mean of the distributions to lower binding energies. The effect is stronger in the bulk-like center of the film (0.7 eV (0.6 eV) for $\alpha$-MADN ($\beta$-MADN)) than at the vacuum surface (0.5 eV). The $z$-coordinate dependence of the distribution of the individual HOMO energies, obtained using the static and polarizable $GW\!/\text{pMM}$ schemes is shown in Sect. S3 of the Supplemental Material sup .

Differences between the surface and bulk energy level structure of organic materials are well-known from UPS studies. For crystalline anthracene, e.g., the experimental binding energy difference for the first and second monolayer was found to be $0.3\pm 0.15\,\mathrm{eV}$ Salaneck (1978). A similar effect is seen in Fig. 4 for the first molecular layer near the vacuum surface, which after subtracting the energy gradient due to dipole orientation shows an increase of the binding energy of about 0.15 eV. Since we simulate freestanding thin films, we also note modifications from the bulk-like behavior at the bottom surface with $z=0$. In view of our interest in analyzing the UPS spectrum after irradiation from the positive $z$ direction, which is only sensitive to the energy level structure in a thin zone near the top vacuum surface, we focus on that region.

III.3 Vertical and Adiabatic Surface Density of States

Fig. 5(a,b) show the SDOS for the two isomers (light-red shaded), obtained using Eq. (12). The SDOS is based on the vertical excitation energies. These include the effect of the electronic polarization of the environment but do not include the effect of the intramolecular structural reorganization process upon charge removal. The ad-SDOS curves in Figure 5 (a,b) (dark-red shaded) show the adiabatic SDOS, obtained from the SDOS by adding an energy shift equal to the total reorganization energy, which is 0.21 eV for both isomers. The actual excitation process is not adiabatic, but is accompanied by the excitation of vibrational modes. Figure 5(c) shows, for a single $\beta$-MADN molecule, the resulting shift of the DOS to a more negative HOMO energy and the resulting broadening to a full-width at half maximum of approximately 0.25 eV. For $\alpha$-MADN, the effect is very similar. The absence of multiple satellite peaks, as usually seen for gas-phase spectra of structurally more simple molecules such as pentacene, is due to coupling of multiple modes, which smears out the structure. The peak of the function $S_{\textrm{el-vib}}^{\textrm{FQ}}$ is almost equal to the vertical HOMO energy. Including vibrational effects thus induces for these systems almost no peak shift with respect to the vertical (polarizable $GW$/pMM) energies, but only a spectral broadening. The resulting calculated UPS spectral intensity function, given by Eq. (1), is shown in Figure 5(a) and Figure 5(b) by the curves labeled “ad-SDOS+FQ”.

It is interesting to compare the full-quantum vibrational response function $S_{\textrm{el-vib}}^{\textrm{FQ}}$ to the classical Marcus function. According to Eq. (13), the peak energy is shifted to lower energies relative to the adiabatic excitation energy by the total reorganization energy of 0.21 eV and the function has a Gaussian shape with a full-width at half maximum of $4\sqrt{\ln(2)\lambda_{\textrm{tot}}k_{\textrm{B}}T}\cong 0.24\,\mathrm{eV}$. A comparison with Fig. 5(c) shows that in this case the full-quantum result differs only weakly from the classical Marcus response function. Apparently, the coupling of the photo-electrons to high-energy vibrational modes is for these systems relatively weak.

III.4 Simulated UPS including Carrier-Vibration Coupling

Figure 6 shows the final UPS spectra of the frontier orbital of $\alpha$- and $\beta$-MADN thin films, simulated for the four different levels of theory, together with the experimental data. Characteristics of the signals (maximum position, onset, and FWHM) are listed in Tab. 1. Following the conventional approach, the onset energy is defined by extrapolating the tangent through the low-binding-energy inflection point of the HOMO peak to zero intensity. The energetic position of the simulated peaks for the different methods reflects the variations discussed for the layer-averages in Fig. 4. Comparison to the reference experimental spectrum now allows assessment of the quality of the various methods and the importance of the individual processes for the analysis of the experiment.

Simulations based on vacuum energies which exclude the effects of inhomogeneous local electric fields and environment polarization either over- (KS) or underestimate ($GW$) the energy of the peak maximum by up to 1.3 eV. The FWHM is nearly exclusively determined by the single-molecule spectral function and results about a third smaller than measured. Inclusion of static local field effects in $GW\!/\text{pMM}$ does not noticeably affect the peak maximum but the additional disorder contributes to further broaden the signal. Accounting for the polarization response of the material upon quasiparticle excitation in $GW\!/\text{pMM}$ we obtain a simulated UPS signal in excellent agreement with the measurement: the largest deviation is found for the peak maximum of $\beta$-MADN and amounts to only 50 meV, which is also the experimental resolution. Most importantly, the comparison emphasizes that it is possible to achieve a predictive modeling of frontier orbital energies at the accuracy needed for an accurate understanding and prediction of device performance. The fact that we achieve the same accuracy by studying two isomers, with different molecular structures and thin film morphologies, supports the robustness of our approach.

In device simulations, the bulk adiabatic ionization energy, $\varepsilon_{\textrm{HOMO,bulk}}^{\text{ad}}$, is needed. This energy may be obtained from a linear extrapolation of the bulk polarizable $GW\!/\text{pMM}$ energies shown in Fig. 4 to the surface plane at $z=z_{0}$, plus the reorganization energy $\lambda_{\text{tot}}$ of 0.21 eV. The resulting values, $\varepsilon_{\textrm{HOMO,bulk}}^{\text{ad}}=-5.91(-5.89)\pm 0.05\,\mathrm{eV}$ for $\alpha(\beta)$-MADN, are located in between the UPS peak and onset energies (see Tab. 1). Using either the peak or the onset value would, in this case, thus introduce an error of about 0.1 eV or more. Our refined protocol for the analysis of UPS measurements provides a methodology for avoiding such an error. One may see from Fig. 5(c) that when carrying out a measurement for a single molecule, the onset energy would provide an excellent approximation to $\varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}}$. However, that coincidence is fortuitous, as the peak shape and onset energies depend on the mode-resolved reorganization energies and the temperature. That may already be seen when considering the peak shape obtained within the semiclassical Marcus-theory (Eq. (13)). The difference between the onset energy and the adiabatic ionization energy is then equal to $(-\lambda_{\textrm{tot}}+\sqrt{8\lambda_{\textrm{tot}}k_{\textrm{B}}T})$. For a system with $\lambda_{\textrm{tot}}=0.2\,\mathrm{eV}$ (close to the value for MADN) and for $k_{\textrm{B}}T=0.025\,\mathrm{eV}$ (close to room temperature), the onset and adiabatic ionization energies then indeed coincide. However, the electron-vibrational mode coupling shows a significant dependence on the molecules considered. The total reorganization energies vary from less than 0.1 eV to more than 0.3 eV, with a tendency to decrease with increasing molecular size Devos and Lannoo (1998); Deng and Goddard (2004); Kera and Ueno (2015). Furthermore, in thin films the spectral broadening due to energetic disorder leads to a shift of the onset value to a smaller binding energy, whereas the reduced screening at the thin film surface leads to shift to a larger binding energy. For the two isomers of MADN, the former effect is largest, so that the absolute value of the onset energy is slightly smaller than $|\varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}}|$.

IV Summary and Conclusions

In summary, we have developed a multiscale approach that provides a prediction of the frontier orbital UPS spectrum of amorphous molecular thin films. The approach includes (i) first-principles calculations of a realistic thin-film morphology, (ii) the electronic properties at a state-of-the-art (many-body Green’s functions) level of quantum chemistry, (iii) embedding in a polarizable molecular mechanics environment, and (iv) the effects of the vibrational modes that are excited in the experiment. We have focused on two isomers of MADN, for which the non-degeneracy of the HOMO state leads to an exceptionally narrow width of the HOMO orbital UPS spectrum. Our work shows how the spectrum is related to the disorder-induced energy level distribution in the bulk of the organic semiconductor and near the surface. The good agreement between the calculated and experimental peak positions and widths of the UPS spectra indicates that our approach provides a route towards accurately predicting the bulk adiabatic ionization energy $\varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}}$, which is the HOMO energy needed in device simulations.

We find that neither the HOMO energy that would follow from the UPS peak energy nor the onset energy coincide with $\varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}}$. For the MADN films studied, both assumptions would introduce an error of about 0.1 eV or more. Instead, the actual value of $\varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}}$ is in this case intermediate between the onset and peak values. The error made by taking $\varepsilon_{\textrm{HOMO,bulk}}^{\textrm{ad}}$ equal to the UPS onset energy is determined by the balance between the effects of energetic disorder and reduced screening at the thin film surface. For device applications, the relative error between different materials is most important. From our study, we expect that the relative error when taking the onset energy is largest for the case of two materials with strongly dissimilar energetic disorder energies or strongly dissimilar molecular polarizabilities.

Acknowledgments

This work has been supported by the Innovational Research Incentives Scheme Vidi of the Netherlands Organisation for Scientific Research (NWO) with project number 723.016.002 and by Horizon-2020 EU project MOSTOPHOS (Project No. 646259, XdV)