The complexity of photocathode designs and detector materials, and the need to model their performance for short pulse durations, the response to high-frequency photons, the presence of coatings and/or thinness of the absorptive layer, necessitates modifications to three-step and moments models of photoemission that are used in simulation codes. In this study, methods to include input from computationally intensive approaches, such as density functional theory to model optical properties and transfer matrix approaches to treat emission from the surface or transport past coatings, by means of parametric models are demonstrated. First, a technique to accurately represent optical behavior so as to model reflectivity and penetration depth is given. Second, modifications to bulk models arising from the usage of thin film architectures, and a means to rapidly calculate them, are provided. Third, a parameterization to model the impact of wells associated with coatings and surface layers on the transmission probably is given. In all cases, the methods are computationally efficient and designed to allow for including input from numerically intensive approaches that would otherwise be unavailable for simulations.

## I. INTRODUCTION

Models of photo-excitation and electron transport are needed to predictively estimate the quantum efficiency $QE$ and emittance $\epsilon n,rms$ of photoemitters with coatings or heterostructures that meet the needs of future x-ray free-electron lasers such as the *Linac Coherent Light Source II* (LCLS-II) and *Dynamic Mesoscale Material Science Capability* (DMMSC) xFELs.^{1–3} The same methods are argued to be useful to respond to the demands of other technologically important applications and developments, such as (i) beam optics codes, such as MICHELLE^{4–8} for the treatment of pulsed or density modulated beams, (ii) simulations treating emission through and from negative electron affinity (NEA) and/or roughened surfaces and semiconductors,^{9–18} (iii) simulations in which processes generate a range of photon energies in exotic bulk materials that leads to internal photo-excitation effects as in detectors, photovoltaics, and optoelectronic devices,^{19,20} (iv) photo-enhanced emission mechanisms for materials with coatings or heterostructures used for energy conversion,^{21,22} and (v) photoemission processes modified by multi-photon emission,^{23} short pulse effects,^{24,25} or rapid heating from high intensity lasers.^{26–28} Such examples share one or more features of light absorption over a range of energies, transport within bulk materials, and/or electron emission through or past surfaces with coatings or heterostructures in a manner that complicates phenomenological models often used to treat photoemission and generally referred to as Three-step Models (TSMs)^{29–34} or Simple Moments models (SMMs).^{35,36}

The present work reformulates a program began in a prior work^{37} to upgrade the SMM to account for increasingly complex physics associated with modern photocathode candidates; to enable a better account of the physics of absorption, transport, and emission; and to implement the theory in simulation methods compatible with beam optics code requirements and the utilization of the theory for characterization and analysis of sources. We first briefly recount the application of SMM to bulk materials such as metals (e.g., copper^{34}) and semiconductors (e.g., multialkali antimonoides with or without cesiated surfaces^{38} or graphene layers^{2,39} or perovskites^{20,40,41}), in Sec. II. It then undertakes specific modifications to extend the applicably of the underlying physics models and thereby address opportunities associated with recent developments. In particular, improvements to the SMM will account for materials physics, optical parameters, and scattering; extend the optical models to frequencies for which measured data are not available; evaluate the transmission probability $D(E)$ using exact methods for surface barriers more complex than the simple triangular and metallic (or “Schottky–Nordheim”^{42}) barriers from which analytical models^{43} can be developed; and accounts for modifications associated with thin photocathodes.^{44} Such modifications are needed for modeling photoemission when conditions are rapidly changing because small duration charge bunches^{37} are drawn from metals and semiconductors with coatings and layers.^{2,45–47} The changes fall into five categories:

parameterizations of the optical constants and laser penetration depth extended to high energy regimes discussed in Sec. III,

material physics based parameters affecting, e.g., density of states factors, effective masses, and scattering terms based on density functional theory (DFT) and other methods outlined in Sec. III C,

replacement of the (infinite extent) bulk semiconductor material by a layer (or a

*thin film*) of finite thickness presented in Sec. IV,more complex transmission probabilities, accounting for applied fields, surface coatings, resonances, and reflectionless well/barrier combinations that require exact evaluation described in Sec. V, and

meso-scopic factors caused by surface contamination by high work function materials such as carbon, different crystal faces, effects of temperature and thermal-field emission, and field enhancement due to surface roughness or geometric emitters

^{28,48}contribute and require additional considerations.

The first four are treated in the present study; the last is treated separately^{17,49} or part of on-going studies (e.g., temperature-dependence of scattering factors for novel materials and their impact on the transport models).

The extensions to SMM are constrained by two requirements: first, their implementation is guided by the needs of simulation and design, in particular, by the needs of beam optics codes and device modeling requirements and second, they remain open to incorporating input from experiments, measurements, and details provided by computationally intensive approaches such as density functional theory^{50} and Monte Carlo simulations^{11} using parametric methods.

The need for computational expediency in contrast to comprehensive theoretical models^{51–54} is made apparent by efforts to model the complications of surface roughness, work function and other emission site statistical variation, delayed emission, and space charge forces^{17,55} that are implicated in emittance growth, halo formation,^{56,57} and non-optimal launch times associated with short pulse length demands. The relentless demand for higher beam brightness necessitates simulations of beam dynamics where fluctuations and initial particle distribution affect the prediction of photoinjector performance.^{28} How the beam optics codes treat the emission process, therefore, has bearing. A 3D beam optics simulation using a modern particle-in-cell (PIC) code for photoemission in an RF photocathode can have upward of $3$–$30\xd7106$ macro-particles or more per beam bunch. Similar large scale macro-particle counts are required when modeling micro- to meso-scale rough surfaces features and/or surface material property variations in order to predict transverse energy spread. For very short pulses, where the transport time of a photoexcited electron to the surface is comparable to the duration of the laser pulse or when scattered electrons contribute significantly to an emission tail existing after the laser pulse duration is over,^{7,9} a history of intensity and field strength conditions (to say nothing of temperature) at earlier times and the repetitive calculation of emission at *each time step* during the beam pulse further complicate the nature of the emission models that can be employed to model the origin of fluctuations. A computationally expedient and simpler methodology that permits generalizations enabling minimizing computation time per emitted particle is needed. The reformulations developed in the present work are restricted by computational expediency in how they fulfill that need yet bring in the desired physics. The additional demands associated with internal photoemission effects at all frequencies associated with x-ray and $\gamma $-ray detectors add to the complexity in a different but no less complicated manner.^{58}

## II. SIMPLE MOMENTS MODEL

A Simple Moments Model (SMM) derives from the same framework responsible for the Fowler Nordheim and Richardson equations^{59,60} and is structurally analogous to Three-Step Models (TSMs) of photoemission,^{30,31,33,61,62} a phenomenological approach in which separate models for absorption of photons, transport of electrons through bulk material, and emission of electrons past a surface barrier are combined. The nomenclature “Moments” is used because both quantum efficiency and emittance are proportional to the moments $\u27e8kz\u27e9$ and $\u27e8k\u22a52\u27e9$ for $QE$ and $\epsilon n,rms$, respectively.^{36,63} In either case, it is the distribution of emitted electrons over which the moments are taken that require the extensions to the SMM developed herein, with which the momentum component (and to what power) is considered being of secondary importance. For that reason, even though the SMM methodology applies equally to $QE$ and $\epsilon n,rms$, for narrative simplicity, it is both more convenient and likely more familiar to focus on the former. In its zero temperature limit, SMM models quantum efficiency (QE) by

where the scattering factor $F\lambda (x,E)$ is given by^{64}

where $\u210f\omega $ is the photon energy, $R(\omega )$ is the reflectivity, *R*(*ω*) ¼ reflectivity, $D(E)$ is the density of states (assumed parabolic), $\mu $ is the chemical potential (or Fermi energy at $T=0$), $\Phi $ is the work function, $E~=E+\u210f\omega $, and $Eg,Ea$ are the bandgap and electron affinity. The height above the Fermi level for metals $\varphi =\Phi \u22124QF$ is the work function reduced by the Schottky barrier lowering factor $4QF$, where $Q\u2261q2/16\pi \epsilon 0=\alpha fs\u210fc/4$ and $\alpha fs\u22481/137.036$ is the fine structure constant.^{65} The factor $xm(E)=cos\u2061\theta m$ defines the “escape cone” or the minimum angle for which $Ez=Ecos2\u2061\theta $ exceeds the barrier height and is

for metals characterized by a work function $\Phi $ or semiconductors characterized by an electron affinity $Ea$. The scattering factor $f\lambda $ depends on the ratio between the laser penetration depth $\delta (\omega )$ [see Eq. (32)] and the mean free path between scattering events $\u2113=(\u210fk/m)\tau $, where $\u210fk/m$ and $\tau $ are the velocity and the associated relaxation time of the photoexcited electron, resulting in

where $p(\omega ,E)=\delta (\omega )/\u2113(E~)$ (the tilde indicating that it is the energy of the *photoexcited* electron that is used). It assumes that electrons that have undergone a scattering event no longer have sufficient energy to surmount the surface barrier as governed by the transmission probability $D(E~z)$, resulting in limits on the angular and energy integrations on which $QE$ depends. The integral in Eq. (3) is analytically given by

### A. Metals

The SMM often makes further approximations to obtain analytic forms of Eqs. (1) and (2). For metals, near threshold, $\u210f\omega $ is close to $\varphi $, and therefore, $f\lambda (cos\u2061\theta ,E)\u2248cos\u2061\theta /[1+p(\omega ,\mu )]$. Doing so results in

where $p=p(\omega ,\mu )$. Replacing the integral with its trapezoidal approximation gives to the leading order,^{59,66}

which recovers the Fowler–DuBridge approximation $QE\u221d(\u210f\omega \u2212\varphi )2$ to leading order. Equation (7) enjoys good correspondence with data: a comparison with the measured data of Dowell *et al.*^{34} using representative values of the various frequency-dependent factors for copper (Cu) gives the correspondence shown in Fig. 1, for which the value of $p=m\delta /\u210fk\tau $ used entails that $\tau \u22481.3$ fs assuming $\delta =12.9$ nm and $\u210fk=2m(\mu +\u210f\omega )$ (evaluated for $\lambda =266$ nm or $\u210f\omega =4.661$ eV): this value is lower than, but nevertheless comparable to, the relaxation times in copper for similar conditions.^{34,66} It does not, however, account for variations in $\delta (\omega )$, $R(\omega )$, or $\tau (\mu +\u210f\omega )$ that will cause changes as a function of frequency.

### B. Semiconductors

A similar analysis for semiconductor parameters can be performed using a triangular barrier model for the transmission probability,^{2} giving

where $s\u2032=[(\u210f\omega \u2212Eg)/Ea]1/2$, $s=[(\u210f\omega \u2212Eg\u2212Ea)/Ea]1/2$, and $Ea$ and $Eg$ are the electron affinity and bandgap, respectively, which recovers the parametric form suggested by Spicer^{29} and which compares favorably to measured values (Figs. 8 and 9 of Ref. 2).

For bulk photocathodes under long illumination conditions, low intensities that do not heat the photocathode, and photon energies comparable to the surface barrier, the simple forms of Eqs. (7) and (8) are well-suited for simulation and characterization. Complications such as heating due to high intensity lasers leading to additional emission mechanisms,^{26,27,35,67} ultrafast pulses^{23,24,68} and delayed emission,^{7,9,37} evolving surface conditions during operation or designed structures,^{2,69} large area emission from thin materials^{44} with irregular surface conditions for which coatings^{70} or layers are present,^{17,39,45,46} or illumination at short wavelengths^{40} for which optical and material parameters are unavailable or inadequate require that the SMM be reformulated and methods developed for it to utilize more varied sources of parameters and models. An auxiliary goal of the present work is to ensure that the reformulated methods do not demand a computationally prohibitive numerical overhead or introduce black-box components that would undermine their utility in beam optics codes, device simulations, and characterization efforts (e.g., optimization routines available in symbolic and computational packages for the determination of parameters).

## III. OPTICAL MODEL

A simple model of the dielectric constant $\epsilon (\omega )=K(\omega )\epsilon 0$ is derivable from a damped harmonic oscillator model applied to a free-electron gas and can account for features of the optical properties of metals and semiconductors with large carrier concentrations,^{71,72} but the same underlying Lorentzian forms naturally emerge from a quantum theory of free-carrier absorption^{73,74} applied to semiconductors for a wide range of photon energies. Rakić *et al.*^{75–77} have developed algorithms for modeling peaks in the imaginary part of the dielectric function using a small number of Lorentzians for metals and semiconductors parameterized by optimization methods. Here, a means to incorporate many Lorentzians is developed, allowing excellent matching with measured data over a large energy range and, importantly, providing a means to parameterize Density Functional Theory (DFT) results for which data do not exist or have not been extended to frequency regimes of interest. For comparisons to experimentally measured dielectric constants in Figs. 2–8, 11, and 12, data have been aggregated (from tabulated data or digitally extracted) from commonly used sources for metals and semiconductors^{77–83} and compared to online databases.^{84}

### A. Lorentz–Drude (LD) model

A dielectric function $K(\omega )$ can be written as^{78}

where $\chi f(\omega )$ and $\chi b(\omega )$ are intraband and interband electric susceptibilities. The intraband susceptibility is described parametrically by the free-electron Drude–Zener model with the oscillator strength $f0$ and the damping rate $\Gamma 0$,^{76}

whereas the interband susceptibility is described parametrically by the simple semi-quantum model resembling the Lorentz result for insulators^{77} given by

where

is the plasma frequency for $\rho o$ electrons per unit volume and $Nb$ is the number of interband transitions with the frequency $\omega j$, oscillator strength $fj$, and damping rate $\Gamma j$.

The total dielectric constant is then the sum of the free and bound components, or $K(\omega )=Kf(\omega )+Kb(\omega )$, both of which are complex. Having obtained the dielectric function, the optical parameters $n$ and $k$ are evaluated from its real and imaginary parts $K(\omega )=Kr(\omega )+iKi(\omega )$ by

In turn, the reflectivity $R(\omega )$ and penetration depth $\delta (\omega )$ are obtained from

[compare Eq. (32)]. For example, by adding two phenomenological Lorenztian terms, the prediction of reflectivity for copper (see Fig. 60 of Ref. 35) is improved. Systematic methods to evaluate and optimize the weighting factors $fj$, damping factors $\Gamma j$, and resonant frequencies $\omega j$ have been treated by Rakić *et al.*^{77} and Adachi *et al.*^{85} and can be applied to measured^{79–82,86} or calculated (e.g., by DFT) values for both metals and semiconductors. In contrast, the rapid and repetitive determination of $R(\omega )$ and $\delta (\omega )$, as demanded by predictive models of photoemission and simulation models used in beam optics codes, requires a reconsideration of how the factors $(fj,\Gamma j,\omega j)$ may be determined, with methods allowing for flexibility and speed taking precedence and (depending on application) extended to frequencies for which measured data may not be available. A computationally expedient and simpler methodology that permits generalizations is needed.

### B. Metal optical model

The method is first demonstrated for metals, beginning with copper as a representative and well understood test case. $\epsilon (\omega )$ is taken to be the sum of three classes of terms: a low frequency *Drude* component based on Eq. (10) and characterized by $fd$ and $\Gamma d$, a high-frequency *Lorentz* component based on Eq. (11) characterized by $fo$, $\Gamma o$, and $\omega o$ such that the latter two are in some sense “large”, and *resonant* components based on the $j\u2208(1,N)$ frequencies $\omega j$ and associated with $fj$ and $\Gamma j$. They are found for metals by matching the imaginary part $Ki(\omega )$ to Eqs. (10) and (11) because the Kramers–Kronig relations^{71} ensure that $Kr(\omega )$ can be determined from the same parameters. By introducing $x=\u210f\omega $ and function $y(x)$, the forms of which are fashioned after the behavior of $Ki(\omega )$, and by defining $\gamma =\Gamma /\u210f$,

Furthermore, let the measured (or calculated) data $Ki(\omega )$ be processed such that it is mapped onto a *regularly* spaced set of values $\u210f\omega i\u2261xi$ (where $\Delta \u2261xi+1\u2212xi$) using, e.g., cubic spline fitting: the “$i$” subscript will refer to that set of data, whereas the “$j$” subscript is reserved to correspond to the resonant levels determined by $\u210f\omega j=xj$. Such a mapping enables methods drawn from finite difference techniques^{87} to be used for approximating derivatives and finding local maxima. Let $yi\u2261y(xi;a)$. The factors $(fj,\Gamma j,\omega j)$ are then determined as follows. Observe in advance that the $fj$ so chosen will be replaced by a fourth step and also that the plasma energy $\u210f\omega p$ evident in Eq. (11) has been folded into its definition. The method of extracting Lorentzian parameters is adapted from techniques familiar to spectral analysis^{88} and is undertaken in four steps:

evaluation of Drude parameters,

evaluation of Lorentz ($\u210f\omega o,\gamma o,fo$) parameters,

evaluation of Lorentz ($\u210f\omega j,\gamma j,fj$) parameters, and

correction to the resonant $fj$ parameters.

These steps are described next.

#### 1. Drude

Set $a\u21920$ and $\gamma \u2192\gamma d$ in Eq. (15) with

for a suitably chosen small $xn$ ($n\u22652$), where $ln\u2061y$ is approximately linear in $ln\u2061x$. Although $fd$ so evaluated gives a good fit to $Ki(\omega )$ for a small $\omega $, the non-optimization of the algorithm and its reliance on potentially noisy data result in estimates of $Kr(\omega )$ being slightly off there. That small offset [which affects $R(\omega )$] can be mitigated by $fd\u2192hfd$ with $h$ close to unity *after* Step 4 is performed: for copper, $h\u22480.9$ is used in Figs. 5 and 6, whereas $h=1.28$ is used for Figs. 7 and 8.

#### 2. Lorentz

Set $a\u2192\u210f\omega o\u2261xo$ and $\Gamma \u2192\gamma o/\u210f$. Approximate $x/y(x)$ by a quadratic in $x2$, or $x/y\u2248A+Bx2+Cx4$ for a range of $xa\u2264xi\u2264xb$ with $xa$ and $xb$ chosen to span a region for which the polynomial approximation to $x/y(x)$ in $x2$ is reasonable. Letting $g(x)=1/(xy(x))$, then

In Fig. 2, $xa=30$ eV and $xb=70$ eV were chosen for the evaluation. While the resonant $\gamma j$ will be small ($\u223cO(1$)), the $\gamma o$ will be large ($\u223cO(10$)).

#### 3. Resonant

Finding the $xj$ first requires removing the Drude and Lorentz components: this is accomplished by subtracting the theory (red line) from the measured data line in Fig. 2. Let $Y0$ be the resulting set of data points, with the “0” indicating that only the Drude and Lorentz components have been subtracted, but no resonant term has as yet been removed. $Y0$ is shown in Fig. 4. The curve $Y0(x)$ shows a series of peaks taken to be near the $xj$, the selection of which is an esthetic choice (e.g., the highest of the peaks and/or a shape that appears to match a Lorentzian and/or progress toward a desired behavior): “esthetic” means that no algorithm is involved, but rather, the $xj$ is selected on the basis of how they affect the shape of $y(x)$ when the $xj$ resonant contribution is removed. This entails that the $xj$ is *not* selected in the order of magnitude, but rather in the order that they best remove peaks and features. This is acceptable for present purposes because although such a process may be automated (in principle), (i) it meets the immediate needs of developing a library of optical parameters that can be called as required by simulation codes, (ii) it is accurate to the extent needed for simulations of photocathode yield at a given wavelength, and (iii) simulations of detector performance sum over a (large) range of frequencies responsible for photoexcited electrons so that small variations as a function of frequency are smoothed. The by-hand approach is, therefore, preferable for present purposes as it focuses on important features and does not preclude future automation of the process. The locations of the $xj$ for copper are indicated by the yellow-filled circles in Fig. 2. Let $xn$ be the local extremum (peak or valley) of a set of five points, and let $y~n$ be the corresponding values of $Y0(x)$, shown in Fig. 3, such that there are five pairs of coordinates ${(xn\u22122,y~n\u22122)\cdots (xn+2,y~n+2)}$. Using the five-point finite difference approximation to the derivatives [Appendix (A1.3.2) of Ref. 59],

then the location of the extremum $xm$ [where $y\u2032(xm)=0$] is approximated by

The use of a 5-point scheme is advantageous when the shapes of the extrema are more hyperbolic than parabolic, but the scheme can occasionally introduce undesired effects for closely spaced extrema due to the overlapping of the resonant terms (as occurs when $\u210f\omega $ is small). Introduce the coefficients $A$ and $B$ defined by

Introduce an initial value of $s\u2261\u2212A/2Bxm2$, which is then iterated according to

with four iterations generally being sufficient. Then, it is found that

This approach is an iterative method to find $a=xj$ in Eq. (15) when the location of its maximum is known. The last step is to remove the resonant term from the $Y0(xi)$ data points, or

after which the process is repeated to remove the next desired local extremum at $xj+1$ and $y~$ drawn from $Yj+1$ until the resulting $YNb$ is of a desired flatness, where $Nb$ is the number of resonant terms required.

#### 4. Weighting correction

The values of $fj$ returned in Step 3 were temporary, only serving to remove the influence of a local maximum or minimum: each subsequent one corrupts the earlier $fj$ because of the finite width of $y(x;xj)$ as governed by $\gamma j$. The final step is, therefore, to renormalize the weighting factors $fj$. With $\gamma j$ and $xj$ determined, the $fj$ are found by inverting the matrix equation

where $Y$ and $f$ are vectors of the length $Nb$ (number of resonant terms) such that $[Y]i=Y0(xi)$, $[f]i=fi$, and $[M]ij=y(xi;xj)/y(xm;xj)$ is an $Nb\xd7Nb$ matrix, the form of which is determined from Eq. (15). Note that $xm$ is evaluated for each associated $xj$. The values of $(\omega j,\delta j,fj)$ as determined by these steps, along with the Drude and Lorentz parameters, are shown in Table I for copper.

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

D | 0 | 0.0888 | 79.710 | L | 45.826 | 210.34 | 2611.5 |

11 | 1.435 | 0.8263 | 0.5554 | 4 | 8.789 | 5.070 | 0.771 8 |

9 | 1.901 | 0.9819 | −2.370 | 8 | 11.96 | 1.602 | 0.112 4 |

1 | 2.695 | 1.223 | 4.158 | 5 | 14.41 | 6.247 | 0.659 3 |

3 | 3.616 | 1.178 | 1.486 | 7 | 17.44 | 4.418 | 0.095 79 |

10 | 4.768 | 2.017 | −0.4978 | 6 | 25.43 | 3.842 | 0.239 7 |

2 | 5.06 | 2.363 | 4.396 |

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

D | 0 | 0.0888 | 79.710 | L | 45.826 | 210.34 | 2611.5 |

11 | 1.435 | 0.8263 | 0.5554 | 4 | 8.789 | 5.070 | 0.771 8 |

9 | 1.901 | 0.9819 | −2.370 | 8 | 11.96 | 1.602 | 0.112 4 |

1 | 2.695 | 1.223 | 4.158 | 5 | 14.41 | 6.247 | 0.659 3 |

3 | 3.616 | 1.178 | 1.486 | 7 | 17.44 | 4.418 | 0.095 79 |

10 | 4.768 | 2.017 | −0.4978 | 6 | 25.43 | 3.842 | 0.239 7 |

2 | 5.06 | 2.363 | 4.396 |

The same exercise can be performed for other metals such as gold, for which the parameters are shown in Table II, the behavior of $K(\omega )$ being depicted in Fig. 7, and the $R(\omega )$ and $\delta (\omega )$ plotted in Fig. 8. The range in $\u210f\omega $ over which $K(\omega )$ values are measured is smaller than that for copper (90 eV for Cu compared to 30 eV for Au), making the determination of the $L$-parameters affected. No effort was made to mitigate the defect in order to show the impact of non-optimized parameters: it is seen that the evaluation remains reasonably robust for the estimation of $R(\omega )$ and $\delta (\omega )$. Optimization beyond the straightforward implementation of Steps 1–4 is deferred to a separate study. The absence of an optimization step results in some metals (such as lead) being more challenging to parameterize well but enables the steps to be implemented far more easily. Specifically, metals such as lead can be characterized as having higher parameter sensitivity, which is a numerical, rather than physical, characteristic of their dielectric response and thus require more refined parameter optimization.

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

D | 0 | 0.0262 | 70.564 | L | 25.691 | 42.307 | 801.98 |

17 | 1.909 | 1.044 | −0.6409 | 6 | 9.32 | 0.5973 | −0.154 6 |

16 | 2.392 | 0.6644 | −0.6887 | 5 | 10.28 | 6.402 | 0.874 5 |

12 | 2.783 | 0.5378 | 2.258 | 8 | 12.3 | 2.601 | −0.515 2 |

13 | 3.108 | 0.681 | 1.451 | 4 | 13.22 | 2.988 | 0.857 8 |

11 | 3.887 | 1.737 | 4.258 | 7 | 17.43 | 2.403 | −0.247 1 |

14 | 4.705 | 1.310 | 0.5908 | 3 | 18.45 | 1.449 | −0.215 1 |

15 | 5.856 | 2.995 | 0.894 | 1 | 20.37 | 13.42 | 1.084 |

9 | 7.329 | 4.783 | 0.2931 | 10 | 20.85 | 1.485 | 0.035 35 |

18 | 7.685 | 1.152 | 0.5327 | 2 | 24.77 | 2.820 | −0.515 3 |

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

D | 0 | 0.0262 | 70.564 | L | 25.691 | 42.307 | 801.98 |

17 | 1.909 | 1.044 | −0.6409 | 6 | 9.32 | 0.5973 | −0.154 6 |

16 | 2.392 | 0.6644 | −0.6887 | 5 | 10.28 | 6.402 | 0.874 5 |

12 | 2.783 | 0.5378 | 2.258 | 8 | 12.3 | 2.601 | −0.515 2 |

13 | 3.108 | 0.681 | 1.451 | 4 | 13.22 | 2.988 | 0.857 8 |

11 | 3.887 | 1.737 | 4.258 | 7 | 17.43 | 2.403 | −0.247 1 |

14 | 4.705 | 1.310 | 0.5908 | 3 | 18.45 | 1.449 | −0.215 1 |

15 | 5.856 | 2.995 | 0.894 | 1 | 20.37 | 13.42 | 1.084 |

9 | 7.329 | 4.783 | 0.2931 | 10 | 20.85 | 1.485 | 0.035 35 |

18 | 7.685 | 1.152 | 0.5327 | 2 | 24.77 | 2.820 | −0.515 3 |

### C. Semiconductor optical model

The Drude–Lorentz (DL) susceptibility model provides a suitable phenomenological parametric representation of experimental values for semiconductor dielectric functions.^{89} The application of Steps 1–4 can be extended to the treatment of semiconductors but subject to modifications. In insulators and semiconductors, where Van Hove singularities [non-smooth points in the density of states (DOS) where $|\u2207k(Ek)|$ vanishes, and thus, $Dj(\u210f\omega )$ is not a smooth function of its argument] are of interest, this model can achieve the more detailed representation of band edges in the optical spectrum by including critical points associated with these singularities in the joint density of states.^{90} For the interband transitions between two bands, $c$ and $v$, with energies $Ec(k)$ and $Ev(k)$, the susceptibility in Eq. (11) can be written as

where $\omega cv(k)=Ecv(k)/\u210f$, $Ecv(k)\u2261Ec(k)\u2212Ev(k)$, $fcv(k)=2|Pcv|2/(mEcv(k))$, and $Pcv$ is the matrix element of the momentum operator. Accordingly, close to the band edge, the imaginary and real parts of the interband dielectric function $\Delta K(\omega )=4\pi \chi b(\omega )$ with susceptibility from Eq.(22) are given, respectively, by^{85,91,92}

where $V$ is the volume and the $\delta $-function ensures energy conservation. The sum over the $\delta $-function in Eq. (23) can be turned to the joint density of states $D$ by

where $Sk$ is the surface defined by constant $Ecv(k)$ and the index $j$ labels critical points known as Van Hove singularities^{93} for which $\u2207kEcv=0$. In three dimensions, the energy near critical points can be expanded as

Four possible Van Hove singularities exist, classified according to the number of negative coefficients $Mj:j\u2208{0,1,2,3}$, where $M0$ and $M3$ represent a maximum and a minimum in the spectrum where all $\alpha i$ are positive (maximum) or negative (minimum), and $M1$ and $M2$ are saddle points where one or two of the $\alpha i$ are negative. Various combinations of critical points can be used to fit the dielectric function for a large number of direct bandgap semiconductors.^{85,94–98} Near the $M0$ critical point, in the range of energy where the expansion above is valid, the imaginary part of $K(\omega )$ can be written as

where $p=2$, $Eg=Ecv(0)$ is the bandgap, and $\Theta (x)$ is the Heaviside step function. For all 3D critical points $Mj$, a general form of the dielectric constant is^{99}

with $\eta $ vanishing. For the model leading to Eq. (26), the energy is assumed to be relatively close to the critical point such that typically $|Eg\u2212\u210f\omega |\u2243Eg$. In the case of a resonance structure possessing features that are of interest, accurate parametric representations exist.^{91,92} The present intent, however, is to construct a dielectric response function over a wide range of energies by means of interpolation between the $M0$ band edge, where $K\u223c\u210f\omega \u2212Eg$, and regions of a high frequency, where $\u210f\omega \u226bEg$ and $K\u223c\omega \u2212p$. In particular, capturing general trends of $K(\omega )$ is desirable, consistent with the inherent sensitivity of integral formulations of $QE$ implicit in, e.g., Eqs. (1) and (2).

The extension of the Drude–Lorentz model to semiconductors has additional complications due to the lower free-carrier density of electrons in comparison with metals. The validity of the (quasiclassical) Boltzmann transport equation (and hence the Drude–Zener theory^{100}) used in the modeling of solids in the determination of the dielectric function is altered by quantum mechanical effects, as can be shown using either the density matrix or the second order perturbation theory (both giving the same result): the relaxation time at high frequencies becomes frequency dependent, and this in turn alters the absorption coefficient. Briefly, the absorption coefficient’s relation to the optical conductivity, which is proportional to the relaxation time $\tau (\omega )$ [Eqs. (114) and (126) of Ref. 101], entail that the frequency dependence of $\tau (\omega )$ transfers to the absorption coefficient. A $1/\omega 3$ dependence describes the absorption coefficient instead of the $1/\omega 2$ predicted by the quasi-classical Boltzmann equation and is revealed by a quantum mechanical calculation of the absorption coefficient for the polar optical scattering mechanisms for, e.g., GaAs, InAs, InP, and CdTe.^{101} At higher carrier concentrations, impurity scattering increases in importance, and the characteristic dependence of the product of the absorption coefficient with the index of refraction varies from $\omega \u22124$ to $\omega \u22122$ and is dependent on the ratio of the photon energy to the initial electron energy.^{73,74,100,102} As a result, a higher power of $\omega $ can appear in the denominator of Eq. (25). Practically, it means that the value of $p$ in Eq. (25) takes into account the quantum effects discussed in Sec. III C, for which $p$ is expected to range from 3 to 4 instead of 2. Because measured data at large $\omega $ can be lacking, a demonstration of the dependence will include a description of the DFT methods that are used to estimate the behavior of the optical constants at higher energies than are available from experimental data so far considered.

### D. Density functional theory

DFT calculations were performed considering an ideal system for calculating optical properties of the alkali antimonides Cs$3$Sb and CsK$2$Sb, as shown in Fig. 9. For present purposes, the parameterization analysis is only conducted for Cs$3$Sb. Electronic structure calculations were done using the Vienna *Ab initio* Simulation Package (VASP),^{103–105} including core state effects via the VASP implementation^{106} of Projector Augmented-Wave (PAW) methods.^{107} We used the Local Density Approximation (LDA)^{108,109} to DFT^{110,111} for basic relaxation of structures, choosing a plane wave cutoff energy of 520 eV. Lattice constants of 0.91 nm were chosen for Cs$3$Sb and 0.86 nm for CsK$2$Sb, in the Fm$3_$m space group, to approximate experimental values. For all structures, a $k$-mesh of $8\xd78\xd78$ was used for non-local exchange, resulting in actual $k$-spacings of $0.3\xd70.3\xd70.3$ per 0.1 nm. The $k$-mesh was forced to be centered on the $\Gamma $ point.

The properties calculated determine the proper choice of density functional, as well as the particular system of interest. Calculations for graphene, as well as other metallic and dielectric systems, have been treated previously.^{112–115} For the present optical calculations, we chose a hybrid functional (HSE06)^{116} that combines a screened Hartree–Fock approach within DFT. The expression for the complex imaginary dielectric function was obtained by summing over conduction bands:^{117} in a formulation reflecting conventions of VASP calculations and at the risk of some confusion, it is commonly represented as

where for this equation only, $e$ is the electron charge, $uck$ and similar are periodic parts of Bloch wave functions, $\Omega $ is the cell volume, $\omega k$ is a multiplicity factor for each k-point, and transitions are made from occupied to unoccupied states within the first Brillouin zone.

Real and imaginary parts of dielectric functions are connected by Kramers–Kronig relations^{118} from which the optical constants $R(\omega )$ and $\delta (\omega )$ are calculated as per Eqs. (13) and (14). A comparison of measured vs calculated values of $Ki(\omega )$ for both Cs$3$Sb and CsK$2$Sb are shown in Fig. 9. The DFT data for $\u210f\omega $ are shifted, given the inability of semi-local density functional models to reproduce the absolute value of the bandgap: such shifts in the $x$-coordinate are not uncommon. It is assumed, however, that for the present study, DFT is sufficiently accurate for qualitatively estimating the functional character of $Ki(\omega )$, which can then be adjusted parametrically to fit a target functionally.

Although experimentally measured $Ki(\omega )$ data exist for the multi-alkali antimonides, the energy range is generally not as extensive as simulation desires. Extending the energy range of $Ki(\omega )$ using DFT, for the purposes of estimating its functional character at higher energies, and performing a parameter fit results in the family of curves shown in Fig. 10, where a conventional $Eg=1.6$ eV value is chosen, and where the $p$-curves are normalized so that they pass through an anchor point chosen to be 90% of $Ki(\omega )$ at $\u210f\omega =1.7478$ eV, although the actual power that is selected depends on how the Lorentz–Drude factors are chosen and used in Fig. 11. The closest fits are for $3<p<4$, compatible with the discussion regarding the Drude–Zener theory. The $p=3.5$ curve passes through the high $Ki(\omega )$ values well so that $p$ is chosen. Performing the $\u210f\omega j$ analysis analogous to Steps 3 and 4 using the parameters shown in Table III, but using Eq. (25) instead of the Drude model of Eq. (10), gives rise to the red “Theory” line and seen to match the data well. From the theory line, the $R(\omega )$ and $\delta (\omega )$ relations can be obtained, as in Fig. 12.

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

35 | 1.336 | 0.322 3 | 2.981 | 32 | 3.306 | 0.096 01 | 0.656 5 |

34 | 1.538 | 0.248 2 | 5.702 | 8 | 3.394 | 0.133 4 | 1.791 |

14 | 1.595 | 0.043 19 | 4.071 | 17 | 3.496 | 0.089 38 | 0.832 6 |

33 | 1.639 | 0.082 28 | −0.2193 | 9 | 3.825 | 0.221 | 1.488 |

36 | 1.655 | 0.758 4 | −1.765 | 25 | 3.975 | 0.088 26 | 0.725 5 |

15 | 1.797 | 0.003 73 | 0.1212 | 5 | 4.087 | 0.114 2 | 2.868 |

16 | 2.099 | 0.111 9 | 1.405 | 18 | 4.24 | 0.056 75 | 0.270 |

7 | 2.224 | 0.165 4 | 2.939 | 6 | 4.418 | 0.291 5 | 2.625 |

1 | 2.337 | 0.170 6 | 12.19 | 10 | 4.583 | 0.116 3 | 1.174 |

29 | 2.394 | 0.116 1 | −1.865 | 19 | 4.703 | 0.102 6 | 0.179 6 |

2 | 2.551 | 0.158 3 | 6.065 | 11 | 4.837 | 0.146 9 | 0.421 4 |

28 | 2.693 | 0.093 57 | −0.8654 | 21 | 4.953 | 0.113 4 | 0.099 34 |

3 | 2.792 | 0.128 4 | 7.723 | 20 | 5.097 | 0.123 5 | 0.158 5 |

27 | 2.839 | 0.171 3 | −1.645 | 12 | 5.441 | 0.091 88 | 0.582 8 |

30 | 2.894 | 0.077 55 | 0.2366 | 23 | 5.599 | 0.098 76 | −0.113 5 |

31 | 3.041 | 0.131 0 | 0.1525 | 24 | 5.709 | 0.122 5 | −0.077 27 |

4 | 3.112 | 0.336 4 | 6.963 | 22 | 5.883 | 0.123 3 | −0.217 9 |

26 | 3.232 | 0.123 1 | −1.847 | 13 | 6.01 | 0.185 9 | 0.775 4 |

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

35 | 1.336 | 0.322 3 | 2.981 | 32 | 3.306 | 0.096 01 | 0.656 5 |

34 | 1.538 | 0.248 2 | 5.702 | 8 | 3.394 | 0.133 4 | 1.791 |

14 | 1.595 | 0.043 19 | 4.071 | 17 | 3.496 | 0.089 38 | 0.832 6 |

33 | 1.639 | 0.082 28 | −0.2193 | 9 | 3.825 | 0.221 | 1.488 |

36 | 1.655 | 0.758 4 | −1.765 | 25 | 3.975 | 0.088 26 | 0.725 5 |

15 | 1.797 | 0.003 73 | 0.1212 | 5 | 4.087 | 0.114 2 | 2.868 |

16 | 2.099 | 0.111 9 | 1.405 | 18 | 4.24 | 0.056 75 | 0.270 |

7 | 2.224 | 0.165 4 | 2.939 | 6 | 4.418 | 0.291 5 | 2.625 |

1 | 2.337 | 0.170 6 | 12.19 | 10 | 4.583 | 0.116 3 | 1.174 |

29 | 2.394 | 0.116 1 | −1.865 | 19 | 4.703 | 0.102 6 | 0.179 6 |

2 | 2.551 | 0.158 3 | 6.065 | 11 | 4.837 | 0.146 9 | 0.421 4 |

28 | 2.693 | 0.093 57 | −0.8654 | 21 | 4.953 | 0.113 4 | 0.099 34 |

3 | 2.792 | 0.128 4 | 7.723 | 20 | 5.097 | 0.123 5 | 0.158 5 |

27 | 2.839 | 0.171 3 | −1.645 | 12 | 5.441 | 0.091 88 | 0.582 8 |

30 | 2.894 | 0.077 55 | 0.2366 | 23 | 5.599 | 0.098 76 | −0.113 5 |

31 | 3.041 | 0.131 0 | 0.1525 | 24 | 5.709 | 0.122 5 | −0.077 27 |

4 | 3.112 | 0.336 4 | 6.963 | 22 | 5.883 | 0.123 3 | −0.217 9 |

26 | 3.232 | 0.123 1 | −1.847 | 13 | 6.01 | 0.185 9 | 0.775 4 |

### E. Perovskites

Experimental data for the dielectric function of cesium lead halide perovskites is available, as is computational data for the cubic phase.^{41,119,120} Cesium lead halide perovskites arrange in lower symmetry phases at room temperature, with the experimental data most likely measured for these phases. In the absence of more conclusive experimental and computational studies, approximate fits analogous to the methods above can be constructed using the available information with regard to the bandgap $Eg$ and the high-frequency dielectric function $K\u221e=\epsilon \u221e/\epsilon 0$.

The bandgap in cesium lead halide perovskites can be tuned by replacing one halogen with another^{120} from the largest in CsPbCl$3$, $Eg\u22483.1$ eV, to a smaller in CsPbBr$3$, $Eg\u22482.4$ eV, to the smallest in CsPbI$3$, $Eg\u22481.7$ eV. These approximate values can be estimated using combinations of experimental^{121,122} and computational^{41,120} data.

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

38 | −0.9919 | 1.434 | 0.071 98 | 6 | 12.94 | 2.627 | 0.497 6 |

56 | 1.940 | 0.5407 | −0.102 2 | 19 | 13.7 | 0.9547 | −0.050 37 |

31 | 2.434 | 0.6001 | 0.509 3 | 18 | 14.54 | 0.7026 | 0.098 8 |

55 | 2.495 | 0.4193 | −0.240 5 | 50 | 14.97 | 0.4653 | 0.134 8 |

54 | 2.649 | 0.1511 | −0.209 7 | 52 | 15.25 | 0.2310 | 0.045 29 |

11 | 2.792 | 0.6273 | 2.601 | 20 | 15.68 | 1.3720 | 0.236 2 |

1 | 2.959 | 0.1687 | 2.922 | 53 | 16.15 | 0.6062 | 0.014 08 |

32 | 3.049 | 0.2493 | −0.821 8 | 51 | 16.54 | 0.4444 | 0.081 02 |

2 | 3.255 | 0.5134 | −1.508 | 7 | 17.02 | 0.7430 | 0.333 2 |

33 | 3.769 | 0.4166 | 0.350 2 | 4 | 17.69 | 0.4592 | 0.591 6 |

12 | 4.002 | 0.539 | 0.535 2 | 21 | 18.74 | 1.518 | 0.238 6 |

34 | 4.168 | 0.4703 | 0.425 9 | 22 | 19.14 | 0.6398 | 0.061 48 |

41 | 4.516 | 0.4247 | −0.220 3 | 28 | 19.94 | 0.5764 | 0.110 7 |

15 | 5.127 | 1.0100 | −1.117 | 10 | 20.98 | 1.570 | 0.279 4 |

35 | 5.300 | 0.5067 | 0.467 4 | 23 | 22.09 | 1.656 | 0.114 4 |

42 | 5.669 | 0.4418 | −0.271 8 | 25 | 22.90 | 1.394 | 0.127 1 |

36 | 6.449 | 0.3569 | 0.013 81 | 24 | 23.92 | 0.8345 | 0.146 1 |

37 | 6.568 | 0.4259 | 0.019 61 | 26 | 25.52 | 1.294 | 0.084 56 |

16 | 6.679 | 2.5940 | −1.198 | 27 | 26.64 | 0.7474 | 0.061 6 |

13 | 7.429 | 0.2069 | 1.064 | 44 | 27.82 | 0.5566 | 0.024 73 |

3 | 8.043 | 2.0560 | 3.225 | 47 | 28.28 | 0.8834 | 0.022 8 |

14 | 8.912 | 0.3308 | 0.392 5 | 48 | 31.20 | 0.9443 | 0.019 05 |

9 | 9.424 | 1.1630 | 1.277 | 46 | 32.12 | 0.8113 | 0.020 17 |

17 | 10.18 | 0.8626 | 0.610 3 | 29 | 33.17 | 0.6884 | 0.010 66 |

39 | 10.36 | 0.1339 | 0.059 74 | 43 | 34.37 | 0.8436 | 0.024 05 |

5 | 10.93 | 0.6245 | 1.087 | 45 | 35.43 | 0.5430 | 0.022 2 |

8 | 11.37 | 0.5368 | 0.325 1 | 49 | 39.89 | 1.589 | 0.014 61 |

40 | 11.75 | 0.1583 | 0.068 8 | 30 | 40.91 | 0.7115 | 0.023 62 |

j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. | j
. | $\u210f\omega j$ . | γ_{j}
. | f_{j}
. |
---|---|---|---|---|---|---|---|

38 | −0.9919 | 1.434 | 0.071 98 | 6 | 12.94 | 2.627 | 0.497 6 |

56 | 1.940 | 0.5407 | −0.102 2 | 19 | 13.7 | 0.9547 | −0.050 37 |

31 | 2.434 | 0.6001 | 0.509 3 | 18 | 14.54 | 0.7026 | 0.098 8 |

55 | 2.495 | 0.4193 | −0.240 5 | 50 | 14.97 | 0.4653 | 0.134 8 |

54 | 2.649 | 0.1511 | −0.209 7 | 52 | 15.25 | 0.2310 | 0.045 29 |

11 | 2.792 | 0.6273 | 2.601 | 20 | 15.68 | 1.3720 | 0.236 2 |

1 | 2.959 | 0.1687 | 2.922 | 53 | 16.15 | 0.6062 | 0.014 08 |

32 | 3.049 | 0.2493 | −0.821 8 | 51 | 16.54 | 0.4444 | 0.081 02 |

2 | 3.255 | 0.5134 | −1.508 | 7 | 17.02 | 0.7430 | 0.333 2 |

33 | 3.769 | 0.4166 | 0.350 2 | 4 | 17.69 | 0.4592 | 0.591 6 |

12 | 4.002 | 0.539 | 0.535 2 | 21 | 18.74 | 1.518 | 0.238 6 |

34 | 4.168 | 0.4703 | 0.425 9 | 22 | 19.14 | 0.6398 | 0.061 48 |

41 | 4.516 | 0.4247 | −0.220 3 | 28 | 19.94 | 0.5764 | 0.110 7 |

15 | 5.127 | 1.0100 | −1.117 | 10 | 20.98 | 1.570 | 0.279 4 |

35 | 5.300 | 0.5067 | 0.467 4 | 23 | 22.09 | 1.656 | 0.114 4 |

42 | 5.669 | 0.4418 | −0.271 8 | 25 | 22.90 | 1.394 | 0.127 1 |

36 | 6.449 | 0.3569 | 0.013 81 | 24 | 23.92 | 0.8345 | 0.146 1 |

37 | 6.568 | 0.4259 | 0.019 61 | 26 | 25.52 | 1.294 | 0.084 56 |

16 | 6.679 | 2.5940 | −1.198 | 27 | 26.64 | 0.7474 | 0.061 6 |

13 | 7.429 | 0.2069 | 1.064 | 44 | 27.82 | 0.5566 | 0.024 73 |

3 | 8.043 | 2.0560 | 3.225 | 47 | 28.28 | 0.8834 | 0.022 8 |

14 | 8.912 | 0.3308 | 0.392 5 | 48 | 31.20 | 0.9443 | 0.019 05 |

9 | 9.424 | 1.1630 | 1.277 | 46 | 32.12 | 0.8113 | 0.020 17 |

17 | 10.18 | 0.8626 | 0.610 3 | 29 | 33.17 | 0.6884 | 0.010 66 |

39 | 10.36 | 0.1339 | 0.059 74 | 43 | 34.37 | 0.8436 | 0.024 05 |

5 | 10.93 | 0.6245 | 1.087 | 45 | 35.43 | 0.5430 | 0.022 2 |

8 | 11.37 | 0.5368 | 0.325 1 | 49 | 39.89 | 1.589 | 0.014 61 |

40 | 11.75 | 0.1583 | 0.068 8 | 30 | 40.91 | 0.7115 | 0.023 62 |

Similarly, the high-frequency dielectric constant can be deduced from experimental measurements of the exciton binding energy, provided the effective masses are known,^{41} or directly calculated^{120} using the density functional perturbation theory.^{123} Again, the exciton binding energy is most likely measured in low symmetry phases: 20 meV (I), 40 meV (Br), and 70 meV (Cl),^{124} which is consistent with binding energies reported earlier.^{122} Estimates based on experiments are given in Ref. 41, where $K\u221e=4.5$ (Cl), 4.8 (Br), and 5.0 (I). For CsPbI$3$, however, the dielectric constant of 5.0 is taken from measurements on hybrid organic–inorganic perovskite MAPbI$3$, which forms a cubic structure similar to CsPbI$3$, with Cs replaced by the organic cation of methylammonium (CH$3$NH$3+$ or MA). MA is one of two widely used organic cations (the other is formamidinium NH$2$CH$=$NH$2+$ or FA) to stabilize the cubic phase. In general, the dielectric constant in the range of $4.5$–$5$ is consistent with the results of first-principles calculations for several halides being 4.1 (Cl), 5.0 (Br), and 6.3 (I).^{120}

A simple fit is based on approximation of the dielectric function near the $M0$ van Hove critical point at the edge of the bandgap, believed to be formed by the valence and conduction band at the high symmetry point $R$ of the Brillouin zone. At the $M0$ point, and introducing $s\u2261\u210f\omega /Eg$, the real and imaginary parts can be approximated as

where $\Theta (x)$ is the Heaviside step function and $a$ and $b$ are numerically determined by matching the LDR model with DFT simulations (or measured data), but which are understood as arising from the coefficients of Eq. (26). For the three materials shown in Fig. 13, $a=1$ for all, and $b=$ 12.4, 16.0, and 21.2 for CsPbCl$3$, CsPbBr$3$, and CsPbI$3$, respectively. For small $\u210f\omega $, $p=2$, but as with the Drude–Zener discussion (Section III C), its value shall be set by fits to the data.

For the purposes of demonstrating the methodology on perovskites, therefore, first consider the more tractable lead halide perovskite PbCI$3$N$2$H$5$, as shown in Fig. 14, to demonstrate performance, with the results shown in Figs. 15 and 16 using Table IV parameters. The method is then applied to cesium lead halide cases CsPbX$3$ for $X=$ (Br, Cs, I), in Fig. 17. In these applications, the method is clearly not optimized but is rapid and flexible and provides a reasonable behavior of the critical parameters $R(\omega )$ and $\delta (\omega )$, particularly for high photon energies $\u210f\omega $. Moreover, it allows for a substantial adjustment of parameters in an algorithmic manner and, therefore, may lead to an automated determination of the resonant and DL parameters, closer to the spirit of Rakić *et al*.,^{75,77} but for many more $\omega j$ components. For example, many of the more rapid variations in $Ki(\omega )$ can be mapped by artificially restricting the size of $\gamma j$ and using more resonant terms ($\omega j$) where a detailed structure is present. It is cautioned that attention to fine structure, however desirable as a matter of accuracy, is not necessary in simulations that span a broad range in the spectrum: variation within that range is of greater interest and is more consequential when a range of frequencies are present (as for detectors) or can be investigated for particular frequencies (as for photocathdoes). Therefore, the LDR method of extracting $K(\omega )$, and thereby $R(\omega )$ and $\delta (\omega )$, is significantly advantageous for the simulation of internal photoemission processes, characterization of photocathode materials, and utilization by beam codes. Further optimizations are not precluded and will be taken up elsewhere. Finally, DFT simulations of perovskites with a top cesium layer cause changes in the work function^{58} and in turn introduce changes in the LDR parameters that can presumably be accommodated: such a study will be reported separately, but an indication of how surface layers affect emission are considered in Sec. V.

## IV. THIN FILM MODEL

A thin film model differs from the bulk calculation: in the former, reflections occur at the substrate (back contact), and such reflections from a metallic substrate have been long suspected to be potentially important.^{125} Specifically, with proper engineering, an increase in the quantum efficiency is expected.^{44} Such “etalon” cathodes^{70} further decrease the response time of a photocathode by curtailing how far into a bulk material electrons are photoexcited and thereby limiting the transport time to the surface.^{7} The methodology to analyze them is first applied to the well-known bulk-vacuum problem and then extended to treat a thin film on a metallic substrate.

### A. Bulk-vacuum emission

From Maxwell’s equations for the monochromatic electromagnetic field of the frequency $\omega $, $c\u2207\u2192\xd7H\u2192+i\omega K(\omega )E\u2192=0$ and $c\u2207\u2192\xd7E\u2192\u2212i\omega \mu (\omega )H\u2192=0$,^{118} where in non-magnetic materials’ permeability is $\mu 0=1/\epsilon 0c2$, it follows that

In general, the dielectric function, $K(\omega )=1+\chi f+\chi b$, includes both intra- and interband susceptibilities. In conductors, however, the imaginary part of the intraband component dominates according to Eq. (10) as $\chi f=i\sigma /\omega $, where $\sigma =f0\omega p2/(4\pi \Gamma 0)$ is the Drude conductivity and Eq. (29) then turns into $[c2\u22072+\omega 2]E\u2192=\u2212i\sigma \omega E\u2192$. Assuming the incident laser is along the $x^$-direction, the electric field $y^\u22c5E\u2192=Ey$ is given by $Ey(x,t)=Eoei\kappa x\u2212i\omega t$, with the associated magnetic field $Hz=Hoexp(i\kappa x\u2212i\omega t)$ with $Ho=(c\kappa /\omega )Eo$, where $\kappa $ is the wave number, and then Eq. (29) for the homogeneous medium entails,

The complex index of refraction $n^\u2261n(\omega )+ik(\omega )$ is then

The imaginary part $\u2111(n^)=k(\omega )$ leads to dampening as light penetrates the bulk material. The penetration depth introduced in Eq. (14) is found from how intensity decays, which is proportional to $|E|2$, and therefore,

Using the convention that terms in vacuum are designated by a “0” subscript and those in bulk by a “1” subscript, then when the wave is incident on the surface, part of the wave is reflected or $Er=rE0exp(\u2212i\kappa 0x\u2212i\omega t)$ and part is transmitted or $Et=tE0exp(i\kappa 1x\u2212i\omega t)$. Because $E(x)=t1Eoei\kappa x$ in bulk, $\delta (\omega )$ is independent of $x$. Demanding continuity of $E$ and $\u2202xE$ across the surface results in equations readily handled by introducing a matrix $M$ and a coefficients vector $C$ defined by

for which the relations of continuity at the surface are compactly expressed as

with the boundary conditions $t0\u22611$ and $r1\u22610$. The solution of the matrix equation results in the commonly known relations

in terms of which the reflection coefficient $R0=|r0/t0|2$ is given by the well-known result,

where $n^=n\u2212ik$ and $n0=1$ have been used.

### B. Thin film on a metallic substrate

Now, let the photocathode material be deposited to a thickness $L$, which can be nanometers to micrometers in thickness, on the surface of a metal substrate. This is schematically shown in Fig. 18. The coefficients $Cn$ are altered by the introduction of another relation that must be satisfied, in addition to Eq. (34), of the form

but now with the boundary conditions $t0\u22611$ and $r2\u22610$. The resulting matrix equation,

with $t0=1$ and $r2=0$, may be solved for $r0$ and $t2$. Introducing the compact forms [recall Eq. (31)]

where $n^0=1$, then it is found

Alexander *et al*.^{44} observe that the simple vacuum-bulk solution of Eq. (35) is recouped in the limit that $k1L\u2192\u221e$, but an equivalent method is to take $n^2\u2192n^1$ for which $A\u21920$, and Eq. (35) is directly recovered. Continuing, the coefficients for the photocathode film are given by

The reflection coefficient at the film-substrate boundary is then

as may have been anticipated by Eq. (36). The electric field $E$ in the film is the top element of $M1(x)\u22c5C1$ or

for which the intensity $I\omega (x)\u221d|E(x)|2$ may be evaluated given the material parameters $n^2$, $n^1$, and $n^0=1$. The interference between the transmitted and reflected waves, however, means that a constant dampening factor similar to Eq. (32) may not always be possible: depending on $n1$ and $k1$, absorption may peak within the film rather than at the surface, with consequences for quantum efficiency $QE$ and emittance $\epsilon n,rms$ due to losses associated with photoexcited electrons transporting back to the surface.

Multi-alkali antimonide photocathode materials (e.g., Cs$3$Sb and CsK$2$Sb) are under active consideration to meet the needs of future x-ray Free-Electron Lasers (xFELs) because they (along with III–V semiconductors such as GaAs) are well-established high quantum efficiency photocathodes once the appropriate surface treatments are implemented (although at the cost of limited lifetime).^{2} Two examples are indicative for parameters drawn from the parameters for Cs$3$Sb-on-Cu for $\lambda =532$ nm and CsK$2$Sb-on-Ag for $\lambda =650$ nm, referred to as cases A and B, respectively (compare to, for example, Cs$2$Te, with $0.8\u2272n\u22721.8$ and $0.3\u2272k\u22720.7$ for $\lambda =254$ nm^{126}), compare to values shown in Table V for common photocathode materials. In both cases, let $L=\lambda /3$. The intensity is proportional to $|E1(x)|2$ with $E1(x)$ given by Eq. (44). Also shown is the $L\u2192\u221e$ (bulk) calculation, where the decay is governed by Eq. (32). Although the values of $n(\omega )$ in both cases are approximately $\lambda /L$, the higher value of $k(\omega )$ in Case A leads to a rapid decline as $x$ increases so that reflection off the substrate and the interference effects it entails in $E1(x)$ are mitigated. In Case B, constructive interference effects can occur so that the intensity has peaks further from the surface. As a result, the escape cone is narrowed, and therefore, the $QE$ and $\epsilon n,rms$ will be affected.

λ (nm)
. | 1064 . | 800 . | 650 . | 532 . | 405 . | 355 . | 266 . | 206 . |
---|---|---|---|---|---|---|---|---|

$\u210f\omega $ (eV) | 1.17 | 1.55 | 1.91 | 2.33 | 3.06 | 3.49 | 4.66 | 6.02 |

Ag n | 0.405 | 0.469 | 0.571 | 0.74 | 1.10 | 1.28 | 1.51 | 1.56 |

Ag k | 7.36 | 5.28 | 4.06 | 3.08 | 2.11 | 1.81 | 1.44 | 1.23 |

Au n | 0.118 | 0.144 | 0.216 | 0.406 | 1.34 | 1.51 | 1.37 | 1.16 |

Au k | 6.73 | 4.75 | 3.53 | 2.43 | 1.89 | 1.49 | 1.64 | 1.06 |

Cu n | 0.414 | 0.328 | 0.365 | 0.582 | 1.26 | 1.25 | 1.54 | 1.13 |

Cu k | 7.44 | 5.31 | 4.01 | 2.87 | 2.20 | 1.93 | 1.77 | 1.66 |

Pb n | 2.97 | 2.19 | 2.76 | 2.79 | 1.56 | 1.23 | 1.05 | 0.708 |

Pb k | 6.50 | 4.48 | 3.72 | 4.05 | 3.76 | 3.25 | 2.50 | 2.16 |

Cs_{3}Sb n | 3.61 | 3.73 | 3.81 | 2.96 | 1.98 | 1.37 | 1.01 | 0.378 |

Cs_{3}Sb k | 0.0434 | 0.512 | 1.26 | 2.12 | 2.24 | 1.98 | 1.70 | 1.16 |

CsK_{2}Sb n | 3.19 | 3.17 | 3.11 | 3.21 | 2.34 | 1.68 | 1.23 | 0.717 |

CsK_{2}Sb k | 0.123 | 0.313 | 0.564 | 0.992 | 2.06 | 1.84 | 1.55 | 1.47 |

GaAs n | … | 3.68 | 3.83 | 4.13 | 4.42 | 3.54 | 3.66 | 1.26 |

GaAs k | … | 0.09 | 0.18 | 0.34 | 2.07 | 2.02 | 3.34 | 2.47 |

λ (nm)
. | 1064 . | 800 . | 650 . | 532 . | 405 . | 355 . | 266 . | 206 . |
---|---|---|---|---|---|---|---|---|

$\u210f\omega $ (eV) | 1.17 | 1.55 | 1.91 | 2.33 | 3.06 | 3.49 | 4.66 | 6.02 |

Ag n | 0.405 | 0.469 | 0.571 | 0.74 | 1.10 | 1.28 | 1.51 | 1.56 |

Ag k | 7.36 | 5.28 | 4.06 | 3.08 | 2.11 | 1.81 | 1.44 | 1.23 |

Au n | 0.118 | 0.144 | 0.216 | 0.406 | 1.34 | 1.51 | 1.37 | 1.16 |

Au k | 6.73 | 4.75 | 3.53 | 2.43 | 1.89 | 1.49 | 1.64 | 1.06 |

Cu n | 0.414 | 0.328 | 0.365 | 0.582 | 1.26 | 1.25 | 1.54 | 1.13 |

Cu k | 7.44 | 5.31 | 4.01 | 2.87 | 2.20 | 1.93 | 1.77 | 1.66 |

Pb n | 2.97 | 2.19 | 2.76 | 2.79 | 1.56 | 1.23 | 1.05 | 0.708 |

Pb k | 6.50 | 4.48 | 3.72 | 4.05 | 3.76 | 3.25 | 2.50 | 2.16 |

Cs_{3}Sb n | 3.61 | 3.73 | 3.81 | 2.96 | 1.98 | 1.37 | 1.01 | 0.378 |

Cs_{3}Sb k | 0.0434 | 0.512 | 1.26 | 2.12 | 2.24 | 1.98 | 1.70 | 1.16 |

CsK_{2}Sb n | 3.19 | 3.17 | 3.11 | 3.21 | 2.34 | 1.68 | 1.23 | 0.717 |

CsK_{2}Sb k | 0.123 | 0.313 | 0.564 | 0.992 | 2.06 | 1.84 | 1.55 | 1.47 |

GaAs n | … | 3.68 | 3.83 | 4.13 | 4.42 | 3.54 | 3.66 | 1.26 |

GaAs k | … | 0.09 | 0.18 | 0.34 | 2.07 | 2.02 | 3.34 | 2.47 |

The emergence of the constructive interference effects as $n^1$ changes can be investigated. Figure 19 considers a hypothetical example for which $\u2111(n^p)=k$ stays fixed at $k=0.3$, and the index of refraction $\u211c(n^p)=n$ changes from 1.0 to 4.5 for the case of a bulk material ($L\u2192\u221e$) and a thin film ($L=\lambda /3$) with $\lambda =600$ nm and $n^s=0.5+5.0i$. When $k$ is increased, the oscillations are damped and the bulk behavior recovered; conversely, when it is smaller, interference effects are stronger and the oscillations, similar to those in Fig. 20, are correspondingly larger, particularly when $L/\lambda $ is close to an integer.

For actual materials, the dependence of both $n$ and $k$ on $\omega $ results in maps of $I\omega (x)$ of greater complexity. The example of a layer of Cs$3$Sb, of thickness $L=160$ nm, on a substrate of copper in Fig. 21 shows both $I\omega (x)$, which better shows the interference effects giving rise to peaks in the intensity, and $ln\u2061(I\omega (x))$ in Fig. 19, which better gives a sense of the decay of the intensity into the film. Practically, the actual behavior of $n^p(\omega )$ affects choices regarding the thickness of the film effect and the optimal choice of the substrate. Theoretically, the departure of $ln\u2061(I\omega (x))$ from linearity affects both TSM and moments models^{2,35,37} of the quantum efficiency $QE$, and the presence of peaks in $I\omega (x)$ within the film has implications for emittance $\epsilon norms$.^{44} Both the exploration of configurations and material choices and the simulation of photoemitter performance are, therefore, enhanced by having a computationally efficient model of the optical properties and an ability to predictively estimate them by DFT as in Sec. III E.

## V. EMISSION PROBABILITY MODEL

Insofar as a high quantum efficiency is a primary metric of photocathode utility for photoinjectors,^{69} a notable feature of the highest QE photocathodes is the incorporation of cesium in the stoichiometry and/or its presence as a surface layer^{127–130} to reduce the work function or eliminate it as for the III–V negative electron affinity, or NEA, photocathodes.^{131,132} For the latter cathodes, Spicer^{131} and Fisher *et al.*^{133} proposed a narrow, or “interfacial,” barrier at the surface, often represented as a tall triangular barrier but of sufficient thinness that tunneling electrons mostly pass through. For other metals and semiconductors for which a submonolayer covering of cesium acts to reduce the barrier height to emission, the manner by which the barrier is reduced, is a consequence of dipole effects described successfully by phenomenological theories.^{66,134,135} Our investigations of coatings on surfaces using DFT have suggested that the hypothetical triangular barrier may, in fact, be better treated as a potential well^{37,43} for which the qualitative features of a triangular barrier are mimicked. The barrier-well structure resulting from those simulations qualitatively resembled a model of photoemission from a metal with a partial alkali or alkali earth coating on the surface proposed by Fowler in his treatment of “selective” photoemission (the transmission of electrons from the conduction band of a metal first through a barrier then past a well with the barriers and wells being square, or constant potential, regions), which, although challenged, matched data.^{136} Here, a quantification of the behavior resulting from a coating layer is undertaken. The evaluation of the transmission probability $D[E(k)]$ demonstrates that the electrons that most contribute to $QE$ are those for which the normal energy $Ez=\u210f2kz2/2mn$ [in contrast to $E\u22a5$ that appears in the mean transverse energy ($MTE$) associated with emittance]. These electrons are sufficiently above the conduction band minimum that effective mass modifications do not arise; that is, the “effective” mass of the electron can be taken as the free-electron mass (*nearly free-electron* approximation). This approximation is made throughout. In addition, for the one-dimensional (1D) transmission problem, it is convenient to use $kz$ rather than $Ez$ and to suppress the $z$-subscript, as shall be done below.

### A. Potential model of a coating

For integer $\nu $ in Eq. (45), $D(k)\u21921$ such that *all* incident electrons for a particular $k$ are transmitted. We have shown previously that analogous behavior occurs for rectangular, triangular, and parabolic wells,^{37,43} wherein for an arc of $k$ values, reflectionless transmission occurs. For the purposes of constructing an analytic representation of $D(k)$ useful for the simulation of $QE$, however, the behavior of $D(k)$ with respect to a particular $k$ for the $sech2$ potential is advantageous, and therefore, it shall be used exclusively here. The justification is because macro-averaging of the potentials associated with coatings, as obtained from DFT, closely conforms to the $sech2$ potential.^{37}

For the non-integer $\nu $ in Eq. (45), then reflection can occur and can be determined using the transfer matrix approach (TMA) following methods described previously.^{43,59,139} For a potential mimicking a coating, $V(x)$ is given by

where $xo$ specifies the location of the monolayer coating and $an$ is ostensibly given by the Bohr radius for the $n$th shell of the coating atom or $an=n\u210f/\alpha fsmc$. For cesium, $n=6$; therefore, $a6=0.3175$ nm. Also, a representative potential is shown in Fig. 22.

As in Ref. 43, $D(k)$ can be plotted with increasing $\Delta V$ to reveal when reflectionless conditions occur: an example for cesium-like parameters is shown in Fig. 23. Clearly, specific values of $\Delta V$ are associated with $D(k)=1$ for all $k$. In between, reflection for small $k$ occurs. $\Delta V$ is correlated to the strength of the dipole associated with sub-monolayer coatings of cesium on surfaces in Gyftopoulous–Levine theory,^{59,66,135} and is, therefore, related to the degree of submonolayer coverage of cesium. Therefore, a representative selection of $\Delta V\u2264Ry/4$ needs to be parameterized. A contour representation of that region in shown in Fig. 24.

The region between the maximum reflection near $\Delta V=0.1Ry$ and the reflectionless condition near $\Delta V=0.164Ry$ shows conditions ranging from significant reflection for a small $k$ to near total transmission for all $k$. For a step-function potential barrier, $D(k)$ is known analytically^{59} to be given by $D(k)=4kk\u2032/(k+k\u2032)2$ with $k\u2032=k2\u2212ko2$ and $\u210f2ko2/2m=Vo$ equal to the height of a step-function barrier, corresponding to the asymptotic limit (as the field vanishes) of Fowler and Nordheim’s triangular barrier $DFN(k)$ [compare Eqs. (15) and (23) of Ref. 37]. Such limits have but one independent variable (the barrier height), whereas the *ad hoc* triangular barrier of Spicer and Fisher *et al.* are further characterized by a field. As previously shown,^{37} the triangular barrier model is a fair representation of the $sech2$ well in that the behavior of $D(k)$ in between the reflectionless (integer $\nu $) regions qualitatively resembles the step potential $D(k)$. Therefore, a rapid two-parameter model of $D(k)$ for regions where reflection is significant for the $sech2$ potential wells is sought so as to quantitatively account for the presence of coatings.

### B. Analytic two-parameter D(k) model

Let each line shown in the waterfall plot of $D(k)$ of Fig. 24 be indexed by $j$ with the depth of the well given by $\Delta Vj\u2261fjRy$ and $fj=jfmax/Nv$, for which Fig. 24 sets $fmax=0.25$ and $Nv=96$. The reflectionless conditions then occur for $j=21$ and $62$ (corresponding to $f21=0.0547,\nu =1$ and $f62=0.1615,\nu =2$). The $Dj(k)$ for $39\u2264j\u226462$ are shown in Fig. 25, approximately corresponding to the region where $1.5\u2264\nu \u22642$. In examining the reflections that occur in that region, a sampling of which are shown in Fig. 25, it is evident that a better approximation to $D(k)$ for the values of well depth that show reflection is reminiscent of the step potential. A two-parameter candidate is

where $ka$ and $r$ are to be determined (and not to be confused with the notation for the optical constants in Secs. III and IV). Define a new $kj$ by the relation $Dj(kj)=1/2$ whereby $ka2r=3kj2r$ and resulting in

The best least-squares values of the remaining parameter $rj$ are straightforward, resulting in the approximations shown by lines in Fig. 25.

Beam simulations^{4,7} demanding evaluations of $QE$ over sub-$\mu $m-scale areas and sub-ps time steps entail a large number of calls to the algorithms for generating $QE$ when modeling an electron bunch: the number scales as $Nx2Nt$, with $Nx=L/\Delta x$ and $Nt=T/\Delta t$, with $L$ and $T$ being the diameter of the emission area and the duration of the pulse, respectively. Device simulations for detectors, for which pair production and Compton scattering can result in energetic electrons that pass a bulk material/a contact barrier,^{40} have analogous demands. Consequently, even though an approximate form of $D(k)$ given by Eq. (48) has become available, the requirements to perform a least-squares fit to the critical $r$-parameter and extract $kj$ from a TMA evaluation, are substantial obstacles to its usage. It is necessary to either pre-load values for representative conditions and interpolate between them or reduce the time required for an evaluation. The latter approach motivates additionally approximating $r$ and $kj$ by simple functions. As shown in Fig. 26, the approximations,

perform adequately. Therefore, a means to rapidly modify $D(k)$ for varying submonolayer coatings is available, or, conversely, to rapidly determine parameters for a given configuration. More importantly here, however, is that an approximation to $D(k)$ that replaces the *ad hoc* triangular barrier model is now made available and accounts for the effect of a coating or a surface layer (like graphene^{37,45,46}) so as to predictively estimate $QE$ from surfaces that have been modified. Rapid methods to treat interface barriers are described separately.^{140}

## VI. CONCLUSION

The calculation of quantum efficiency and photoexcitation in complex multicomponent materials and structures requires numerically efficient algorithms for beam optics and device simulation codes to model the generation of electron beams from a photocathode, predict detector behavior, or characterize performance. Three-step models (TSMs) and simple moments models (SMMs) used to treat bulk materials provide an effective and widely used basis for simulation. However, these approaches require extensions to their underlying physics models to treat novel and coated materials with surface modifications, particularly if exploiting wave interference (an “etalon” thin film cathode) to enable fast responding structures with improved quantum efficiency is desired or if extending predictive capabilities to frequency regimes where optical or material properties are lacking is needed. In this contribution, modifications to the SMM model to allow for the effects of thin films were developed and demonstrated across a number of technologically relevant materials.

The predictive capabilities of DFT to treat the optical properties of materials, and of TMA to model emission past barrier structures, have been shown here to satisfy the physics modeling needs. Their computational cost (ill-afforded in device simulation and beam optics codes) has been mitigated by developing parametric models. For DFT, optical parameters such as the reflectivity $R(\omega )$ and penetration depth $\delta (\omega )$ of Eq. (14) are well-described by parameterizing a Drude–Lorentz fit with additional resonant terms for both metals and semiconductors. A method to identify and evaluate the required frequency, damping term, and weight, or $(\omega j,\gamma j,fj)$ parameters was presented, vetted on metals and semiconductors used in photocathodes, and projected to treat perovskites where the optical constants are extended to large frequencies using DFT methods. For TMA, the types of barrier modifications associated with coatings were parameterized by a Poschl–Teller well representation and their relation to the transmission probability associated with triangular barriers for non-reflectionless conditions. These techniques (in conjunction with methods treating interface transport and temperature effects) are being adapted for inclusion into characterization and simulation codes, the performance of which shall be reported in a separate contribution.

## ACKNOWLEDGMENTS

This work was supported in part by the Los Alamos National Laboratory (LANL) Directed Research and Development Funds (LDRD) and partly conducted at the Center for Nonlinear Studies (CNLS) and the Center for Integrated Nanotechnologies (CINT), U.S. Department of Energy (DOE). This research used resources provided by the LANL Institutional Computing (IC) Program. LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218NCA000001). This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. We thank Dr. Oksana Chubenko (Arizona State University) for discussions regarding the optical constants and their modeling.

## DATA AVAILABILITY

The data that supports the findings of this study are available within the article.

## REFERENCES

*Matter-Radiation Interactions in Extremes*or MaRIE.

*IEEE International Vacuum Electronics Conference/International Vacuum Electron Sources Conference (IVEC/IVESC) 2012*(Monterey, CA, 2012), pp. 341–342.

*IEEE International Vacuum Electronics Conference (IVEC), Monterey, California, 2014*(IEEE, 2014), pp. 79–80.

*IEEE International Vacuum Electronics Conference IVEC*(IEEE, London, UK, 2017).

*Advances in Imaging and Electron Physics*, edited by P. Hawkes (Elsevier, San Diego, CA, 2007), Vol. 149, pp. 1–338.

*Solid-State Photoemission and Related Methods: Theory and Experiment*, edited by W. Schattke and M. A. Van Hove (Wiley-VCH, Weinheim, 2002), p. 51.

*An Engineering Guide to Photoinjectors Photocathode Theory*, edited by T. Rao and D. H. Dowell (CreateSpace Independent Publishing Platform, 2016), p. 350.

*Solid-State Physics: An Introduction to Principles of Materials Science*, Advanced Texts in Physics (Springer, Berlin, 2003).

*Handbook of Optical Constants of Solids*, edited by E. D. Palik (Academic Press, Orlando, FL, 1985), pp. 169–188.

*American Institute of Physics Handbook*, 3rd ed., edited by D. E. Gray (McGraw-Hill, New York, 1972).

*CRC Handbook of Chemistry and Physics*, 88th ed., edited by D. R. Lide (CRC, 2007), p. 12.126.

*Numerical Solution of Partial Differential Equations: Finite Difference Methods*, Oxford Applied Mathematics and Computing Science Series (Clarendon Press/Oxford University Press, Oxford, UK/New York, 1985).

*Reviews of Infrared and Millimeter Waves*, edited by K. J. Button (Academic Press, New York, 1983), Vol. 8, pp. 127–171.

*Electron Emission and Adsorption Phenomena*, Cambridge Series of Physical Chemistry (Cambridge University Press, 1935).