Frontier Orbitals and Quasiparticle Energy Levels in Ionic Liquids

Room temperature ionic liquids play an important role in many technological applications and a detailed understanding of their frontier molecular orbitals is required to optimize interfacial barriers, reactivity and stability with respect to electron injection and removal. In this work, we calculate quasiparticle energy levels of ionic liquids using first-principles many-body perturbation theory within the GW approximation and compare our results to various mean-field approaches, including semilocal and hybrid density-functional theory and Hartree-Fock. We find that the mean-field results depend qualitatively and quantitatively on the treatment of exchange-correlation effects, while GW calculations produce results that are in excellent agreement with experimental photoelectron spectra of gas phase ion pairs and ionic liquids. These results establish the GW approach as a valuable tool for understanding the electronic structures of ionic liquids.


I. INTRODUCTION
Room temperature ionic liquids (RTILs) are salts formed of molecular cations and anions that exist in the liquid state at or near room temperature. They find widespread use as solvents, [1,2] dispersants, [3][4][5] and electrolytes, [6][7][8] and exhibit several unusual properties, including high electrochemical stability windows [9,10] and very low equilibrium vapour pressures. [11,12] From a fundamental point of view, it is important to understand the character of the frontier molecular orbitals and determine their quasiparticle energy levels in RTILs, as these determine technologically important properties such as band alignment at interfaces, reactivity, and stability with respect to electron injection or removal. In recent years, a number of experimental and theoretical investigations of the electronic structure of ionic liquids have been reported. For example, photoelectron spectroscopy has been used to study the valence band electronic structure of liquid RTILs and ionic liquid vapours consisting of neutral cation-anion pairs. [13][14][15][16][17][18][19] Computational studies of RTILs have mostly been carried out in the framework of density functional theory (DFT). [13,[19][20][21] Advantages of DFT include its relatively modest computational cost and its ability to predict ground state geometries with good accuracy, as long as dispersion interactions are taken into account. [22] DFT is also often used to gain insights into the electronic structures of materials by comparing Kohn-Sham (KS) eigenvalues to measured photoelectron spectra. However, KS eigenvalues cannot be rigorously interpreted as quasiparticle energies (with the exception of the energy of the highest occupied molecular orbital (HOMO) [23]) * j.lischner@imperial.ac.uk which are measured in photemission spectroscopy. This is the origin of the famous band gap problem of DFT. [24] True quasiparticle energies can be obtained from Green's function techniques, such as the GW approach. In the GW approach, the one-electron Green's function G is obtained by solving the Dyson equation with a self-energy which is given by the product of the Green's function and the screened interaction W . In principle, the GW self-energy should be evaluated using the fully interacting Green's function and screened interaction. In practice, however, a mean-field Green's function G 0 and a mean-field screened interaction W 0 obtained from a DFT or Hartree-Fock (HF) calculation are often used. This approximation, termed G0W0, has been demonstrated to produce highly accurate quasiparticle energies for a wide range of materials. For example, previous work has shown that G0W0 calculations can predict band gaps in solids and first ionization energies of small molecules with high accuracy. [25][26][27][28][29][30][31] Similarly, G0W0 yields accurate results for the position of the d-bands in noble metals relative to the Fermi level, [32][33][34] whereas standard DFT functionals do not. A downside of the G0W0 method is that the results can depend on the mean-field starting point. To overcome this problem, partially and fully selfconsistent GW schemes have been introduced. [35][36][37] In this work, the GW method is used to study the electronic structures of room temperature ionic liquids (RTILs). As a case study, the electronic structure of the 1-Ethyl-3-methylimidazolium tetrafluoroborate ([EMIM][BF 4 ]) ion pair is analyzed in detail with a focus on the nature of the frontier molecular orbitals in this system. Calculated quasiparticle energies from G0W0 calculations are also compared against recent photoemission measurements of several different ionic liquids. In particular, gas phase spectra of ionic liquid vapours are compared against simulated spectra of free ion pairs, and liquid phase spectra of RTILs are com- pared against theoretical calculations of periodic crystalline RTILs. In all cases, excellent agreement between measured photoemission spectra and GW calculations is found, while DFT results depend sensitively on the treatment of exchange-correlation effects.

II. RESULTS
We first consider the electronic structure of the [EMIM][BF4] ion pair ( Fig. 1). Fig. 2 shows the calculated densities of states (DOS) of the [EMIM][BF4] ion pair from different levels of theory. The leftmost column contains results from three different mean-field methods: Hartree-Fock, DFT with the PBE0 functional, [38] and DFT with the PBE functional. [39] The Mulliken decomposition [40] of the total DOS into cation and anion contributions is also shown. The three curves exhibit significant quantitative and qualitative differences. For example, whilst PBE predicts that the HOMO orbital is centered on the anion, PBE0 and HF place the HOMO orbital on the cation and the associated HOMO energies differ by several electron volts among the different approaches. To illustrate this point further, isosurface plots of the HOMO-1, HOMO, and LUMO orbitals are shown in Fig. 3. The three leftmost columns show that PBE, PBE0 and HF predict three different sets of frontier orbitals in this system. In particular, all three frontier orbitals are localzied on the [EMIM] ion in HF, while the HOMO-1 in PBE0 and PBE is on the [BF 4 ] ion. In PBE, the HOMO is also localized on the [BF 4 ]. These results clearly demonstrate that standard mean-field methods are not able to unambigiously answer questions about the nature and energies of the frontier molecular orbitals in the [EMIM][BF 4 ] ion pair.
Calculated densities of states from G0W0 and eigenvalue self-consistent GW (evSCGW) calculations are shown in the middle and rightmost columns of Figure 2. Already at the G0W0 level, the dependence on the mean-field starting point is significantly reduced and all G0W0 results predict that the HOMO orbital lies on the [EMIM] cation. The starting point dependence is even weaker in the evSCGW results. The frontier orbitals from the GW calculations are shown in the rightmost three columns of Fig. 3 and are in qualitative agreement with each other. In particular, all frontier orbitals are localized on the [EMIM] ion. Note that the meanfield wavefunctions are not updated in either G0W0 or evSCGW, explaining the different shapes of the LUMO state of evSCGW@HF compared to evSCGW@PBE and evSCGW@PBE0. Instead, changes in the frontier orbitals arise due to changes in the energy ordering of the one-electron eigenstates when the eigenvalues are recalculated using the GW method. In [EMIM][BF 4 ], G0W0 and evSCGW change the ordering of the frontier orbitals when using a PBE or PBE0 starting point, but not when using a HF starting point.
It is also instructive to consider the absolute energy levels of the frontier orbitals. The calculated energies of the HOMO and the LUMO of the [EMIM][BF 4 ] ion pair from different levels of theory are given in Table  I. The HOMO energies from different mean-field approaches differ by almost 4 eV with HF giving the lowest value (-10.60 eV) and PBE giving the highest (-6.80 eV). This spread is significantly reduced by the one-shot G0W0 correction with G0W0@HF still giving the lowest value (-11.03 eV) and G0W0@PBE giving the highest (-10.44 eV). Eigenvalue self-consistency does not change the G0W0@HF result, but shifts the G0W0@PBE result down by 0.6 eV. Considering next the LUMO level, we find that HF predicts a positive LUMO energy and therefore an unbound state, while the LUMO is bound in PBE and PBE0. G0W0 and evSCGW calculations confirm that the LUMO level indeed lies above the vacuum level and is unbound. Finally, we also compare the mean-field and GW results to ∆-self-consistent-field (∆SCF) calculations, see Table I. The ∆SCF method has been previously used to predict electrochemical stability windows in ionic liquids. [20,41] We find that in the [EMIM][BF 4 ] ion pair, like G0W0 and evSCGW, ∆SCF calculations predict that the LUMO lies above the vacuum level. If PBE0 is used as the mean field theory, G0W0, evSCGW and ∆SCF yield similar results for both the HOMO and the LUMO energies. This agreement indicates that PBE0 is a reliable mean-field starting point for GW calculations in these systems.
We next compare GW results for different ion pairs to experimental photoelectron spectra of ionic liquid vapours. The simulated spectra are constructed from G0W0 calculations with a PBE0 starting point based on the "Gelius approximation", i.e. the spectrum is a sum of atomic orbital projected density of states (pDOS) curves, each weighted by the per-electron photoionization crosssection of that subshell at the relevant photon energy. [42] Uniform Gaussian broadening has been applied to each theoretical spectrum.
Experimental and theoretical gas phase spectra of the 1-Ethyl-3-methylimidazolium trifluoromethanesulfonate ([EMIM][OTf]) and 1-Ethyl-2,3-dimethylimidazolium bis(trifluoro-methylsulfonyl)imide ([EMMIM][NTf 2 ]) ion pairs are shown in Figure 4. In both cases, excellent agreement between theory and experiment is observed. We emphasize that no shifts or calibrations of any kind have been applied to the theoretical spectra, i.e. both the absolute and the relative binding energies of valence electrons in these ion pairs are predicted with excellent accuracy by the G0W0@PBE0 approach.
Interestingly, the agreement between the experimental photoelectron spectrum of vaporized [EMIM] [BF 4 ] and the G0W0@PBE0 result for the free ion pair is somewhat worse than for the ion pairs discussed above, see Figure 6. In particular, peaks A, B, and D' are missing from the simulated spectrum, and the intensity ratios of peaks D, E, and F are different from the experimental ones. In previous studies, it has been observed that the two ions of the [EMIM][BF 4 ] ion pair can react to form an adduct upon heating. [44,45] Figure 5 shows the structure of the adduct. To assess if adduct formation is responsible for the differences between the simulated and the measured spectra, we performed GW calculations on the adduct. We then added the adduct spectrum to the ion pair spectrum assuming that the vapour is a 1.5:1 mixture of ion  pairs and adducts. Figure 6 shows that the resulting spectrum is in much better agreement with the measurement. In particular, peaks B and D' are present and the intensity ratios of peaks D, E, and F are correct, but peak A is still missing. Including eigenvalue self-consistency or using a different mean-field starting point was also not found to reproduce peak A; see the Supplementary Information. We therefore hypothesize that this missing peak originates from a different decomposition product or an ion pair dimer. Finally, we also carry out GW calculations of ionic liquids in the condensed phase and compare them to experimental photoelectron spectra. In principle, the simulated spectrum of the ionic liquid should be obtained by averaging results of different liquid configurations. However, performing many GW calculations of large unit cells is computationally extremely challenging. Instead, we instead carry out GW calculations on ionic liquids in a solid, crystalline phase. This approximation is justified as the internal structures of the ions and their average coordination environments are similar in the solid and the liq- uid phases. Figure 7 shows the unit cells of the three crystalline ionic liquids 1-Butyl-3-methylimidazolium hexafluorophosphate ([BMIM][PF 6 ]), [EMIM][BF 4 ], and 1-Butyl-3-methylimidazolium chloride ([BMIM]Cl) and also compares the simulated G0W0@PBE0 spectra to experimental photoelectron spectra taken in the liquid phase. In photoemission measurements of liquids and solids, the experimental binding energies are given relative to the Fermi level, but since the position of the Fermi level relative to the band edges is not known a priori, the calculated spectra have been shifted by a constant amount to best match the experiment. Excellent agreement between theory and experiment is found for [ [15]. In the spectrum of [BMIM]Cl, the separation between the two most intense peaks (peaks I and IV) is overestimated by approximately 1.3 eV, but otherwise the measured spectrum is reproduced with good accuracy. The results shown in Figures 4 and 7 demonstrate that the GW method is very well suited for predicting quasiparticle energy levels in ionic liquids and free ion pairs.

III. DISCUSSION
An alternative method for modelling photoelectron spectra of ionic liquids based on DFT was proposed in reference [46]. In this study, it was shown that experimental spectra of liquid RTILs can be reconstructed from DFT partial density of states (pDOS) curves of free ion pairs by shifting the cation and anion pDOS curves relative to each other by an amount that is determined on a case-by-case basis. The size of these shifts was originally interpreted as the difference between the average electrostatic potentials experienced by the cation and the anion. Our GW results, however, suggest that this interpretation needs to be revised. In particular, Figure 2 shows that the GW self-energy corrections give rise to a significant relative shift of the anion and cation pDOS curves. This shift does not arise from electrostatic effects, but instead from a more accurate treatment of exchange and correlation effects. Therefore, the shifts applied in reference [46] do not arise solely from electrostatic effects and should be interpreted as empirical corrections that contain contributions from both self-energy effects and changes in average electrostatic potential.
Detailed knowledge of the character of the frontier molecular orbitals and their quasiparticle energies is of crucial importance for understanding the electronic structures of ionic liquids. In this study we have shown that interpreting DFT Kohn-Sham eigenvalues as true quasiparticle energies can lead to qualitatively and quantitatively inaccurate results. This problem can be overcome by the GW method which produces results that are in excellent agreement with state-of-the-art photoemission data. These results suggest that the GW method is a useful tool for studying the electronic structure of ionic liquids and can be used to gain insights into electronic properties that are relevant to ionic liquid devices, such as band alignment at interfaces and stability with respect to electron injection and removal.

IV. METHODS
All calculations reported in this work were performed using the FHI-aims electronic structure program, [47][48][49] that uses atom-centered local basis functions defined on a numerical grid. The geometries of the free ion pairs were relaxed using DFT with the PBE0 exchange-correlation functional until the forces on the atoms were less than 0.005 eV/Å. Van der Waals interactions were accounted for using the Tkatchenko-Scheffler method. [50] For each ion pair, a number of different configurations were manually constructed, and in the end the relaxed geometry with the lowest energy was used for the density of states calculations. The default "tight" numerical basis sets were used during the geometry optimizations. For the bulk crystals, the calculations were performed at experimental geometries from X-ray crystallography. [51,52] The implementation of the GW method in FHI-aims is described in reference [53]. The self-energy was calculated on the imaginary frequency axis with 100 frequency points, and the Pade approximation with 16 fitting parameters was used for the analytical continuation of the self-energy onto the real axis. For the G0W0 and evSCGW calculations, the NAO-VCC-nZ basis sets were used (NAO-VCC-4Z for the ion pairs and NAO-VCC-3Z for the bulk solids). [54] All of the occupied and empty electronic states spanned by the basis sets were included in the GW calculations. A graph showing basis set convergence is included in the Supplementary Information. The bulk calculations were performed at the Gamma point only.

V. DATA AVAILABILITY
The structures of all of the ion pairs and solids considered in this work are given in the supplementary information.
VI. ACKNOWLEDGEMENTS J.M.K. and J.L. acknowledge support from EPRSC under Grant No. EP/R002010/1 and from a Royal Society University Research Fellowship (URF\R\191004). Via J.L.'s membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work used the ARCHER UK National Supercomputing Service. I.K. and V.K. acknowledge Estonian Centre of Excellence in Research project "Advanced materials and high-technology devices for sustainable en-