Conference-ISSS-8-A One-Dimensional Model for Photoemission Calculations from Plane-Wave Band Structure Codes

A computational method is devised for calculating angle-resolved photoelectron spectra using a plane-wave basis and repeated-slab geometry. The method is tested with a one-dimensional model with a rectangular potential well. The model parameters are adjusted to reproduce electronic structure of graphene at the Γ point of the Brillouin zone, as obtained from density functional theory. Photoemission final states with proper boundary conditions are constructed from linear combinations of supercell eigenstates such as to match the free wave at the center of the vacuum region. These states, calculated with a moderate plane-wave cut-off, agree very well with the exact wave function. The computed photoemission intensity differs strongly from the popular (single) plane-wave approximation and shows a pronounced energy dependence. [DOI: 10.1380/ejssnt.2018.49]


I. INTRODUCTION
Angle-resolved photoelectron spectroscopy (ARPES) is the most direct method for probing the band structure of crystalline solids [1].It can also be used to measure wave functions of molecules through the so-called orbital tomography method [2].While the theory of photoemission is well-established [3], computation of ARPES spectra remains a challenge for many systems.In practice, the final state plane wave (PW) approximation [2,4] is often used despite its known shortcomings [5].Photoemission final state calculations beyond the PW approximation have been done with multiple scattering and band structure methods.For metals, good results are obtained with multiple scattering theory, where the proper photoemission boundary conditions can be easily implemented, either in the layered KKR approach [3,6] or in a fully real-space scheme [7].However, the multiple scattering method is not sufficiently accurate for molecules and other covalently bonded open structures including light elements.ARPES can also be calculated by matching bulk bands to the free electron waves in vacuum via some surface potential model [8].This approach yields also good results for metal surfaces [9,10] but it seems difficult to generalize it to complex surfaces and low-dimensional systems.On the other hand, the supercell (or "repeatedslab") method with a PW basis set, can describe accurately the electronic structure of very diverse systems, including molecules, complex surfaces and adsorbates.The supercell approach has been applied to ARPES only very recently by Kobayashi et al. [11] who computed the spinpolarization of low energy ARPES from a Bi(111) surface.
Here we outline a method for obtaining final state wave functions with the proper photoemission boundary condi-tions.Our method is similar to that of Ref. [11], but there are some differences which we shall discuss below.The final states waves are expanded over a set of supercell eigenstates with different perpendicular momentum and with energy approximately equal to the chosen photoelectron energy.These waves are matched with the free photoelectron wave at the center of the vacuum region, i.e. between repeated slabs.We test this scheme with a onedimensional model which mimics the electronic structure of graphene at the Γ-point of the Brillouin zone, as shown by direct comparison with band structures obtained from density functional theory (DFT).Using a moderate PW cut-off, we obtain virtually exact wave functions for the rectangular potential well.Photoemission from the highest occupied band is calculated and found to display a pronounced energy dependence.The computed ARPES spectrum differs dramatically from that obtained in the (single) PW approximation, despite the fact that the photoemission final state bands resemble free electron bands.

II. THEORY
In the independent particle approximation, the ARPES intensity for a transition from initial state wave ψ i to final state wave ψ f is proportional to the square of the transition matrix element where A is the vector potential of the light and P is the electron momentum operator.ARPES boundary conditions are such that for points r far from the surface, the final state wave ψ f becomes a PW, ϕ p (r) = exp(ip • r), with energy ϵ p = p 2 /2m, where p is the measured photoelectron momentum.Such a final state with proper boundary conditions will be denoted ψ p .We compute both ψ i and ψ p in a supercell approach with PW basis set.Since the supercell has periodic boundary conditions in all three dimensions, which is wrong for the perpendicular direction, the supercell eigenstates must be matched to the photoelectron PW in the vacuum region.For crystal surfaces the parallel component of the crystal momentum, k || , is conserved in the photoemission process.The perpendicular component k z is not conserved, so we develop the final state over a number of Bloch waves with different k z .In this paper we introduce a matching technique for the one-dimensional case, where the wave functions are only functions of z, and we put k = k z .The general method in three dimensions will be presented in a forthcoming publication.Inside the supercell of length c, the photoelectron wave ψ p (z) is expanded over a set of Bloch eigenfunctions ψ nk of energy ϵ nk , where the prime indicates that the sum is restricted to a few states with ϵ nk ≈ ϵ p .The coefficients α nk are found by matching the wave amplitude and the first derivative of ψ p , to the free wave ϕ p (z) = exp(ipz) at the cell boundary z 0 is chosen at the center of the vacuum region between repeated slabs.In one dimension there are only two matching equations ( 3), and so we include only two band states (nk = 1, 2) in the sum (2).The supercell periodicity introduces small gaps in the band structure, as it is well known from the Kronig-Penney model [12] and confirmed in our calculations below.Moreover, in numerical practice, the Brillouin zone is sampled, and thus the energy spectrum ϵ nk is discrete.As a consequence, energy conservation between inner and outer waves can in general only be satisfied approximately, i.e. the states included in the sum (2), have energies ϵ nk slightly different from ϵ p .Here we select the two band states with energies just below and just above ϵ p .The precision with which the condition ϵ nk = ϵ p can be met, depends on the number of k-points and the width of the band gaps, which scale inversely with the supercell size c.We use a PW basis, so the eigenfunctions ψ nk are written as where G is a reciprocal lattice vector and C nk G are the PW amplitudes.Inserting (4) into (3) yields the following 2×2 linear system for the coefficients α 1 , α 2 of the two selected band states. ( ) The wave function ψ p in Eq. ( 2) with coefficients α 1,2 so obtained, is a continuum wave with proper photoemission boundary conditions.In the application below, we use a simple rectangular potential well and choose parameters appropriate for normal emission from graphene.The present method is similar to the one recently proposed by Kobayashi et al. [11] but there are some clear differences.While Kobayashi et al. determine the Bloch wave coefficients by numerically minimizing the difference between the combined Bloch wave and the plane wave in a finite range of the vacuum region, we perform a direct matching of wave functions on a single plane in vacuum.With our method, the exact wave function can in principle be recovered as we show in Fig. 3 for a simple model potential.In the method of Ref. [11], when the matching region is chosen too small, the result will be inaccurate.If it is chosen very large, the supercell becomes very large and thus the numerical calculation becomes much more heavy.This problem is absent in our method.The matching can be done at any plane in the vacuum region and it does not increase the supercell size.Both methods should have advantages and shortcomings and so it is interesting to pursue both approaches in the future.

A. Density functional theory calculation of graphene
In order to find realistic parameters for graphene, we have performed a DFT calculation in the repeated slab geometry using the plane wave code VASP [13].The supercell lattice parameter c is set to 20 Å, which is large enough to avoid interactions between periodic images and to get free electron-like states at high energy.The origin of the energy scale is chosen as the vacuum level, determined at the center of the supercell between graphene layers.Figure 1  (E > 0) follows closely the free electron result with √ E dependence.This shows that the vacuum is large enough and the supercell continuum states are sufficiently complete.

B. One-dimensional model with rectangular potential well
We test our method with a 1D model and a rectangular potential well.The width (2.1 Å) and depth (−52 eV) of the potential well are chosen such as to roughly fit the k z -bands of the DFT graphene calculation at the Γpoint.The model is solved by numerical diagonalization of the hamiltonian, with a PW basis and periodic boundary conditions with a period of c = 20 Å as in the DFT calculation.The PW cut-off was set to 1200 eV and the one-dimensional 1st Brillouin zone was sampled with 34 k-points.As seen in Fig. 2, both the continuum bands and the bound state at −7 eV fit the DFT bands very well.In the DFT band structure, there are a few additional continuum bands with vanishing k z dispersion (thin lines in Fig. 2 a).As these states have zero group velocity along z, they do not contribute to the photocurrent and can be ignored for the present problem.
Figure 3 shows the wave functions obtained with the present matching method (red lines) for (a) the bound state ψ i at E i = −7.2eV and (b,c) the final continuum state ψ f at E f = 14.0 eV = E i + ℏω(He I).The wave ψ f is compared with the free wave (thin black line) and the exact solution (dashed blue line) of the (non-periodic) single well problem solved in the usual way by matching the plane wave eigenstates of exact kinetic energy in each region.The numerical wave function is almost identical with the exact wave.This shows that the present match- FIG. 3. Bound state (a, ψi at Ei = −7.2eV) and continuum state (b,c ψ f at E f = 14.0 eV) wave functions (red lines) obtained by numerical solution of the one-dimensional periodic potential model with PW basis set.ψ f is compared with the exact wave function of the non-periodic potential well (dashed blue lines) and to the free wave (thin black lines).The vertical dashed lines indicate the borders of the potential well.
ing method yields accurate photoemission wave functions with a moderate PW cut-off (∼ 1000 eV) comparable to ground state calculations.However, the number of k zpoints must be taken much larger than in ground state calculations (where a single k z is sufficient for well separated slabs) in order to obtain photoelectron waves of all possible perpendicular momenta.We also see from Fig. 3, that inside the potential well, the exact final state differs strongly from the free wave (thin black line), which highlights the shortcoming of the PW approximation.Volume 16 (2018) Ono, et al.

C. Photoemission intensity calculation
From the wave functions obtained with the matching method, we have calculated the photoemission intensity using Eq. ( 1) with A=e z .The intensity is plotted in Fig. 4 (red lines) as a function of final state energy.The curve shows some numerical noise below 20 eV, indicating that for some photoelectron energies, our simple scheme of selecting two band states in Eq. ( 2) is not very accurate.We are currently working on an improved method with a larger set of trial Bloch waves.For comparison, Fig. 4 also shows the photoemission intensity obtained in the (single) PW approximation (blue lines).At any given photon energy, the calculated intensities differ dramatically.In the important ultraviolet region below 40 eV, the energy dependence of the PW approximation is opposite to that of the matching method.This clearly shows that the PW approximation cannot be trusted for ARPES intensity calculations [5], especially at low energy.Interestingly, however, on a larger energy scale (Fig. 4, inset) the energy dependence looks similar, except for a shift of about 60 eV.This can be understood from the fact that the kinetic energy inside the potential well is 52 eV larger than outside.As the initial state is essentially localized inside the potential well, only the inner region matters for the photoemission calculation.Thus the correct intensity can roughly be reproduced using the PW approximation with a final state energy increased by about 50 eV.At present we do not know whether this is a particular features of the rectangular well potential, or whether such energy scaling also holds for other potential shapes.

IV. CONCLUSIONS
We have outlined a method for ARPES calculations from PW band structure codes and have successfully tested it with a simple one-dimensional model.The model parameters have been chosen to fit the k z band dispersion of graphene at the Γ-point obtained from a DFT calculation in the repeated slab mode.The DFT density of states above the vacuum level roughly follows the electron curve which shows that the PW basis is sufficiently complete.In the 1D model, photoemission final states are obtained by matching the proper combination of Bloch waves with the outer free wave.These numerical waves agree very well with exact solutions of the rectangular well potential.The photoemission intensity has a strong energy dependence and differs dramatically from the (single) PW approximation for any given energy.The final state band structure of the 1D model reproduces quite well the DFT bands at the chosen k || point (Γ), indicating that the present model calculation is realistic, and that the matching method can be extended to wave functions obtained from DFT supercell calculations in real 3D systems.

FIG. 1 .
FIG.1.Density of states (DOS) of graphene obtained from density functional theory (red solid line) compared with free electron DOS (dashed blue line).

FIG. 2 .
FIG. 2. Comparison of band structures.(a) kz band dispersion at the Γ-point of graphene from DFT. Thin lines indicate dispersionless bands.EF is the Fermi-level.(b) Bands obtained in the 1D model.Ei and E f indicate the levels used in Fig. 3. (c) Free electron bands.

FIG. 4 .
FIG. 4. Photoemission intensity from bound state at Ei = −7.2eV as a function of final state kinetic energy.Present wave function matching method (red) versus PW approximation (blue).