One-dimensional potential for image-potential states on graphene

In the framework of dielectric theory the static non-local self-energy of an electron near an ultra-thin polarizable layer has been calculated and applied to study binding energies of image-states near free-standing graphene. The corresponding series of eigenvalues and eigenfunctions have been obtained by solving numerically the one-dimensional Schr{\"o}dinger equation. Image-potential-state wave functions accumulate most of their probability outside the slab. We find that a Random Phase Approximation (RPA) for the non-local dielectric function yields a superior description for the potential inside the slab, but a simple Fermi-Thomas theory can be used to get a reasonable quasi-analytical approximation to the full RPA result that can be computed very economically. Binding energies of the image-potential states follow a pattern close to the Rydberg series for a perfect metal with the addition of intermediate states due to the added symmetry of the potential. The formalism only requires a minimal set of free parameters; the slab width and the electronic density. The theoretical calculations are compared to experimental results for work function and image-potential states obtained by two-photon photoemission.


Introduction
Graphene layers display a number of interesting properties and potential applications owing to the linearly dispersing bands found near the K point in the Brillouin zone [1,2]. However, in order to get complete characterization of graphene, other regions in the Brillouin zone need to be considered. In particular, unoccupied states in the vicinity of the point can play a significant role in the transport of currents [3]. Indeed, it is well known the importance of conduction band minima valleys located around in ballistic electron emission processes, where currents are injected on a substrate under voltage-dependent matching restrictions (e.g. k -conservation [4]) that are relevant to design field emission transistors [5]. On the other hand, the transport of heat also involves other regions of the Brillouin zone apart from the K point. Recently it has been shown how the thermal conductivity of a few layers of graphene supported on silicon dioxide depends crucially on flexural modes [6] that are sensitive to the electronic structure around [7]. Therefore, the dielectric response of very few layers of graphene is a key physical element to effectively design devices based on graphene. The dielectric response is directly related, and therefore can be investigated by trapping electrons in the region of unoccupied states between the Fermi level and the vacuum level. These states, bound by their self-induced long-range image potential, are called image-potential states [8]. The experimental [9,10] and theoretical [11] studies of image-potential states constitute an ideal probe to better understand the properties of graphene layers.
The image force is a non-local effect asymptotically dominated by correlation effects [12]. In order to study the infinite Rydberg series arising from the image potential, one needs to compute the effective one-dimensional potential, (z), representing the real part of the quasistatic self-energy for an external probe charge. This self-induced potential is a continuous function spanning from inside the material, where it represents the exchange and correlation energy, to the vacuum region, where it should have the correct Coulomb-like asymptotic behavior ∝ − 1 4z . Such a potential can only be obtained from a non-local spatial formalism, since a local approach results in a correlation potential decaying exponentially in the vacuum region, following the density behavior outside the solid [13]. For a self-consistent first-principles theory such a non-local functional dependence requires costly numerical calculations [14]. Therefore, it is useful and natural to search for simpler ways to obtain such an effective potential, which is the basic ingredient needed to understand the physics of image-potential states bound by an ultra-thin polarizable layer like graphene. The simplest of these alternatives is to introduce a set of fitting parameters to continuously join solutions valid either inside or outside the solid. This point of view has been taken, e.g., by Silkin et al to study image-potential states in freestanding graphene, joining a function with the correct classical asymptotic behavior outside to a first-principles calculation inside a graphene layer [11]. Such ab initio calculations depend on choosing a model for the exchange and correlation potential; local density approximation (LDA) has been employed by Silkin et al [11]. Furthermore, a very large unit cell needs setting to minimize effects between charged periodic images (i.e. 85 Å vacuum separator was used to reach convergence, implying a large number of plane-waves and a serious computational effort). Finally, such a calculation needs to be supplemented by a few adjustable parameters, e.g., the choosing of a matching point to join the inside potential to the asymptotic classical potential.
In this paper, we analyze an alternative that makes use of a minimal free parameter set and makes a simple, flexible and accurate basis for interpreting experimental results. We use a well-known model for the reflection of electromagnetic fields at the surface (infinite barrier specular model [15]) coupled with a non-local static dielectric response so the desired selfenergy can be obtained [16]. For the electrodynamics model, two free parameters are introduced, i.e. the electronic polarizability of the thin slab, which is determined by the Fermi-Thomas (FT) screening wavelength inside the slab (electron density of the material), and a geometrical dimension given by the layer thickness. This approach leads in a natural way to a potential with proper physical features: it is continuous and finite over the full spatial domain and it has the right asymptotic behavior toward the vacuum region. Recently, Ghaznavi et al [17] have studied the nonlinear screening of an external charge near a doped graphene layer by solving a FT model via a nonlinear integral equation. The nonlinear image potential shows changes up to 0.2 eV with respect to the classical image potential, and supports the use of the random-phase approximation (RPA) dielectric response for graphene. In the case of a graphene layer lying on a metallic support, another free parameter may be introduced in this model: the wave function penetration in the material. Image-potential states are supported in many metallic surfaces for energies between the vacuum level and the Fermi level due to the existence of surface band gaps that prevent the penetration of the wave function toward the bulk and allow the existence of bound states [8]. Therefore, penetration of wave functions is determined by the band structure of particular materials and surface orientations in a way that cannot be incorporated into our model, except by including a free parameter that globally determines this penetration. Here we shall first discuss two limiting cases: 0-100% penetration of the graphene layer, and then the continuous evolution from one limit to the other. Since for image-potential states we are interested in regions in reciprocal space with k near and energies between the vacuum level and the Fermi energy, it is well justified to model graphene as a polarizable electron gas with quadratic energy dispersion, as seen from the relevant bands of graphite and graphene [11,[18][19][20]. A connection between the interlayer states and the image-potential states in graphite, graphene and other carbon-based materials has been established by different groups, e.g. Silkin et al [11], Feng et al [21] and Hu et al [22]. These authors have pointed out that such states exist universally in two-dimensional materials.
The static dielectric response, ( k), has been modeled by the RPA and by its small k expansion, the FT approximation [23]. While the RPA yields a more accurate description of excitation in the material, introducing the FT allows us to write the potentials as analytic expressions or quasi-analytical ones which merely depend on the final numerical step involving the simple integration of a function decaying quickly for large values of the argument. A single parameter, the screening constant related to the density of states at the Fermi level, k FT ∝ ∂n 0 ∂µ , fixes both the scales for energies and lengths facilitating the rescaling of results for different materials and sizes (atomic units are used throughout the paper, except where explicitly stated otherwise). Taking graphite as a model (2 g cm −3 , 2s 2 2p 2 ), a typical value for graphene is found to be k FT ≈ 1 (r s ≈ 2.5), although its precise value should depend on factors such as doping, external potentials, etc; this is accommodated in our results through the simple aforementioned scaling with k FT . The other parameter used to characterize a thin slab is its width, 2d. For a single atom thick layer of graphene a reasonable value for d should be related to the spatial extension of π carbon orbitals, d ≈ 1 au. The value of this parameter turns out not to be critical for this work because wave functions spread over regions much larger than d.
The theoretical results are compared with the experimental results for graphene on various substrates. The presence of the substrate leads to charge transfer from or to the graphene which can be described also as doping. The resulting work function change due to graphene has been modeled by Giovanetti et al [24,25]. We find good agreement with the experimental data. The doping changes the available screening charge and leads in turn to a change in the energy of the image-potential states. The experimental data from two-photon photoemission (2PPE) are in agreement with the calculated dependence.

Self-induced potential by an ultra-thin slab
For an external probe charge near a slab (Q = 1) we seek the potential acting on Q by the polarization charges induced in the medium by Q itself. This is obtained by computing the total potential and subtracting the charge's own naked potential. To ensure suitable boundary conditions, and according to the specular reflection model at the surface, auxiliary pseudo-media are introduced for the polarizable slab and the vacuum that reduce the calculation to matching solutions obtained in different regions of space for the homogeneous media everywhere [12]. Details of the calculation for the thin slab are given in appendix A. The resulting potential depends on a number of integrals that include the dielectric response of the system; for the RPA these are computed numerically. On the other hand, within the FT approximation an expression that only depends on a single numerical integration can be obtained, , where χ = κ 2 + k 2 FT . This is an useful expression that can be computed very efficiently. We shall see that for the purpose of computing the energy levels of image-potential states it makes an excellent approximation to the more costly RPA calculation.
In figure 1, we show the potential for a slab occupying the region −d z d, both in the FT approximation (black continuous line), and in the RPA one (black dots). In the region determining the Rydberg series (|z| d), both approaches yield similar values and are in agreement with the correct asymptotic power law. Near the center of the slab, FT overestimates the interaction over RPA by about 30-40%, a difference that is reflected mainly in the lowest state (node-less) that has a significant weight in the central part of the slab where the difference between RPA and FT is larger. This state is located outside the window between the vacuum level and the Fermi energy and is not part  (1) and (2) can be evaluated analytically to obtain the following useful particular values (see appendix B): We remark that the latter value corresponds to the Coulomb hole [26]. Furthermore, for z > 1 k FT (vacuum region) the potential is well approximated by a classical law, 1 4(z−z 0 ) , corrected by an image plane, z 0 . The value of z 0 can be obtained by expanding the integrals for κ ≈ 0, which fixes the position of the image plane in this model: figure 1, see appendix B). Finally, equation (3) suggests a simple procedure to find a good approximation to the RPA potential. Given a physically representative value for k FT (r s ), one can find another value, is a good approximation to the full RPA potential that can be computed quickly and easily.

Eigenvalues
We solve numerically the Schrödinger equation [27,28] to compute the eigenvalues and eigenfunctions corresponding to the RPA and FT model potentials described above. Quite generally, bound states can be unequivocally labeled by the number of nodes n, with energies increasing as the number of nodes increases. Furthermore, as long as the potential is symmetric, eigenfunctions have either even or odd parity, for an even or odd number of nodes, n.
We show in figure 1 the first five eigenvalues for (z); dashed (red) and dotted (blue) horizontal lines for even and odd parities. It is worth noting the structure of this series; there is an isolated eigenvalue (the lowest one), while the remaining states cluster near the vacuum level.  (2)) are given for a free-standing graphene ultra-thin layer and compared with similar results obtained from an ab initio LDA calculation matched to an asymptotic expression for the image potential by Silkin et al [11] (states are labeled according to their ordering in the image-potential series of unoccupied states lying above the Fermi energy and below the vacuum level). The effect of a repulsive barrier located at z ≈ −10 Å is presented under R b−20 (notice that in this case, the symmetry is broken and the parity is no longer a good quantum number). For comparison, the experimental values measured on Gr/Ir are quoted [9], and compared with the eigenvalues for the classical Whittaker's problem [8]. Finally, we show the effect of approaching the repulsive barrier to z ≈ 0.5 Å, R b=1 , that can also be computed by introducing an appropriate quantum defect (δ = 0.2). To guide the eye, we highlight in bold face numbers that can be compared across different calculations or experiments and have been given an accompanying interpretation in the text.
For standard densities (e.g. 1 < r s < 10) this first level appears below ∼ −4.5 eV, an estimate for the work function in graphene (thick line). Therefore, this n = 0 eigenstate does not fit in the standard definition for a Rydberg state, necessarily located between the vacuum and Fermi levels to be observable in a standard experiment.
In table 1, we compare eigenenergies calculated for the RPA and FT models with similar theoretical values obtained from a formalism based in density functional theory (DFT) [11], and with experimental values reported for Gr/Ir [9]. We have also included, as a reference, the limiting case of the Rydberg series for a perfect metal E m+1 = − 1 32(m+1) 2 , m = 0, 1, 2, . . ., where m refers to the number of nodes for each state (note that these wave functions only extend to z > 0 half-space, and that the zero at the origin is not counted as a node since it derives from the boundary conditions). Figure 2 shows on a log-log scale the scaling law for the first few members of the series. We remark that for the purpose of computing binding energies for image-potential states the FT model with a suitable value for k * FT yields values of the same quality as those obtained from the RPA at a much lower computational cost. For the potential created by free standing graphene, Silkin et al [11] have computed binding energies for image states by matching an ab initio DFT model for the total energy to an asymptotic expression for the image potential beyond some distance, z 0 . A pseudopotential with valence electrons 2s 2 2p 2 and a pseudoized core for 1s 2 electrons has been used. In order to compare the results from these different theoretical approaches, we align the states attending to their energies in the same order as they appear in the image series. This labeling of states is straightforward and unambiguous, and most importantly it facilitates a comparison between the different theories and with the experimental values. On the other hand, the theoretical eigenvalues can be classified by looking at the properties of their eigenfunctions, notably the number of nodes and for a symmetric potential the parity. Eigenvalues corresponding to a symmetric potential like the one plotted in figure 1 are easily classified in the same order as they appear, n, by the number of nodes (n), and their parity (+ for even n, and − for odd n). However, the parity might become not a good quantum number in the presence of defects or a perturbing potential modifying the symmetry of the potential (e.g. a supporting substrate for the graphene layer), while the ordering dictated by the number of nodes is related to the orthonormality of wave functions and makes a robust labeling scheme. Within this ordering, the first image-potential state is simply the first unoccupied one to lie above the Fermi energy and below the vacuum level determined by the work function, W . In the RPA/FT model we are introducing in this paper, the first image-potential state has been labeled n = 1 in table 1; its wave function has one node and it is antisymmetric, E(1) RPA = −0.68 eV. In the DFT model, the first image-potential state we find is E(1) DFT = −1.29 eV, and its wave function has two nodes and it is symmetric (labeled in [11] as 1 + ). The number of nodes on this wave function has been correlated in [11] with the occupied σ and π states below and the required orthonormality condition between them. The discrepancy between the number of nodes and the parity on wave functions from these two theoretical approaches is not a fundamental one. Presumably a pseudopotential including a different core, like a simpler one constructed with only 2p 2 electrons, or a more complete one including the 1s 2 electrons, would alter the number of nodes in the first unoccupied state because the orthonormality condition for the low lying states would be different. Similarly, in our RPA model one could imagine a scenario, where the value of k FT is so small as to make the n = 0 state the first unoccupied state, or inversely, so large as to get the n = 1 state below the Fermi level to make the n = 2 state the first member of the series of unoccupied states. Therefore, we conclude that the more convenient way to compare different image series from different theories or experiments is the ordering attending their energy, because it does not depend on the fine details of the model being used. This criterion brings good agreement between DFT and RPA/FT, and with experiments, which is remarkable taking into account how different both theories are, which is something to be highlighted. Perhaps the exception is the first value predicted by DFT, −1.29 eV, that is a bit too low. We do not take this small discrepancy seriously since both theories involve a number of approximations and parameters (the pseudopotential and the exchange-correlation potential for DFT; k FT and d for RPA), and it is well known the difficulties in DFT of getting accurate values for empty states (e.g. band gaps) without including a more sophisticated description for electron correlation, which is the physical effect responsible for the asymptotic behavior, and one of the assets of the RPA/FT formalism.

Eigenfunctions
The similarity of the values found for the antisymmetric (odd, n − ) members of the series of states for (z) and the classical Rydberg series is striking, and draws some attention (third row in table 1). Such similarity can be understood by looking at the corresponding wave functions ( figure 3). Although the wave functions for Rydberg states are spatially located mostly in the vacuum region, the symmetric members of the series have ψ(z = 0) = 0 at the origin of the slab (figure 3, the continuous line in upper-left panel). This is most conspicuous for the ground state wave function ψ 0 that is more like the ground state of a harmonic oscillator fitted to the bottom of the potential well ( (z) ≈ −0.5 + 0.13z 2 au, E 0 = −6.2 eV) than to states in the one-dimensional hydrogen-like series for a semi-infinite metal that are assumed to go to zero at the image plane. Boundary conditions force all Whittaker wave functions to go to zero at the origin, a condition that in the case of a symmetric well can only be fulfilled by the odd wave functions. Moreover, if n − and m give the number of nodes for odd wave functions for the symmetric potential and the Rydberg one, respectively, we can make a one-to-one correspondence, n − −1 2 = m, that simply tells us that both sets of wave functions have the same number of nodes if n − is divided by two (only half-space) and the node at the origin is discounted. To compare with wave functions corresponding to Whittaker's classical problem defined only over half the space, one might use Silkin et al's labeling for image-potential states: a quantum number related to the nodes of the wave functions for the symmetrical problem only on half of the space supplied with its parity to take into account that antisymmetric wave functions add an extra node to the count. From a physical point of view, we can envisage two relevant limits: a free standing slab producing a symmetric potential with states labeled by the number of nodes and their parity, and a slab on a substrate where a particular surface gap may prevent penetration of wave functions inside the material, leaving only half the space accessible for image-potential states. In this latter scenario, parity would not be an useful quantum number. It is reasonable to describe the case of a graphene layer grown and supported on a particular substrate as having a place somewhere in between these two limits depending on details such as the interaction between the graphene layer and the support, the electronic surface structure of the combined system, work function, etc. In what follows, we shall assume that all these details can be taken into account in the simplest terms by a single free parameter giving the amount of penetration of wave functions. This electrodynamics formalism is not intended to describe these details, but it can bring useful quantitative information on the observed levels, and as a consequence, further conceptual understanding.
This formalism predicts, for a free-standing graphene layer, the appearance of new states associated with the even parity (e.g. eigenvalues between −0.44 and −0.36 eV for n = 2 and between −0.15 and −0.14 eV for n = 4 in figure 2). These new states have been obtained numerically and fit well into the classical scaling law proportional to n −2 for n up to 7. Obviously, the symmetric potential can be perturbed by external ones (e.g. the cases of graphene on a support), and these states would be affected accordingly.

Modeling the substrate
So far we have discussed a model that effectively represents a free-standing graphene layer. For those cases where the graphene layer has been deposited on a metallic support the substrate is expected to manifest itself in two main physical ways: (i) the wave functions may be constrained to be outside some spatial region where the substrate enforces an electronic gap; and (ii) as a consequence of the interaction between the layer and the support, some charge may be transferred to/from the slab, modifying the density of states at the Fermi level, i.e. the value of k FT .
To assess how sensitive the eigenvalues are to the penetration of wave functions into the material we have added to (z) a repulsive term modeled as an exponential wall, R b (z) = e −a(z−b) . Since we are only interested in creating a decaying state inside the support, we fix the parameter a to a large value, a = 20 au, akin to the infinite hard-wall limit; its effective role is to expel states from the z < b region, ensuring the exponential decay of wave functions inside an electronic band gap. The resulting potential for the repulsive barrier located near the slab surface (b = 1 au, green dashed line in figure 1) is similar to the classical series with an image plane, 1 4(z+z 0 ) , and can be easily solved by introducing a quantum defect in the Rydberg formula (compare energies for the same number of nodes in the eigenfunctions for R b=1 and the Whittaker series with quantum defect δ = 0.2 in table 1). The barrier, on the other hand, can be introduced below the surface, mimicking the effect of an electronic band gap due to a supporting substrate. The evolution of the first few eigenvalues with the position of the barrier has been shown in the two limits in table 1. As long as the barrier is located far away from the ultra-thin slab (e.g. b 20 au) we get values reminiscent from the original members for the unperturbed symmetrical potential. At the other limit, a barrier located just on the surface very much recalls the classical solution. For large m the eigenvalues are determined by the potential in the vacuum region and the existence of such a barrier distorts the states less and less as they approach the vacuum level. Table 2 gives the expectation mean values in Å, z = ∞ 0 ψ(z) z ψ(z) dz, for Whittaker wave functions compared with the ones obtained for the RPA or FT potentials (e.g. (B.2)). These values compare well with each other, which reflects the manifest similarity between wave functions shown in figure 3, and show how the important region for the potential moves quickly away from the surface as m grows. The fact that z d for image-potential states with n > 1 implies that wave functions are quite insensitive to the potential inside or near the layer and they are mostly influenced by the asymptotic region where FT and RPA are equivalent. This suggests that higher k-corrections to the dielectric function arising from the RPA are not very important, at least for n > 1 states. Taking away the first level, largely affected by the details near the bottom of the potential, the rest of the series is only modified by a percentage comparable with the differences found in table 1 between similar entries.
The effect of doping may be explored by looking at the k FT dependence of the energies of the image series. The important parameter of the dielectric model is the charge density of the graphene layer, which can be related to the doping level and the work function. This is summarized in the current formalism via a single parameter, the FT wave-vector k FT . While k FT cannot be easily extracted in our formalism from doping levels of graphene, we can compare with experiments by exploiting the linear dependence predicted by this theory between the screening wave-vector and the work function of the thin slab. To compare with available experimental results, it suffices to establish a connection between k FT and a relevant energy scale, e.g. the distance between the vacuum level and the Fermi level, i.e. the work function, W . To this end, we note that in our approach the work function for a semi-infinite surface, or a slab with d 1 Å, is − k FT 2 in the FT model and ≈ − k FT 3 in the RPA one. We can see that RPA, so far our best approach, overestimates W when compared with experimental values by about a factor ≈ 2. The linear dependence between W and k FT , on the other hand, is solidly anchored in the theory, having its origin in the net attractive interaction between the external electrons  [9,29,30]. For Ru the data [10] are plotted by filled circles (green) for two different work functions: 4.00 and 4.24 eV. Data for graphite have been taken from [31]. and the polarization charges created inside the polarizable material. Therefore, W should be proportional to αk FT , albeit the proportionality constant should be corrected to compare with the experimental value. We take α = 1 6 , which corresponds to the empirical work function for graphene, W 0 = 4.5 eV. Therefore, the empirical dependence W = k FT 6 allows us to translate k FT into an experimental energy scale, and to set up an energy origin at the value W 0 . Figure 4 shows the results for the first three members of the series by straight lines. The experimental data will be discussed in section 3.
An interesting feature of this result is the fact that the higher members of the series are less and less affected by changes in k FT (W and/or ρ(E F )). This is clearly seen in figure 4 where the slope of ∂ E n (W )/∂ W decreases for the higher n values. In figure 5 we plot this result, and show that it is well fitted by the empirical function 1 6n 2 . Since k FT is linearly proportional to the density of states at the Fermi level, ρ(E F ), this gives us the empirical scaling law expected for the different terms of the image-potential state series with doping.

Experimental results
The theoretical results of the preceding section can be tested experimentally by 2PPE. This technique is able to measure image-potential states with high accuracy and results for various graphene-covered surfaces have been reported [9,10,29,30]. Photoelectron spectroscopy can also provide precise values for the work function of the surfaces, which is correlated with the doping level and charge density of the graphene layer. We first discuss the work function before turning to the image-potential states.

Work function
The amount of charge transfer from the substrate to the graphene layer is determined by the work function difference between the substrate and graphene. This situation has been modeled with a simple capacitor model by Giovanetti et al and validated by comparing with results from calculations for various surfaces [24,25]. The curve plotted in figure 6 shows the calculated work function of the graphene-covered surfaces versus the work function of the clean substrates for the capacitor model. The calculated values plotted by green open squares fit the curve quite well [24,25]. The available experimental data are shown by blue open circles [9,10,29,30,32,33]. The size of the circles represents approximately typical experimental error bars. The curve was fitted to the experimental results for the noble-metal surfaces [29] and monolayer graphene on SiC [30]. The fit was done with the work function of graphene of 4.50 eV compared to 4.48 eV in the original work [24,25]. The chemical shift was reduced from 0.90 to 0.89 eV. Overall the fit of this work describes also the calculated values (green open squares in figure 6) very well. The surfaces included in the fit all have a graphene-metal distance of around 3.4 Å [34][35][36][37][38]. The surfaces of Ni, Ru and Pd were not included in the fit and are shown by blue dashed circles. They have shorter graphene-metal distances [24,25,39] and would require a larger chemical shift for a satisfactory description. The capacitor model provides a meaningful measure to compare different surfaces with comparable graphene-metal distances via the work function. In addition, the work function for graphene-covered surfaces can be estimated from the work function of the substrate at least for weakly coupled graphene.
The charges transferred between the substrate and the graphene layer originate from the Dirac cone on the graphene side. These might be electrons or holes depending on the doping of the graphene. The shift of the work function due to the charge transfer can therefore be directly  [36,[40][41][42][43] and are plotted in figure 6 as filled red circles. The agreement is perfect for SiC, but for the other surfaces larger deviations are observed. One has to keep in mind that for the noble-metal substrates the graphene layer is p-doped and the Dirac point cannot be observed directly in photoemission spectroscopy. Its energy is extrapolated from the dispersion of the Dirac cone below the Fermi energy. The experimental results are in good agreement with the calculations, which show the constant difference between the work function and the energy of the Dirac point as predicted by calculations [24,25]. This has been found also for different modifications of graphene on SiC [30].

Binding energies
The important parameter of the dielectric model is the charge density of the graphene layer, which can be related to the doping level and work function. The energies of the image-potential states as a function of work function have been plotted as red solid lines in figure 4. Large open circles (blue) show the available experimental values [9,29,30]. The symbol size represents typical experimental uncertainties. The data points for graphite were taken from [31], which observed the lowest three image-potential states. Other groups reported binding energies of 0.85 ± 0.10 eV for the n = 1 image-potential state [44,45]. The experimental values are slightly above the calculated lines, but the slope is in fairly good agreement. The agreement could be improved by using a smaller value for the parameter α relating the work function and FT wave vector, W = αk FT . The data point for the n = 1 state of SiC lies significantly too high. This might be an effect of different screening properties of the dielectric substrate or residual binding to the buffer layer. The accuracy and work function range of the experimental data for the noble-metal surfaces is not sufficient to determine the slope independently. Therefore, we only conclude that the experimental data are in reasonable agreement with the predictions of the dielectric model calculations.
Finally, we discuss the interpretation of experimental data measured on Gr/Ru [10]. Three bound states have been measured with respect to the Fermi energy (all energies in eV): 3.44 (1 ), 3.59 (1) and 3.82 (2). The work function was measured as 4.24 eV. On the corrugated graphene on the Ru (0001) surface also lower areas exist which have a work function of 4.00 eV. On all other surfaces only the n = 1 and 3 image-potential states have been found. It is therefore worthwhile to check whether on Ru the additional state might be the n = 2 image-potential state. The data are plotted as green solid circles with error bars in figure 4 for the two different values of the work function. For a work function of 4.00 eV, the experimental binding energies are in reasonably good agreement with the calculated lines. This assignment would also be compatible with the observed monotonic decrease of the lifetime with the binding energy [10]. However, we cannot rule out the consistent interpretation based on different local work functions on the corrugated graphene on Ru as proposed by Armbrust et al [10].

Conclusions
Using standard models for the dielectric response and the reflection of electromagnetic waves at a surface we have computed the static self-energy for an ultra-thin slab mimicking a graphene layer. The self-induced potential goes continuously from the exchange and correlation energy inside the material to the classical asymptotic image potential in the vacuum. For the purpose of obtaining image-potential state binding energies we find that FT makes an excellent and convenient approximation to the accurate RPA. Eigenvalues and eigenfunctions have been compared with the Whittaker classical series and recent experiments on Gr/Ir. A free standing graphene ultra-thin layer produces a spatially symmetric self-energy that induces an image potential series with even and odd states. The odd members of the series show a remarkable resemblance to the solution of the Schrödinger equation for the classical image potential (Whittaker wave functions). On the other hand, even wave functions arise as new states that differ from Whittaker in several key respects, e.g. their non-zero density probability at the origin. While the qualitative aspects of the image series supported by a thin slab can be described in terms of quite general physical properties, the detailed quantitative values depend crucially on the penetration of wave functions into the substrate and the thin graphene film supported on it. We have considered the limiting cases of full and no penetration. In the case of films weakly interacting with a support and wave functions penetrating well inside the system, some new states may consequently appear in between the classical ones that can be traced back to the even states in a free-standing slab. In cases where the interaction is strong and the surface electronic structure prevents the penetration of wave functions, behavior more similar to a standard metallic surface is expected. We note that the formalism used here can be easily generalized to include the effect of surface plasmons via a frequency-dependent dielectric function [46].
The experimental results for graphene on various substrates compare well with theoretical predictions. The measured work function change due to graphene is in agreement with the capacitor model of Giovanetti et al [24,25]. This opens the possibility to predict the work function of graphene-covered surfaces from the substrate work function. The work function difference from an isolated graphene sheet is related to the doping of the graphene layer and determines in turn the available screening charge. The resulting change in the energy of the image-potential states calculated with the theoretical model is in agreement with the experimental data from 2PPE. Therefore, the total potential for the symmetric (antisymmetric) cases are (±) n e iq(z 1 +2dn) ± e iq(2d+z 1 +2dn) . (A. 2) The surface charges, σ S and σ A , are separately obtained from the condition Convergence of the different series in these expressions is guaranteed by the zeros of the dielectric function, i.e. the normal modes dressed by the interaction in the medium [12]. For the simple Thomas-Fermi dielectric function we have k 2 FT (k) = (q 2 + κ 2 + k 2 FT ) = (q + iχ)(q − iχ), and the integration over the perpendicular moment, q, can be performed to obtain analytic or semi-analytic expressions that are given below.

A.2. Q inside the vacuum (I and IV, |z| d)
A similar procedure yields for the case of Q in the vacuum region Equation (A.1) yields a numerical procedure that allows us to obtain the self-induced potential. Using the FT approximation, quasi-explicit expressions that only depend on a single numerical integration have been obtained and given in (1). We remark that this is an excellent approximation to the full RPA result for the purpose of computing binding energies of imagepotential states.

Appendix B. Useful analytical results in the Fermi-Thomas approximation
The semi-infinite system (d → ∞), along with the use of the FT dielectric function, brings some simplifications to the expression for the potential that can be exploited to compute the exact values at the surface, well inside the material and in the asymptotic vacuum region. For convenience, in this appendix we move the origin to the surface (z = d).

B.1. Induced potential in the vacuum region
Using the equations in appendix A we obtain the induced potential in the vacuum (now z > 0): that can be further simplified by the use of k 2 (k) = k 2 FT (k) = κ 2 + q 2 + k 2 where 0F1 (; b; z) is the regularized confluent hypergeometric function, H 2 (z) is the Struve function and J n (z) are Bessel functions of the first kind.

B.2. Asymptotic behavior and position of the image plane
We can study the asymptotic behavior in the long wave-length region (κ ≈ 0) by expanding the fraction in the integrand in equation (B.2), κ−χ κ+χ ≈ −1 + 2κ k FT + O(κ 2 ), to obtain

B.4. Potential in a vacuum gap
It is handy to apply the same techniques to the inverse problem: the potential in a vacuum gap between two semi-infinite media at z < ±d. The theoretical procedure proceeds along similar lines as we have analyzed in this paper, and we simply give the FT result here  .
(B.9) ( 0 (z) is the digamma function and γ = −0.577 is the Euler constant), (ii) the classical value corrected by an image plane determined by expanding the full non-local potential calculated in the FT approximation and (iii) the independent non-local potential for each of the two surfaces delimiting the vacuum gap.