Modeling X-ray Photoelectron Spectroscopy of Macromolecules Using GW

We propose a simple additive approach to simulate X-ray photoelectron spectra (XPS) of macromolecules based on the $GW$ method. Single-shot $GW$ ($G_0W_0$) is a promising technique to compute accurate core-electron binding energies (BEs). However, its application to large molecules is still unfeasible. To circumvent the computational cost of $G_0W_0$, we break the macromolecule into tractable building blocks, such as isolated monomers, and sum up the theoretical spectra of each component, weighted by their molar ratio. In this work, we provide a first proof of concept by applying the method to four test polymers and one copolymer, and show that it leads to an excellent agreement with experiments. The method could be used to retrieve the composition of unknown materials and study chemical reactions, by comparing the simulated spectra with experimental ones.


TOC Graphic
C 1s X-ray − X-ray photoelectron spectroscopy (XPS) is a well-established experimental technique to study the chemical properties of solid, liquid, and gaseous materials [1,2]. In XPS, core-electron binding energies (BEs) are measured by means of the photoelectric effect. A typical XPS spectrum shows peaks at various BEs, with the peak height related to the number of electrons emitted from core orbitals at those BEs. The BEs are sensitive to the local chemical environment, such as the type of adjacent atoms, bonds, oxidation states. Therefore, absolute and relative peak positions can be used to retrieve information on the chemical composition of the material. The standard procedure to interpret XPS spectra consists in fitting the peaks to an envelope of Voigt functions with varying amounts of Gaussian and Lorentzian character [3]. However, the fitting procedure becomes increasingly challenging when more complicated compounds are considered, as the peaks tend to overlap and become indistinguishable due to the limited experimental resolution. The development of theoretical techniques to model XPS spectra is there-fore important to support the interpretation of experiments [4,5,6].
A common method to simulate XPS spectra is the Delta self-consistent field method (∆SCF) [5,4]. In this technique, the BEs are calculated as the energy difference between the neutral species and the ionized species, where the energy of the latter is obtained after fixing a hole in a core orbital and relaxing the outer orbitals. Most of the ∆SCF schemes make use of Hartree-Fock (∆HF) or density functional theory (∆DFT) [5,4]. Similar techniques based on coupled-cluster calculations (∆CC) are also used on small molecules [7]. The disadvantage of these approaches is the difficulty of artificially creating a hole in each core orbital, which can easily lead to convergence issues and become cumbersome for larger molecules.
Another promising technique for calculating core-electron binding energy is GW , which is based on many body perturbation theory [8,9,10,11]. GW allows to calculate BEs directly from the one-particle Green's function, without explicitly generating holes in the core orbitals. The one-particle Green's function G(xt, x t ) is defined as follows: if t > t , G is the probability amplitude to find an electron at x at time t after addition of an electron at x at time t , where x and x include both spatial and spin coordinates; if t < t , G is the probability amplitude to find a hole at x at time t after removal of an electron at x at time t. By definition, the poles of the Green's function correspond to the electron removal (attachment) energies as measured by (inverse) photoelectron spectroscopy. In principle, calculating the poles of the Green's function requires a fully self-consistent iterative solution of the Hedin equations [12]. To reduce the computational cost, the iterative procedure can be limited to a subset of variables, yielding different flavours of GW , such as single-shot GW (G 0 W 0 ), evGW , scGW , evGW 0 , and scGW 0 [10]. In the following, we consider only the least expensive G 0 W 0 , which corresponds to the first iteration of Hedin equations. As a starting point for this single-shot perturbation calculation, we consider Kohn-Sham states and eigenvalues calculated using DFT.
For valence electrons, the BEs calculated from G 0 W 0 are in good agreement with experiments [13,14]. Only recently, the accuracy of the G 0 W 0 method has been confirmed also for coreelectron BEs [15,16,17]. Although G 0 W 0 can potentially overcome the limitations of ∆SCF and ∆CC, its application to macromolecules is still unfeasible due to the high computational cost associated with these calculations. Therefore, approximations should be devised to use this method for larger molecular systems. Here, we propose a simple additive approach to extend the range of applicability of G 0 W 0 to macromolecules. As a proof of concept, we test the methodology on the four non-conjugated polymers in Scheme 1. These polymers have important applications in microelectronics as photoresists for extreme ultra-violet lithography [18], besides being good examples of C, H, O materials containing carboxyl, methyl, hydroxyl, and benzene groups. In first approximation, we can simulate a XPS spectrum as the sum of Gaussian peaks of width σ centered at the BEs obtained from G 0 W 0 calculations performed on top of DFT (G 0 W 0 @DFT). In reality, additional effects such as electron scattering in the material or electronphonon interactions can add features to the experimental spectra, such as a broad background and satellite peaks, which are neglected here. We assume that the dipole matrix elements are constant, or, equivalently, that all core electrons within the BE region of interest have the same photoionization cross-section [9,10]. This is justified by the fact that all core orbitals under study are expected to show an approximately equivalent coupling with the X-ray photons, due to their 1s-like spherical symmetry. Each orbital thus contributes to the photocurrent in proportion to the number of electrons occupying that orbital. Therefore, the peaks are weighted by the spin degeneracy of the corresponding orbital, which corresponds to a factor 2 for the molecules under study, as all orbitals are doubly-occupied. Since the wavefunctions of core electrons are strongly localized near the atomic nuclei, we can neglect the interaction between core electrons on distant atoms. This allows us to approximate  the macromolecule as a collection of N mol independent building blocks of molar ratio w i . The sum over all core-electron BEs can then be split over separate parts of the macromolecule, and the theoretical photocurrent simplifies to: where C is a normalization factor to match the peak height with the experimental data, the first sum runs over the N mol building blocks of the macromolecule, and the second sum runs over the N i core orbitals of the i-th building block characterized by binding energy n i and spin degeneracy g n i . In this approximation, the calculation of the BEs can be performed separately on each building block, saving computational time.
Of course, the building blocks should be chosen of the appropriate size in order to include the local chemistry while remaining tractable. In the case of polymers, this can be usually done by considering isolated monomers. A schematic of the calculation process is depicted in Scheme 2. Eq (1) is based on the assumption that the interactions between the building blocks are negligible. This assumption is expected to be generally valid for systems consisting of molecules that are not bonded together, and for non-conjugated polymeric materials where the core-electron BEs are not significantly affected by the monomers being connected to the polymer backbone. For conjugated polymers, however, this approximation could fail due to the strong electron delocalization. On top of this, the spectra of conjugated polymers are also complicated by additional effects such as oxidation and polaron states formation [19], which are currently not included in our approach.
For the four polymers under study, we investigate the validity of the independent-blocks assumption by computing the spectra of DFT optimized geometries of increasingly long poly- To assess the impact of adjacent backbone atoms, we also consider methyl-terminated monomers (CH 3 -M -CH 3 ) and dimers (CH 3 -M 2 -CH 3 ), where the adjacent backbone atoms in the polymer chain are replaced with methyl groups. In the following, we will refer to these five types of cutting as "backbone corrections". Of course, the addition of two methyl groups in CH 3 -M -CH 3 and CH 3 -M 2 -CH 3 introduces two extra core levels in the calculated C 1s spectra. To avoid this artifact, we used a visualization software to identify the two C 1s orbitals localized on the two methyl groups, and removed the contributions of the corresponding BEs before computing the spectra. The core-electron Kohn-Sham orbitals and calculated BEs for the isolated monomers are depicted in Figures S1-S4.  tatively, different approximations lead to differences in peak positions of up to 0.4 eV (Figure S5). Interestingly, the calculated peaks for styrene-based polymers, namely PS and PHS, tend to systematically move to smaller BEs when longer chains are considered. On the contrary, the effect of different backbone corrections on acrylate-based polymers, i.e. PMMA and PBMA, does not seem to follow a clear trend.
Another factor that might affect the calculated BEs is the molecular geometry. The presence of conformers in a sample could, in principle, be one of the causes of experimental broadening. To assess the impact of the molecular geometry on the core-electron BEs, we compare the spectra computed on nine different conformers of isolated monomers of the polymers under study. The conformers were extracted from a model polymer matrix of 50 repeating units optimized with DFT (see Computational details). Hydrogen atoms were added to the extracted monomers along the backbone direction and their position was optimized with DFT while keeping the other atoms fixed. G 0 W 0 was then used to calculate the BEs. Figures 1(e-h) show that the peak variation over different conformers is smaller for styrene-based than for acrylate-based polymers, probably due to the rigidity of the aromatic group, which leads to a smaller conformational space. Overall, the variation due to the presence of conformers is generally below 0.04 eV, except for PBMA, where it is up to 0.16 eV (Tables S1-S4 and Figures S6-S9). These values are smaller than a typical experimental broadening (0.5 eV). Therefore, the peak broadening in XPS spectra most likely arises from other sources, such as finite excitation lifetime and limited experimental resolution.
The results shown above suggest that the XPS spectra of the polymers under study can be conveniently simulated by performing the calculations on much shorter chains. To assess the validity of this approximation, we compare our results with reported experimental spectra of two polymers: one acrylate-based (PMMA [20]) and one styrene-based (PHS [21]) polymer. Both polymers were chosen for the existence of multiple features in the spectra and for the relatively small size of their monomers, which allow for a more complete assessment of various backbone corrections. For PMMA, both C 1s and O 1s edges were considered, whereas only the C 1s edge was investigated for PHS, as the O 1s edge shows only one peak and is thus less informative.
A fitting procedure was performed to compare the calculated spectra with experiments. A rigid shift ∆ was introduced to match the theoretical BEs with the experiment. The origin of the shift will be discussed below. The shift and the Gaussian broadening σ were then fitted to minimize the differential area fraction A between the theoretical spectrum (I) calculated with Eq. (1) and the experimental data (I exp ) after subtraction of a Tougaard [22,23] background (I bg ): where E 1 and E 2 are chosen to contain all relevant features and I max = 1 as the spectra were normalized. Figure 2 shows the fitted results on isolated monomers using G 0 W 0 @BH-LYP, revealing an excellent agreement with the experiment. The comparison was performed for all possible combinations of three basis sets, eight hybrid functionals, and five backbone corrections (Figures S10-S15). Overall, the best results could be achieved by using hybrid functionals with approximately 50% of exact exchange, such as BH-LYP (50%), B2-PLYP (53%), and M06-2X (54%), similarly to the results of previous benchmarks on gasphase molecules, where an optimal fraction of 45% was found [17]. Remarkably, the isolated monomer approximation is in general sufficient to reproduce the XPS profile of the full polymer. For the two polymers under study, improving on the basis set seems to be a better strategy than improving on the molecular cut, e.g. def2-TZVPP on isolated monomers yields a better agreement with the experiment than def2-SVP on dimers or trimers.
Before extending the methodology to copolymers, we briefly comment on the origin of the shift ∆ between the theoretical and experimental BEs. First of all, we would like to emphasize that most of the chemical information in XPS spectra, i.e. the fingerprint of a material, is retrieved from relative, not absolute, binding energies. Therefore, the presence of a rigid shift, although undesirable, does not necessarily limit the applicability of the theoretical method. The results of this work show that ∆ can be any-where between −2 eV and +8 eV, depending on the polymer under study (PMMA or PHS), on the core orbital (C 1s or O 1s), and on the level of theory, particularly the choice of basis sets and hybrid functionals for the underlying DFT calculations (Figures S10, S12 and S14). Although an absolute shift of 8 eV is in general not negligible, in comparison with the absolute BE this corresponds to a relative error of only 1.5% for O 1s, and up to 2.8% for C 1s.
The presence of a shift can be attributed to several factors on both theoretical and experimental sides. On the computational side, (i) the use of incomplete basis sets (def2-SVP, def2-TZVPP, and def2-QZVPP), and (ii) the amount of exact exchange in hybrid functionals (10% to 50%) can induce shifts of up to 2 eV and 6 eV, respectively. For solid materials, another source of discrepancy is (iii) the inconsistency between the zero energy references in theory and experiment, corresponding to the vacuum level and to the Fermi level, respectively. This inconsistency introduces a shift equal to the work function of the material. For non-conductive polymers as those investigated in this work, additional factors such as (iv) calibration issues due to the fact that the Fermi level is not visible in the spectra as it lies within the band gap, and (v) charging effects in the sample can contribute to shifts of several eVs.
In principle, (i) and (ii) can be tackled by extrapolation to the complete basis set limit and by using more expensive self-consistent GW approaches. When only small gas-phase molecules are considered, solving (i) and (ii) can effectively reduce the overall shift to only a few tenths of eV, as reported by Golze et al. [17]. However, when treating solids, the additional factors (iiiv) will overshadow any attempts to eliminate (i) and (ii). In fact, referencing the BEs to HOMO energies would ideally solve (iii) and (iv), but this is not commonly done on the experimental side, also because of the decreased sensitivity of XPS in the valence region. Moreover, (v) cannot be easily eliminated. The calibration issues (iv) and (v) have been known for decades in the XPS community and have been historically tackled by shifting the experimental BEs to set the main C 1s peak at 285 eV. This approach, although very practical, introduces in fact an additional unknown shift (vi), which is what ultimately makes it impossible to retrieve the absolute BEs from most published datasets. There-  [20] and [21] after subtraction of a Tougaard [22,23] background. The numbers in red identify sets of core-electron levels that are almost degenerate in energy and their location(s) on the molecule. Spectra were calculated on isolated monomers with G 0 W 0 @BH-LYP/def2-TZVPP. Fitted values: σ = 0.61, 0.63, and 0.55 eV, ∆ = 6.13, 5.58, and 5.71 eV, A = 0.014, 0.018, and 0.012, respectively; fit ranges: 291-282, 536-530, 289-282 eV. The peak at 291.4 eV in (c) is a π → π * shakeup satellite [21] and cannot be modeled with the technique presented here.
fore, for non-conductive polymers, only theoretical relative, not absolute BEs can be compared with experiments.
So far, we discussed the case of simple polymers. However, the same additive approach can, in principle, be applied also to copolymers, blends of polymers, as well as mixtures of polymers with non-bonded molecules. As explained above, in all these cases, we sum the spectra calculated for each isolated component, after multiplying the intensity by the molar ratio. To investigate the validity of our methodology for macromolecules consisting of multiple building blocks, we consider a poly[(t-butyl methacrylate)-co-(phydroxystyrene)] copolymer (PBMA-PHS) with monomer ratios w P BM A = 52% and w P HS = 48%. The fitting procedure was repeated using Eq. (2) considering a Tougaard background model [22,23]. Once again, the results are in good agreement with the experiment, as shown in Figure 3. The full benchmark results for different basis sets, hybrid functionals, and backbone corrections are reported in Figures S16 and S17.
Finally, we show that our methodology can also be used to extract the monomer ratio from the experimental spectra, with an error of about 10%. The fit was repeated by optimizing the parameter w P BM A (with w P HS = 100% − w P BM A ) as well as ∆ and σ. The results for different levels of theory are reported in Figure S18. Overall, hybrid functionals with around 50% exact exchange yield fitted PBMA ratio within 10% from the nominal ratio. However, the extracted ratio is very sensitive to the choice of the background model, resulting in a variation of about 10% between linear [22], Shirley [22,24,25], and Tougaard [22,23] models (Table S5 and Figure  S19). Further benchmark studies involving different copolymers and a range of monomer ratios are still required to assess the range of validity of this fitting procedure. Overall, our bottomup approach can be a complementary method to standard fitting procedures to retrieve quantitative chemical composition from experimental spectra.
In summary, we showed that accurate theoretical XPS spectra of macromolecules can be calculated using G 0 W 0 by decomposing a macromolecule into tractable building blocks, such as isolated monomers, and summing up all contributions, weighted by their molar ratio. The calculated spectra are in excellent agreement with the experiments, provided that (i) hybrid functionals with around 50% of exact exchange are used for the underlying DFT calculations, and (ii) a rigid shift is applied, to account for the inconsistent energy reference between theory and experiments. Interestingly, the isolated monomers approximation is sufficient to reproduce XPS spectra of common acrylate-and styrene-based polymers and copolymers. The method may also work for other non-conjugated polymeric systems, pep- Theoretical spectra of a poly[(tbutyl methacrylate)-co-(p-hydroxystyrene)] copolymer (PBMA-PHS) fitted to experimental data from this work after removal of a Tougaard background [22,23]. The spectra were calculated on isolated monomers with G 0 W 0 @BH-LYP/def2-TZVPP and nominal monomer ratio. Fitted values: σ = 0.59 eV, ∆ = 6.17 eV, A = 0.011. The peak around 291 eV is a π → π * shakeup satellite and cannot be modeled with the technique presented here. tides, mixtures of molecules, and other solid materials exhibiting weak inter-monomer and intermolecular interactions. Due to the strong localization of core orbitals, the additive approach might in principle work also in combination with non-GW methods, although we haven't tested this possibility. However, GW is in general preferable to cumbersome methods such as ∆SCF and ∆CC, as GW does not require the explicit creation of core holes. Overall, G 0 W 0 seems to be the best approach in terms of accuracy and computational cost: to give an example of the efficiency of G 0 W 0 , consider that computing the BEs of the isolated monomers of the four polymers in Scheme 1 with G 0 W 0 @BH-LYP takes less than 3 min with a def2-SVP basis set, and up to 50 min with a def2-TZVPP basis set on a single core. Although not discussed in this work, our ongoing calculations show that this approach is also promising for mixtures of polymers with solvent-like molecules. Further studies are still needed to assess if the method can be extended to conjugated polymers, and how impor-tant screening effects are in extended systems, e.g. by comparing the proposed approach with electrostatic embedding methods. In the future, the strategy presented here could be applied to interpret XPS spectra of unknown materials and to study chemical reactions, by comparing the simulated spectra of various candidate compositions with the experimental one.

Experimental details
A 50 nm thick film of poly[(t-butyl methacrylate)-co-(p-hydroxystyrene)] copolymer (52% PBMA, 48% PHS) was spin coated on a 5 × 5 cm 2 silicon wafer in a clean room environment followed by a post application bake at 90 • C for 60 s.
The XPS measurements were performed in a VersaProbe III instrument from Ulvac-PHI using a monochromatized Al Kα (1486.6 eV) pho-ton beam. The sample was kept in the vacuum chamber for one night to ensure that all material outgassing was complete and that the vacuum level in the chamber was sufficiently low (in the range of 10 −8 Pa) to carry out XPS. To reduce the charging on the sample surface and minimize possible degradation due to X-rays, a large exposure spot size of 1000 × 500 µm 2 was used. The measurements were performed at a takeoff angle of 45 • .