Field emission is an important process with a variety of applications. Quantitative predictions of such electron emission need to include details of the internal potentials that shape the electronic wavefunctions (and hence the tunneling probability), predictive analysis of the work function barrier (ΦB), and knowledge of the electron distribution at the surface that constitutes the supply function. Here, these various factors were all collectively considered based on a combined Monte Carlo-density functional theory approach. Results were obtained for both the field-dependent cold electron emission current density as well as photoemission from a short laser pulse. The method also allows for calculations of field-dependent emittance. The technique is general and could be extended to include plasmon–polariton modes, different thicknesses of coatings, and role of surface adsorbates and defects.

Field emission hinges on the quantum-mechanical tunneling phenomenon, wherein a strong electric field applied at a surface extracts electrons.1–10 Soon after the landmark work by Fowler and Nordheim,1,2 tunneling-related α-particle decay11,12 was discovered, followed by tunneling-assisted transport in solids.13 Progress in this area continued, with the electronics industry pushing its efforts toward the nanoscale with innovative work relating to conductive polymers,14,15 molecular wires,16 and molecular electronic devices.17,18 The tunneling concept has also been applied to numerous other fields such as the study of bio-chemistry and reactions between molecular species.19,20 Applications of field emission include sources for flat and thin displays, x-ray tubes, charged particle generation for microwave devices, microscopy, vacuum electronics, lithography, electric propulsion, free-electron lasers, and medical therapies.21–30 A potential benefit of field-emission displays over conventional cathode tubes is its lower power consumption. For many applications, the ability of field emitters to provide large current densities with high pulse repetition rates and brightness are all desirable attributes. In this regard, combined photoemission and field emission enable high brightness and current densities in excess of 100 A/cm2, while allowing for much smaller active cathode areas.31 Large area emitter arrays with a small spread in energy and high power output is another aspect being pursued in this area.32–36 For many applications, metal cathodes have an advantage over their semiconductor-based counterparts, as they are more robust, can better tolerate poor vacuum conditions, and have a lower contamination risk. So here we choose to study electron emission from copper cathodes as a representative example.

Techniques for calculations of electron emission current density are well-established starting with the Fowler–Nordheim (FN) theory1,2 for field emission and the Richardson–Laue–Dushman (RLD) approach37 for the treatment of thermionic emission. Related asymptotic forms were developed by Murphy and Good38 and Herring and Nichols.39 More recently, a generalized thermal-field emission equation was derived,40–42 and studies based on full three-dimensional considerations were reported.43 Analytic expressions for the current density have typically been based on the free-electron model of metals with Fermi–Dirac statistics, the Wentzel–Kramers–Brillouin (WKB) approximation for the tunneling coefficient, and an integration over electron states using a linearization of the Gamow factor.44–47 This factor (G) is given as G = [(8m*)½/ℏ] s1s2 [V(s)-E]½ ds, where m* = electronic effective mass, s1 and s2 are the classical turning points at the potential barrier V(s), and E is the normal energy of electrons incident on the barrier. However, approximations have been used for the Gamow factor starting from an exact triangular barrier solution,44 or alternatively, shape factors have been introduced in the evaluation of G and approximations made for practical and easy use.8,48 It may be mentioned that in some reports based on the above approaches,49,50 the tunneling coefficient T has been expressed in terms of G as T(E) = 1/[1 + exp{G(E)}], which may have limited accuracy at larger energies.

Though these contributions represent significant advances, some aspects could benefit from a careful re-examination for a more comprehensive analysis. In this context, a few issues that merit analysis are discussed below.
  1. First, the tunneling probability for emission from the surface depends on the details of the electronic wavefunction that is influenced by the potential close to the emitting surface. This potential can vary based on the crystal orientation and/or the presence of defects and adsorbates.51–54 Here, rather than assuming a simplistic form of the internal potential such as a triangular barrier, a more accurate analysis is carried out based on the density functional theory (DFT). Similarly, the work function barrier (ΦB) plays a role55 and is influenced by the material, its crystal orientation, the presence of possible defects, etc. Toward this end, ΦB is carefully evaluated here, taking account of the crystal orientation, on the basis of the DFT approach.

  2. The electron emission current density (Jz) at the surface is related to the electron distribution function f(k||,kz), where k|| and kz denote the components of the electronic wavevector along the transverse (parallel to the emitting surface) and normal directions, respectively. This surface electron distribution constitutes the relevant supply function. The relation between Jz, k||, and kz for a three-dimensional (3D) structure without any quantum confinement is given by

J z = ( q ) ( 2 π 2 m ) k z f ( k , k z ) T C ( k z ) k d k d k z ,
(1)
where TC(kz) denotes the transmission coefficient, and f ( k , k z ) denotes the electron distribution as a function of the wavevectors k z and k . Most treatments of the current density, including the FN theory, have relied on the assumption of an equilibrium Fermi–Dirac distribution for the electron ensemble at the surface. However, this assumption will likely be invalid for two reasons: (a) The application of a strong external electric field would mean that a region (on the order of the Debye or Thomas–Fermi screening length) near the surface would be subject to some non-negligible electric field. Consequently, the electrons would be subject to a net drift, effectively resulting in a drifted-heated type of distribution rather than an equilibrium, thermalized one. (b) Furthermore, in the case of photoexcitation (e.g., laser-generated free electrons), the distribution would clearly deviate from the usual Fermi–Dirac form. In view of the above, details of the electron distribution are computed here based on a Monte Carlo transport model with inclusion of electron–electron and electron–phonon scattering events and an internal electric field. Details of carrier generation in the case of laser excitation and the subsequent momentum and energy exchanges within the free-carrier ensemble (consisting of both the cold equilibrium sea of electrons and the photogenerated population) are also appropriately included. As a result of the above, it is likely that the distribution would be skewed in k-space, requiring the use of the full wavevector-dependent distribution function, rather than applying the simpler energy-dependent f(E) approach for Jz calculations.
  1. The wavevector-dependent tunneling coefficient TC(kz) has often been computed on the basis of the WKB approximation, though corrections and adjustments exist in the literature.56 In moving beyond the free-electron approximation, Jensen57 used a triangular representation of the barrier with an adjustable parameter with solution of Schrödinger's equation in terms of Airy functions. In other approaches, field emission currents were calculated self-consistently within the DFT.58,59 A recent contribution by our group did adopt a coupled Schrödinger-DFT numerical analysis60 for numerical calculations of the electronic wavefunctions, with a jellium model representing the cathode interior,61 as the basis for obtaining TC(kz).

  2. Recently, it has been shown that hot electrons in graphene can mediate photoemission that requires optical power densities that are considerably lower than for metal tips.62,63 In a similar vein, the idea of surface photoemission enhancements from metals by resonance tunneling of photoexcited electrons through discrete quantum states through coatings of finite thicknesses was probed by Shuklin et al.64 and Zhou et al.65 Here, we use graphene coating on the metal cathode to assess potential photoemission current density enhancements based on our coupled MC-DFT approach. In particular, two graphene monolayers are used to obtain predictions of the emitted photocurrent.

Taking the above aspects into consideration, the electron current densities have been calculated in this contribution as a function of applied electric fields for both cold field-emission and photogenerated electron emission. In the context of the latter, it may be mentioned in passing that short laser pulses incident on cathodes can help produce coherent electron bunches which would have applications to free electron lasers,66 microscopy and spectroscopy,67,68 function as single-electron sources,69 or be useful for carrier envelope phase detection.70 Here for simplicity, smooth, flat copper surfaces have been assumed, though our group has simulated field emission from non-planar surfaces including nano-emitter arrays,71,72 while also accounting for dynamic electron screening and geometric field enhancements.73–75 Since real emitters have a variety of shapes,76–81 and given the possibility of thermo-mechanical instabilities or uneven surface degradation82 during device operation, a full three-dimensional (3D) analysis or numerical techniques83–87 may be needed as the next step. Such extensions also may be warranted given the well-known inadequacy of consistently predicting field enhancement factors that compare well with experimental data in nano-emitters.88 However, such treatment is beyond the present scope and would be probed elsewhere.

Calculations of electron current densities based on field emission are carried out here for copper coated with graphene. A Monte Carlo scheme for electron transport is used to comprehensively include details of the electron phase-space distribution function at the emitting surface with relevant collisions taken into account. Evaluations of the transmission coefficient through a composite structure consisting of the copper cathode with two monolayers of graphene rely on DFT to provide the internal potential profile V(z). This potential is then used to solve for the wavefunctions based on solutions to the time-independent Schrödinger Wave Equation (SWE) to finally yield the wavevector-dependent transmission coefficient. For completeness, we present calculations for both laser-assisted electron field emission and that from the sea of unexcited cold electrons.

The organization of this contribution is as follows: The simulation method is discussed in Sec. II, and the overall approach has been broken down into several parts. A Monte Carlo implementation is first discussed for electron transport within the metal toward the emitting surface in the presence of an external field. Our approach to the wavevector-dependent tunneling coefficient is then discussed based on a numerical solution of the SWE that utilizes the Numerov technique to yield the wavefunctions. Specifics of the analysis for obtaining the potential profile, a necessary input to the SWE, are also then outlined based on DFT. Section III reports on the results obtained and pertinent discussions on the predicted behaviors. Finally, the conclusions and future research directions are provided in Sec. IV.

The distribution function f(k||,kz) is important for calculations of the emission current density since it characterizes the supply function, and was evaluated based on Monte Carlo (MC) techniques for electron transport in solids.89–91 The MC allows for the inclusion of all relevant scattering processes (optical and acoustic phonons, electron–electron, electron–plasmon, etc.) in k-space and accounts for the momentum and energy exchanges for the photogenerated electrons driven by the external bias. Inclusion of scattering requires proper considerations of the phonon modes in the system. In general, for material with p atoms per primitive cell, one has 3p acoustic modes and 3p-3 optical modes. In the present case of copper, the primitive cell of the lattice has one atom, and hence, optical phonons were not considered. There have been several reports in the past on the treatment of electron scattering with acoustic phonons.92–95 Here, the expressions developed by Jensen et al.96,97 were used, with the energy-dependent electron–(acoustic) phonon scattering time τep[E(k)] given as
τ e p ( E ( k ) ) = π 3 ρ v s 2 4 m ( k B T D ) Ξ 2 k ( T D T ) 5 × ( W ( 5 , T D T ) ) 1 ,
(2)
where W represents the Bloch–Gruneisen function that was approximated to an analytic form for convenience,98, kB denotes the Boltzmann constant, vs denotes the sound velocity (taken to be 4.76 × 103 m/s), TD denotes the Debye temperature (set at 343 K for copper), ρ denotes the mass density (=8.95 × 103 kg/m3), and Ξ denotes the acoustic deformation potential (=6 eV). The values for these parameters were taken from Jensen.97 The electronic scattering process with acoustic phonons was taken to be isotropic, with the polar and azimuthal scattering angles θ and φ (relative to the velocity direction prior to collision), chosen based on randomly generated numbers ri and rj as cos(θ) = 1–2ri and φ = 2πrj. Electron–electron (e–e) scattering involves collisions between two electrons with initial wavevectors k1 and k2 (and corresponding kinetic energies E1 and E2) yielding final states k′1 and k′2. Though numerous reports on the e–e process exist in the literature,99–102 the following relation97 was chosen for the energy-dependent scattering time:
τ e e ( E ( k ) ) = 8 K s 2 α f s 2 π m c 2 ( E ( k ) k B T ) 2 × ( ( 1 + [ E ( k ) E fermi π k B T ] 2 ) × γ ( 2 k q 0 ) ) 1 ,
(3a)
where
γ ( x ) = x 3 4 ( tan 1 ( x ) + { x 1 + x 2 } tan 1 ( x 2 + x 2 ) 2 + x 2 ) .
(3b)
with αfs being the fine structure constant, c being the velocity of light, and Ks being the screening factor with a value of 1.684. Choice of the final state ( k 1 and k 2) is somewhat involved but has been discussed in the literature.102–104 Briefly, the procedure involves defining the vectors g (=k2 − k1) and g ( = k 2 k 1 ), of which g (and hence |g|) are completely known prior to a collision. The polar and azimuthal angles θ0 and φ0 associated with the vector g are then also determined. Next, new angles θ and φ upon scattering are generated as follows: θ = 1 – (2ri)/[1+|g|2 (1 − r)/kd2], where kd = ({ne e2}/{ɛ0kBT})1/2, and φ = 2πrj, where ri and rj denote two random numbers between 0 and 1. The components of g are then obtained based on the following relations:
g x = [ c x p cos ( φ 0 ) cos ( θ 0 ) c y p sin ( φ 0 ) + c z p cos ( φ 0 ) sin ( θ 0 ) | g | ,
(4a)
g y = [ c x p sin ( φ 0 ) cos ( θ 0 ) + c y p cos ( φ 0 ) + c z p sin ( φ 0 ) sin ( θ 0 ) | g | ,
(4b)
g z = [ c x p sin ( θ 0 ) + c z p cos ( θ 0 ) | g | ,
(4c)
with cxp = sin(θ) cos(φ), cyp = sin(θ) sin(φ), and czp = cos(θ). Finally, knowing g based on its components as given in Eq. (4), the requisite wavevectors k 2 and k 1 are obtained from the relations: k 1 = k 1 ( g g ) / 2 and k 2 = g + k 1.

The method first proposed by Lugli and Ferry103 for the Monte Carlo implementation of electron–electron scattering is used here. In this scheme (LF for short), when the free flight of a primary electron is interrupted by such an e–e event, a secondary electron is chosen at random for the two-body scattering process. The total momentum and energy are then conserved for such paired collisions. However, the problem with the method, as outlined by Mosko and Moskova,105 is that in this process the second partner electron is scattered earlier as compared to its originally designated time-of-flight. This aspect can then effectively lead to an over-counting of the e–e scattering rate. Furthermore, since the scattering rates increase monotonically with energy, for a low energy electron that might be chosen as the secondary, the differential between its effectively chosen scattering rate and the true energy-dependent destined collision frequency will increase. Rigorously, therefore, the LF scheme should be modified for better accuracy. However, doing so would be a huge step and would involve completely new Monte Carlo re-calculations, which are beyond the present scope. Potential implementation of a scheme of the type proposed by Mosko and Moskova105 could be attempted in subsequent treatments. Instead, here we chose to retain the current results based on the following reasoning: (i) The inclusion of MC analysis for obtaining the electron distribution at the emission site is a new aspect that has not been included previously to the best of our knowledge. Hence, it still offers an improvement over the conventional scheme of computing electron emission currents. (ii) In cases involving field emission without an external excitation source (e.g., optical illumination), most of the electrons likely to be emitted from the metal would all be from near the Fermi energy. Hence, the differential in energy and thus scattering rates might not be very large. (iii) The scenario for overestimating the e–e scattering rates is most likely in situations involving electrons excited by an external source such as a laser. However, due to the Beer–Lambert law and the exponential decay in excitation intensity with depth, many of the excited electrons would be created near the surface. Given the shorter distance for emission, the probability for ballistic transport and a quick exit from the metal increases. Those excited electrons having the right outward velocity component for an exit would also have a higher energy and velocity. Hence, the effect of a change in e–e scattering rate on the population of outgoing electrons and their k-space distribution, many of them created near the surface, may not be very strong. (iv) The e–e process represents single-particle interactions alone, and this treatment can be improved by including many-body effects through electron–plasmon scattering. However, one would have to carefully treat possible coupled plasmon–phonon modes. Alternatively, one could employ molecular dynamics (MD) treatments to overcome such issues of scattering over-estimation, which our group has done in the past106,107 for electron transport and other systems.82,108 However, the MD aspects are again beyond the present scope.

It must be mentioned that the Pauli exclusion principle was taken into account while treating electron scattering. The implementation was based on a numerical scheme first suggested by Lugli and Ferry.109 Considerations of Pauli exclusion (PE) become necessary when simulating electron transport at high densities, as is the case in the present calculations. This scheme has been used by our group routinely in MC simulation work,110,111 and so all of the details are not provided. Very briefly though, for implementing PE, the normalized distribution function fc at the discretized k-space corresponding to the chosen final state is compared to a random number r between 0 and 1. For r > fc, the Monte Carlo transition to the chosen final state was accepted only for r > fc or else was treated as a self-scattering event. A large occurrence of such PE-related scattering restrictions tends to reduce scattering and increase the average electron ensemble velocity in the presence of an external electric field.

For completeness, a quick point about electron–plasmon interactions is perhaps necessary. Details of the electron–plasmon scattering have been reported, for example, by Lugli and Ferry.112 The scattering rate S e p l basically has two terms corresponding to an absorption or emission of a plasmon. These two terms depend on the factors Nq and Nq +1, respectively, with Nq representing the Bose–Einstein distribution characterizing the plasmons. Given the large value of the plasmon frequency (ωpl = [ne2/(m*ɛ)]½), where n denotes the electron density in copper, m* denotes the electron effective mass, and ɛ denotes the permittivity; the rate of plasmon emission is much higher compared to that for plasmon absorption. However, the emission process can only be permissible for electrons with energies higher than ℏωpl, which is only possible with the use of high-energy lasers. Since such high incident photon energies were not used in our simulations, the electron–plasmon process was ignored. Phonon–plasmon coupling may be an issue that merits some attention and will be examined elsewhere.

For pulsed laser photoexcitation of the surface, the intensity I(t) was assumed to have the usual secant-hyperbolic shape113 as I(t)  = I0 sech2{(t − Tp)/τ}, with Tp representing pulse width, and τ related to the full width at half maxima (tFWHM) through the relation τ = tFWHM/[2 cosh−1{(2)½}]. Thus, either Tp and τ or Tp and tFWHM could be the user-specified parameters. The laser excitation would then lead to an evolution in the number of electrons, N(t), within the MC simulations. For the temporal laser profile chosen, the electron number was related to the pulse width Tp and τ as N(t)  = N * [tanh{(t – Tp/2)/τ} + tanh{Tp/(2τ)}]/[2 tanh{Tp/(2τ)}], with N representing the total user-specified number of photoexcited electrons used in any calculation run of the MC. The total number of electrons simulated was Ntotal, with (Ntotal − N) representing the cold electrons within the Fermi sea. In implementing electron–electron scattering for any particle, a random number (ri) was used to select the secondary interacting species, which could either be a particle from the photoexcited group or an electron from the original Fermi sea. Thus, for example, with ri < N/Ntotal, the secondary would be chosen to be a photogenerated entity.

In the present calculation, an approach similar to that reported previously60,61 was used. In this scheme, the planar emitting region was divided into three zones with the z-axis taken to be normal to the emitting surface. The first region (Zone 1) was furthest from the emitting surface and spanned z ≤ zmin. This zone was assumed to be a charge neutral semi-infinite jellium. Zone 2 covered zmin < z  zmax with a net internal potential V(z) present and was divided into 300 points. Depending on the situation considered, the zmax boundary was either at the copper–vacuum interface or at the end of the second graphene monolayer. The relevant potential V(z) in this region was obtained from DFT calculations based on a suitable averaging procedure over the lateral directions as discussed later. Non-periodic deviations in structure along the lateral direction were ignored, and consequently, inclusion of defects and/or surface irregularities was beyond the present scope. Finally, Zone 3 spanned zmax< z and encompassed the vacuum region over which an external electric field (F0) was assumed to be applied.

Determining the tunneling coefficient [TC(kz)] for electron entry into Zone 3 requires calculations of the appropriate wavefunction ψ(z). Here, the single-particle wavefunction ψ(r) was taken to be ψ(r) = ψ(z) exp(ik׀׀r׀׀), with k׀׀ denoting the wavevector parallel to the surface and ψ(z) obtained from a solution of the one-dimensional Schrödinger wave equation (SWE),
[ 2 ( 2 m ) ] d 2 ψ ( z ) d z 2 + V ( z ) ψ ( z ) = [ E 2 | k | 2 ( 2 m ) ] ψ ( z ) .
(5)
In the above expression, V(z) denotes the internal potential, E denotes the electron energy, with and m* being the usual reduced Planck's constant and effective electron mass, respectively. In Zone 1, ψ(z) was taken to have the following form: ψ(z)  = exp(ikzz) + R exp(−ikzz), with R denoting a reflection coefficient to be appropriately determined. In Zone 2, the Schrödinger equation was solved numerically using the sixth-order Numerov scheme60,114–116 with the appropriate potential V(z) obtained from DFT calculations. A total of 300 points were used. Finally, in Zone 3 with a constant external electric field F0 directed along the negative z-axis, the wavefunction had the form: ψ(z)  = T exp(ikzz) [Ai(−η0) − i Bi(−η0)], with T being an unknown of the problem representing the transmission coefficient, and Ai(−η0) and Bi(−η0) the modified Airy functions with η0= {(2qm E0)/ħ2}1/3 [(E – EF – W)/(qE0) + z]. Solution to the complete wavefunction in the three zones was obtained numerically by discretizing the entire region of interest (zmin ≤ z ≤ zmax) into a set of N finite points in Zone 2. The objective was to obtain values ψ(z2i) of the piecewise continuous wavefunction at the ith point (=z2i) in Zone 2. This led to a total of N + 2 unknowns upon factoring in the unknowns R and T. Matching the wavefunction and its spatial derivative at the zmin and zmax boundaries yielded four equations, and similarly matching ψ(z) at the N − 2 internal boundaries of Zone 2 produced a complete numerical solution for the wavefunction at all discrete points as well as values for R and T. Finally, the requisite tunneling coefficient TC(kz) was obtained as TC(kz) = [{(2qm*F0)/2}1/3/kz]|T|2. It may be pointed out that here we used the traditional matching of wavefunctions between regions of different effective mass [i.e., matching Ψ(z) and m−1 [δΨ(z)/δz]] as pointed out by BenDaniel and Duke.117 More recent work by Zhu and Kroemer118 and then Harrison,119 however, establishes that the traditional conditions for matching at boundaries with discontinuous electron masses are not entirely correct. Instead, matching both m−1/2 Ψ(z) and m−1/2 [δΨ(z)/δz] is more accurate and physically correct. Also, for the graphene, the treatment by Ariel and Natan120 was used for calculating the effective mass.

The internal potential V(z) energies, work function ΦB, and charge densities were all calculated from first principles. The Vienna Ab initio Simulation Package (VASP), a density-functional pseudopotential software tool121–123 that relies on the projector augmented-wave method,124,125 was used. The calculations were performed within the generalized gradient approximation (GGA) to the exchange-correlation potential as parameterized by PerdewBurke–Ernzerhof.126 VASP was used to simulate a perpendicular electric field outside of the slab in the vacuum region.127,128 The basic structure studied in this work consisted of 10-layer slabs of Cu atoms for the (111) orientation. A lattice constant of 3.63 Å was used for all structures and was calculated from the fcc Cu bulk structure. For other simulations, also carried out and reported here, two monolayers of graphene were added to the Cu surface. Details on the DFT application can be found in one of our previous publications.60 

Using VASP, the internal potentials (including the electrostatic, Hartree, and exchange-correlation terms) were obtained. From the three-dimensional potential, an average planar potential V ¯ ( z ) was extracted using the following relation:
V ¯ ( z ) = 1 S V ( r ) d x d y ,
(6)
where S is the area of the unit cell. Results obtained for the internal potential profile for copper one and two monolayers of graphene deposition are shown in Fig. 1, with the three zones also shown for clarity. The C–C bond length calculated from DFT was 1.48 A, while the distances between the two graphene layers and the Cu (111) surface were 3.86 and 3.35 A, respectively. For the results shown in Fig. 1, an external electric field of 5 × 109 V/m was assumed to exist in the vacuum region, leading to the downward slope in Zone 3. The z = 0 point was set at the boundary between Zone 1 and Zone 2, and thus, zmin = 0. A similar approach was followed with other external field values chosen as desired for the different cases studied. Values for zmax can be seen to be at about 26.7 and 29.2 Å for the single and double graphene structures, respectively.
FIG. 1.

Results obtained for the internal potential profile for copper coated with one and two monolayers of graphene. The various zones alluded to for the calculation are shown. An external electric field of 5 × 109 V/m was assumed to exist in the vacuum region (Zone 3). The C–C bond length calculated from DFT was 1.48 A, while the distances between the two graphene layers and the Cu (111) surface were 3.86 and 3.35 A, respectively. (a) Copper with a single graphene layer on top and (b) copper capped with two graphene monolayers.

FIG. 1.

Results obtained for the internal potential profile for copper coated with one and two monolayers of graphene. The various zones alluded to for the calculation are shown. An external electric field of 5 × 109 V/m was assumed to exist in the vacuum region (Zone 3). The C–C bond length calculated from DFT was 1.48 A, while the distances between the two graphene layers and the Cu (111) surface were 3.86 and 3.35 A, respectively. (a) Copper with a single graphene layer on top and (b) copper capped with two graphene monolayers.

Close modal

In terms of Monte Carlo simulations for electron transport in copper, the energy-dependent scattering rates had to be computed first. Results for the rates of the electron interactions with acoustic phonons via the deformation potential and binary electron–electron collisions obtained from Eqs. (2) and (3), respectively, are given in Fig. 2 as a function of the electron energy above the Fermi level. The electron–phonon scattering shows a slight monotonic increase with energy, going up to about ∼6.6 × 1013 Hz at 3 eV. In contrast, the electron–electron scattering increases more sharply and is predicted to reach about 6 × 1013 Hz at 3 eV and then go past the electron–phonon values. Though not shown, the rate 1/τee becomes larger than the 1/τep values for energies beyond 3.4 eV. This clearly indicates that laser excitation will produce rapid momentum relaxation in its photogeneration path.

FIG. 2.

Electron scattering rates as a function of kinetic energy in copper. Both the electron–phonon and binary electron–electron (e–e) interaction events were taken into account. While both show monotonic increases, the e–e rate is predicted to surpass the frequency of electronic scattering with acoustic phonons at about 3.4 eV.

FIG. 2.

Electron scattering rates as a function of kinetic energy in copper. Both the electron–phonon and binary electron–electron (e–e) interaction events were taken into account. While both show monotonic increases, the e–e rate is predicted to surpass the frequency of electronic scattering with acoustic phonons at about 3.4 eV.

Close modal

Using the MC scheme, the cold electron distribution at the emitting surface was computed as a function of the wavevectors k|| and kz. A total of 30 000 electrons were used for the MC simulations, with uniform initial spatial distributions. The starting energy distribution was chosen to be a thermalized Maxwellian, with the energy Ei of the ith electron chosen according to Ei = −(kbT) ln(ri), with ri being a random number between zero and unity, T being the lattice temperature, and kb being the Boltzmann constant. For an externally applied F0, the internal field was assumed to have an exponentially decaying profile into the copper with increasing distance from the copper surface, based on ThomasFermi screening. The more rigorous approach of solving the coupled Poisson-SWE system was ignored here but could be an improvement for the future. Results obtained from the Monte Carlo simulations for the electron distribution function at the copper surface are shown in Figs. 3(a)3(c) for different applied fields. Two aspects can be seen from the plots as expected. First, the distribution roughly follows an exponential type of decay along the k|| axis. There is no shift as the binary e–e scattering conserves momentum and the interaction with acoustic phonons is randomizing. However, a distinct shift in the distribution along the kz direction is evident due to the applied electric field. Figure 3(a) shows the distribution function f(k||,kz) at an external field of 108 V/m, while Figs. 3(b) and 3(c) are obtained at fields of 5 × 108 and 5 × 109 V/m, respectively. An increasing asymmetric shift along the kz wavevector is apparent with progressively higher applied fields. Since the overall electron emission current as given in Eq. (1) depends on the integrated result with influence from both the TC(kz) and f(k||,kz), such an asymmetric shift in f(k||,kz) along kz can be expected to affect Jz.

FIG. 3.

Electron distribution function f(k||,kz) at the copper surface obtained from Monte Carlo simulations in copper at different external fields. An increasing asymmetric shift along the kz wavevector is apparent with higher applied fields. (a) Distribution function f(k||,kz) at an external field of 108 V/m, (b) distribution function profile for a 5 × 108 V/m field, and (c) same k-space distribution obtained at a field of 5 × 109 V/m.

FIG. 3.

Electron distribution function f(k||,kz) at the copper surface obtained from Monte Carlo simulations in copper at different external fields. An increasing asymmetric shift along the kz wavevector is apparent with higher applied fields. (a) Distribution function f(k||,kz) at an external field of 108 V/m, (b) distribution function profile for a 5 × 108 V/m field, and (c) same k-space distribution obtained at a field of 5 × 109 V/m.

Close modal

Next, the transmission coefficients TC(kz) for electron tunneling through the graphene-covered copper cathode were calculated as a function of kz. These values were computed by first obtaining the internal potential V(z) as shown in Fig. 1. Use of V(z) in the SWE to obtain a full solution based on the Numerov technique for all three zones using the details already outlined then led to the final TC(kz) values. All calculations were carried out assuming pure copper with the (111) orientation. The results obtained for a single overlayer of graphene on copper as a function of kz is shown in Fig. 4(a), with the external electric field F0 as the parameter. The transmission coefficients (TCs) for copper with two successive monolayers of graphene are given in Fig. 4(b), again for different applied fields. The presence of an additional monolayer reduces the TCs somewhat, as might qualitatively be expected. Due to the small thickness of the graphene monolayer, its confined state was at a high energy and so did not appear to lead to resonant tunneling. Use of thicker graphene layers could be an interesting calculation for the future. In any case, as expected, the transmission coefficient is predicted to increase with energy, have higher values at stronger applied fields, and to be in the [0,1] range.

FIG. 4.

Electron transmission coefficients (TCs) obtained as a function of energy with applied field as the parameter. (a) Case of a single graphene monolayer on copper, and (b) the TCs at various applied fields for two successive monolayers coated on copper.

FIG. 4.

Electron transmission coefficients (TCs) obtained as a function of energy with applied field as the parameter. (a) Case of a single graphene monolayer on copper, and (b) the TCs at various applied fields for two successive monolayers coated on copper.

Close modal
Next, the TC values as well as the distribution function transmission f(k||,kz) were used to obtain the electron cold emission current density based on the double integrations inherent in Eq (1). As a simple comparison, results for an equilibrium Fermi electron distribution, and assuming parabolic energy–wavevector relationships [ℏ2k||2/(2 m) = E|| and ℏ2kz2/(2 m) = Ez], were also used for computations of the current density Jz. Under simpler conditions, Jz can be expressed as
J z = ( q m k T ) ( 4 π 2 3 ) 0 d E z T C ( E z ) Ln [ 1 + exp ( ( E f E z ) ( k T ) ) ] .
(7)

Results obtained for Jz taking account of both the wavevector-dependent distribution f(k||,kz) from MC calculations as shown previously in Fig. 3 and TC(kz) as given in Fig. 4. For comparison, values for Jz from Eq (7) were also obtained, and the field-dependent current density plots are shown in Fig. 5 for graphene-coated copper. For the MC-based approach, a finite set of distribution functions at various corresponding fields (as shown in Fig. 3) were first calculated and then applied. The values for Jz based on thermalized electrons obeying the equilibrium Fermi–Dirac statistics leads to much lower currents as compared to those predicted from the MC calculations. The presence of the field penetrating the copper shifts the distribution f(k||,kz) to non-zero average values along the kz wavevector with increasing field. That enhances the overall transmission coefficient and also imparts energy into the electron system. The net result is the higher current densities predicted over the simpler thermalized case. At high values of applied electric fields, the disparity reduces as the tunneling of lower energy electrons is strengthened. For comparison, data from the work reported by Toijala et al.129 have been added. In Fig. 5, the data values for a clean copper surface (diamond symbols) are seen to lie between the results for the emission current density obtained on the basis of an equilibrium Fermi–Dirac (FD) electron distribution, and those predicted from a Monte Carlo (MC) generated nonequilibrium distribution function. Thus, clearly, the usual equilibrium FD underestimates the emission current. Though the data by Toijala et al.129 are somewhat lower than the MC predictions, the values seem to converge toward the nonequilibrium distribution result with increasing electric field. The discrepancy at the low fields between the values of Toijala et al. and the calculations might arise from the slight overcounting of the e–e scattering rate as discussed in Sec. II A and reported previously by Mosko and Moskova.105 The presence of graphene coating in the present case might also lead to slight variations.

FIG. 5.

Calculated dependence of current density Jz on the external electric field. In one case, the distribution f(k||,kz) obtained from MC output was used, while the other curve was obtained by using a thermalized Fermi–Dirac electron distribution [Eq. (7)] at the emitting copper surface. For the MC approach, a finite set of distribution functions at various corresponding fields, as shown in Fig. 3, were first calculated and then applied. Results from the work by Toijala et al.129 have also been included for comparison.

FIG. 5.

Calculated dependence of current density Jz on the external electric field. In one case, the distribution f(k||,kz) obtained from MC output was used, while the other curve was obtained by using a thermalized Fermi–Dirac electron distribution [Eq. (7)] at the emitting copper surface. For the MC approach, a finite set of distribution functions at various corresponding fields, as shown in Fig. 3, were first calculated and then applied. Results from the work by Toijala et al.129 have also been included for comparison.

Close modal

Next, the field-dependent output current density for the composite copper with a graphene structure was evaluated for time-dependent photogeneration, as a function of the applied field. A single pulse with a full width at half maximum of 1.5 ps was assumed with a 1.16 eV energy corresponding to a Nd-YAG laser at 1064 nm wavelength. The usual secant hyperbolic laser pulse shape was assumed (i.e., I(t)  = I0 sech2{(t − tpeak)/tp}), with tpeak = 1 ps denoting the time at peak intensity and tp = tFWHM/{2 cosh1[(2)½]} related to the full width at half maxima (=tFWHM) that was taken to be 1.5 ps in simulations. The penetration depth parameter (α) was taken to be 6.129 63 × 107 m−1, which is appropriate for copper at this wavelength. The time-dependent emission of photogenerated electrons obtained from MC-based simulations is shown in Fig 6(a). A total of 30 000 electrons were simulated. As might be expected, the duration for all the photogenerated electrons to fully exit the electrode is longer than the pulse width. The photocurrent density computed from such a laser pulse excitation as a function of applied field is given in Fig. 6(b). For comparison, the predicted current density from field emission alone based on a thermalized electron population is also shown in the figure. Between the two cases, the Monte Carlo-generated photocurrent density is always higher than that from cold emission and is predicted to be roughly constant. This result is similar to the field-independent form predicted by the Fowler–Dubridge equation,130,131 as also discussed by Jensen,42 Zhou and Zhang,132 and others. Though not shown for brevity, the current did increase at higher incident laser energy. Finally to facilitate some rough comparisons to other reports, data presented by Martini et al.133 are also shown in Fig. 6(b). The data were obtained at electric fields that were relatively modest, and so only values at the far right end of the plot were available. Nonetheless, the available data are seen to match the predicted MC results for the photoexcitation case.

FIG. 6.

(a) Monte Carlo simulation results for the time-dependent emission of electrons from the copper material at an external field of 3 × 109 V/m. (b) Predicted results of the electron current density variation with applied electric field for field emission of a thermalized electron population and Monte Carlo results for photoexcitation by a 1.16 eV laser. For comparison, two data points (shown by blue squares) as reported by Martini et al.133 are also included.

FIG. 6.

(a) Monte Carlo simulation results for the time-dependent emission of electrons from the copper material at an external field of 3 × 109 V/m. (b) Predicted results of the electron current density variation with applied electric field for field emission of a thermalized electron population and Monte Carlo results for photoexcitation by a 1.16 eV laser. For comparison, two data points (shown by blue squares) as reported by Martini et al.133 are also included.

Close modal

Finally, results for the electron distribution f(k||,kz) made available from MC simulations at different fields made it possible to obtain quantitative measures of the root-mean-square emittance ɛx.134 The parameter is an important metric for assessing the divergence of a beam or, alternatively, the ability to focus it while minimizing lateral spread.135,136 For uniform and rotationally symmetric emission from a flat surface (as assumed in the present calculations), the normalized emittance ɛnx,rms depends on the spread in the transverse momentum137 and scales as ɛnx,rms ∼ ρ[⟨kx2⟩]½, where kx represents the transverse wavevector (∼ky due to rotational symmetry) and ρ denotes the effective radius of the flat emitter. It can be assumed that values of the transverse wavector in vacuum would be maintained based on momentum conservation at the barrier. Evaluation of the distribution function, as performed through MC simulations here, enables numerical assessments of the emittance factor by taking suitable ensemble averages as a function of applied bias, the operating temperature, and even in situations of external photoexcitation. Results from calculations at different external fields for both the field-emission and photogeneration scenarios were obtained for the relevant [⟨kx2⟩]½ term. The field-dependent behavior normalized to the lowest field chosen was plotted and is shown in Fig. 7. As might be expected, the kinetic energy of photogenerated electron increases due to laser excitation. Hence, the emittance is higher with photoexcitation though the output current density is also much larger as seen from Fig. 6(b). The emittance is predicted to increase with field. Even though the field is only applied in the longitudinal direction (z-axis) direction, the internal scattering mixes the momenta and helps build up the transverse kx component.

FIG. 7.

Ratio of normalized emittance obtained from Monte Carlo simulations as a function of applied field. The ratio was obtained upon dividing by the value at the lowest electric field. Both cold field emission and photoexcitation processes were considered.

FIG. 7.

Ratio of normalized emittance obtained from Monte Carlo simulations as a function of applied field. The ratio was obtained upon dividing by the value at the lowest electric field. Both cold field emission and photoexcitation processes were considered.

Close modal

For completeness, it must be mentioned that here the classical particle-based Monte Carlo simulation was chosen, without any inclusion of quantum approaches, other than the use of the Fermi golden rule of time-dependent perturbation theory for calculations of the electron scattering rates. In contrast to classical processes that can be viewed as comprising of elementary events associated with probabilities, the interplay of phases and amplitudes gives rise to interference effects in the quantum cases that cannot be described as a cumulative sum of probabilities. A convenient approach to a seamless inclusion of both the classical and quantum descriptions is offered by the transport model based on the Wigner function138 and implemented by a signed particle approach.139 The Wigner approach can be reduced to the Boltzmann transport model when the quantum rules for particle evolution are replaced by classical ones. Monte Carlo techniques incorporating such Wigner function-based quantum treatments have successfully demonstrated.140,141 Since here our focus is not on a unified treatment of scattering with inclusive quantum tunneling into the vacuum, we have not attempted to incorporate the Wigner-based method. This is certainly a task for future work and will be reported elsewhere. Furthermore, here we have chosen to use the density functional theory to obtain the potential barrier profile V(z) and the work function ϕB. So unlike the Monte Carlo–Wigner function methods that rely on an assumed tunneling barrier, we are able to include more accurate representations of V(z).

Field emission is an important process with applications toward the development of bright electron sources for high-resolution electron microscopes, for charged particle generation from microwave devices, lithography, electronic displays, and some medical therapeutics. Quantitative predictions need to include details of the internal potentials that shape the electronic wavefunctions (and hence the tunneling probability); details of the work function barrier (ΦB) at the boundary that is dependent on material quality, and values of the electron distribution function at the surface that constitutes the supply function. Here, these various factors were all collectively considered based on a combined MC-DFT approach. Results were obtained for the field-dependent cold electron emission current density Jz(F0). In addition, analysis is also carried out for electron field emission following photogeneration by an external laser pulse. In this context of photogeneration, short laser pulses incident on cathodes can help produce coherent electron bunches with applications to free electron lasers, microscopy, and as single-electron sources.

Though relevant results were presented for copper coated with graphene, the calculations should be viewed as a basic first step since a number of interesting and relevant effects could still be considered. For example, though plasmon scattering was justifiably neglected in the present calculation, polariton–plasmon coupled modes could play a role and need to be evaluated. This is especially germane to photoemission. Also, materials other than copper would require the additional inclusion of electron scattering with optical phonons via the polar and deformation interactions. Furthermore, though the role of defects was not studied here, the DFT-based method used could also be used to analyze systems containing vacancies and/or other defects. The only challenge would be to simulate large volumes in order to eliminate artificially high defect densities, making it computationally cumbersome for numerical calculations.

This work was supported in part by a grant from the Air Force Office of Scientific Research (No. FA9550-22-1-0499). One of us (R.P.J.) acknowledges useful discussions with K. Jensen (University of Maryland). We thank an unknown reviewer for some insightful comments.

The authors have no conflicts to disclose.

Y. M. Pokhrel: Data curation (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Y. Iqbal: Data curation (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). S. C. Shrestha: Data curation (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). M. Sanati: Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal). R. P. Joshi: Funding acquisition (equal); Methodology (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
R. H.
Fowler
and
L.
Nordheim
,
Proc. R. Soc. London A
119
,
173
(
1928
).
2.
R. H.
Fowler
,
Proc. R. Soc. London A
117
,
549
(
1928
).
3.
Modern Developments in Vacuum Electron Sources
(Topics in Applied Physics Book 135), edited by G. Gaertner, W. Knapp, and R. G. Forbes
(Springer, New York,
2020
).
4.
W. W.
Doland
and
W.
Dyke
,
Phys. Rev.
95
,
327
(
1954
).
5.
R. G.
Forbes
and
J.
Deane
,
Proc. R. Soc. A
467
,
2927
(
2011
).
7.
R. G.
Forbes
,
Surf. Interface Anal.
36
,
395
(
2004
).
8.
K. L.
Jensen
,
J. Appl. Phys.
126
,
065302
(
2019
).
9.
K. L.
Jensen
,
Introduction to the Physics of Electron Emission
(
John Wiley & Sons, Inc.
,
Hoboken
,
NJ
,
2017
).
10.
K. L.
Jensen
,
IEEE Trans. Plasma Sci.
46
,
1881
(
2018
).
11.
R. W.
Gurney
and
E. U.
Condon
,
Nature
122
,
439
(
1928
).
13.
C.
Zener
,
Proc. R. Soc. London A
145
,
523
(
1934
).
15.
G. R.
Hutchison
,
M. A.
Ratner
, and
T. J.
Marks
,
J. Phys. Chem. B
109
,
3126
(
2005
).
16.
E. A.
Weiss
,
M. R.
Wasielewski
, and
M. A.
Ratner
, “
Molecular wires: From design to properties
,” in
Topics in Current Chemistry
, edited by
L.
DeCola
(
Springer
,
Berlin
,
2005
), Vol. 257, pp. 103–133.
18.
N.
Renaud
,
M.
Hliwa
, and
C.
Joachim
, “
Unimolecular and supramolecular electronics II: Chemistry and physics meet at metal–molecule interfaces
,” in
Topics in Current Chemistry
, edited by
R. M.
Mezger
(
Springer
,
Berlin
,
2012
), Vol. 313, pp. 217–268.
19.
D. N.
Beratan
,
J. N.
Onuchic
,
J. R.
Winkler
, and
H. B.
Gray
,
Science
258
,
1740
(
1992
).
20.
J. R.
Winkler
and
H. B.
Gray
,
J. Am. Chem. Soc.
136
,
2930
(
2014
).
21.
X.
Calderon-Colon
,
H.
Geng
,
B.
Gao
,
L.
An
,
G.
Cao
, and
O.
Zhou
,
Nanotechnology
20
,
325707
(
2009
).
22.
C. A.
Spindt
,
C. E.
Holland
,
P. R.
Schwoebel
, and
I.
Brodie
,
J. Vac. Sci. Technol. B
16
,
758
(
1998
).
23.
G.
Cao
,
Y. Z.
Lee
,
R.
Peng
,
Z.
Liu
,
R.
Rajaram
,
X.
Calderon-Colon
,
L.
An
,
P.
Wang
,
T.
Phan
,
S.
Sultana
,
D. S.
Lalush
,
J. P.
Lu
, and
O.
Zhou
,
Phys. Med. Biol.
54
,
2323
(
2009
).
24.
R. J.
Barker
and
E.
Schamiloglu
,
High-Power Microwave Sources and Technologies
(
Wiley-IEEE Press
,
New York
,
2001
).
25.
J. H.
Booske
,
R. J.
Dobbs
,
C. D.
Joye
,
C. L.
Kory
,
G. R.
Neil
,
G. S.
Park
,
J.
Park
, and
R. J.
Temkin
,
IEEE Trans. Terahertz Sci. Technol.
1
,
54
(
2011
).
26.
C.
Brau
,
Free-Electron Lasers
(
Academic Press
,
Boston
,
MA
,
1990
).
27.
28.
J. I.
Goldstein
,
D. E.
Newbury
,
P.
Echlin
,
D. C.
Joy
,
C.
Fiori
, and
E.
Lifshin
, “
Scanning electron microscopy and X-ray microanalysis
,” in
A Text for Biologists, Materials Scientists and Geologists
(
Plenum Press
,
New York
,
1981
).
29.
D. M.
Goebel
and
I.
Katz
,
Fundamentals of Electric Propulsion: Ion and Hall Thrusters
(
John Wiley & Sons
,
New York
,
2008
).
30.
J. E.
Polk
,
M. J.
Sekerak
,
J. K.
Ziemer
,
J.
Schein
,
N.
Qi
, and
A.
Anders
,
IEEE Trans. Plasma Sci.
36
,
2167
(
2008
).
31.
K. L.
Jensen
,
P. G.
O’Shea
,
D. W.
Feldman
, and
J. L.
Shaw
,
J. Appl. Phys.
107
,
014903
(
2010
).
32.
L.
Nilsson
,
O.
Groning
,
C.
Emmenegger
,
O.
Kuettel
,
E.
Schaller
,
L.
Schlapback
,
H.
Kind
,
J. M.
Bonard
, and
K.
Kern
,
Appl. Phys. Lett.
76
,
2071
(
2000
).
33.
J. S.
Suh
,
K. S.
Jeong
,
J. S.
Lee
, and
I.
Han
,
Appl. Phys. Lett.
80
,
2392
(
2002
).
34.
Y.
Li
,
Y.
Sun
, and
J. T. W.
Yeow
,
Nanotechnology
26
,
242001
(
2015
).
35.
Z.
Zhang
,
G.
Meng
,
Q.
Wu
,
Z.
Hu
,
J.
Chen
,
Q.
Xu
, and
F.
Zhou
,
Sci. Rep.
4
,
4676
(
2014
).
36.
L.
Wang
,
Y.
Zhao
,
K.
Zheng
,
J.
She
,
S.
Deng
,
N.
Xu
, and
J.
Chen
,
Appl. Surf. Sci.
484
,
966
(
2019
).
37.
O. W.
Richardson
and
A. F. A.
Young
,
Proc. R. Soc. London A
107
,
377
(
1925
).
38.
E. L.
Murphy
and
R. H.
Good
,
Phys. Rev.
102
,
1464
(
1956
).
39.
C.
Herring
and
M. H.
Nichols
,
Rev. Mod. Phys.
21
,
185
(
1949
).
40.
K.
Jensen
,
J. Vac. Sci. Technol. B
21
,
1528
(
2003
).
41.
K.
Jensen
and
M.
Cahay
,
Appl. Phys. Lett.
88
,
154105
(
2006
).
42.
K. L.
Jensen
,
J. Appl. Phys.
102
,
024911
(
2007
).
43.
K. L.
Jensen
,
J. Appl. Phys.
107
,
014905
(
2009
).
44.
R. G.
Forbes
,
Appl. Phys. Lett.
89
,
113122
(
2006
).
45.
J.
Ludwick
,
M.
Cahay
,
N.
Hernandez
,
H.
Hall
,
J.
O’Mara
,
K. L.
Jensen
,
J. H. B.
Deane
,
R. G.
Forbes
, and
T. C.
Back
,
J. Appl. Phys.
130
,
144302
(
2021
).
46.
D.
Biswas
,
J. Appl. Phys.
131
,
154301
(
2022
);
D.
Biswas
and
R.
Ramachandran
,
J. Appl. Phys.
129
,
194303
(
2021
).
47.
R. G.
Forbes
and
J. H. B.
Deane
,
Proc. R. Soc. A
463
,
2907
(
2007
).
48.
K. L.
Jensen
,
J. Appl. Phys.
111
,
054916
(
2012
), see https://doi-org.lib-e2.lib.ttu.edu/10.1063/1.3692571
50.
R. G.
Forbes
,
J. Appl. Phys.
103
,
114911
(
2008
).
51.
L.
Diaz
,
A. A.
Karkash
,
S.
Alshahri
,
R. P.
Joshi
,
E.
Schamiloglu
, and
M.
Sanati
,
Sci. Rep.
13
,
8260
(
2023
).
52.
M.
Maille
,
N. C.
Dennis
,
Y. M.
Pokhrel
,
M.
Sanati
, and
R. P.
Joshi
,
Front. Mater.
11
,
1145425
(
2023
).
53.
M.
Brown
,
L.
Diaz
,
A.
Aslan
,
M.
Sanati
,
S.
Portillo
,
E.
Schamiloglu
, and
R. P.
Joshi
,
Sci. Rep.
12
,
15808
(
2022
).
54.
M.
Brown
,
M.
Sanati
, and
R. P.
Joshi
,
J. Appl. Phys.
131
,
103301
(
2022
).
55.
L.
Liu
,
R.
Liang
,
J.
Wang
, and
J.
Xu
,
AIP Adv.
6
,
015102
(
2016
).
56.
P.
Zhang
and
Y. Y.
Lau
,
Sci. Rep.
6
,
19894
(
2016
).
57.
K. L.
Jensen
,
J. Appl. Phys.
85
,
2667
(
1999
).
58.
N. D.
Lang
,
A.
Yacoby
, and
Y.
Imry
,
Phys. Rev. Lett.
63
,
1499
(
1989
).
59.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
);
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
60.
S. N.
Sami
,
L.
Diaz
,
M.
Sanati
, and
R. P.
Joshi
,
J. Appl. Phys.
128
,
223302
(
2020
).
61.
Y.
Gohda
,
Y.
Nakamura
,
K.
Watanabe
, and
S.
Watanabe
,
Phys. Rev. Lett.
85
,
1750
(
2000
).
62.
F.
Rezaeifar
,
R.
Ahsan
,
Q.
Lin
,
H. U.
Chae
, and
R.
Kapadia
,
Nat. Photonics
13
,
843
(
2019
).
63.
R.
Ahsan
,
M. A.
Sakib
,
H. U.
Chae
, and
R.
Kapadia
,
Phys. Rev. Appl.
13
,
024060
(
2020
).
64.
F. A.
Shuklin
,
I. V.
Smetanin
,
I. E.
Protsenko
,
J. B.
Khurgin
,
N. V.
Nikonorov
, and
A. V.
Uskov
,
Appl. Phys. Lett.
118
,
181104
(
2021
).
65.
Y.
Zhou
,
R.
Ahsan
,
H.
Uk Chae
,
R.
Kapadia
, and
P.
Zhang
,
Phys. Rev. Appl.
20
,
014043
(
2023
).
66.
P. G.
O’Shea
and
H. P.
Freund
,
Science
292
,
1853
(
2001
).
67.
Y.
Zhu
and
H.
Durr
,
Phys. Today
68
(4),
32
(
2015
).
68.
H.
Yanagisawa
,
M.
Hengsberger
,
D.
Leuenberger
,
M.
Klöckner
,
C.
Hafner
,
T.
Greber
, and
J.
Osterwalder
,
Phys. Rev. Lett.
107
,
087601
(
2011
).
69.
G.
Fève
,
A.
Mahé
,
J. M.
Berroir
,
T.
Kontos
,
B.
Plaçais
,
D. C.
Glattli
,
A.
Cavanna
,
B.
Etienne
, and
Y.
Jin
,
Science
316
,
1169
, (
2007
).
70.
B.
Piglosiewicz
,
S.
Schmidt
,
D. J.
Park
,
J.
Vogelsang
,
P.
Groß
,
C.
Manzoni
,
P.
Farinello
,
G.
Cerullo
, and
C.
Lienau
,
Nat. Photonics
8
,
37
(
2014
).
71.
D.
Guo
,
S. N.
Sami
,
L.
Diaz
,
S.
Sanati
, and
R. P.
Joshi
,
J. Vac. Sci. Technol. B
39
,
054201
(
2021
).
72.
D.
Guo
,
W.
Milestone
, and
R. P.
Joshi
,
J. Appl. Phys.
129
,
173301
(
2021
).
73.
J. R.
Harris
,
K. L.
Jensen
,
D. A.
Shiffler
, and
J. J.
Petillo
,
Appl. Phys. Lett.
106
,
201603
(
2015
).
74.
E. M.
Vinogradova
,
E. N.
Egorov
, and
D. S.
Televnyy
,
Vacuum
127
,
45
(
2016
).
75.
D.
Biswas
,
G.
Singh
, and
R.
Kumar
,
J. Appl. Phys.
120
,
124307
(
2016
).
76.
C. A.
Spindt
,
I.
Brodie
,
L.
Humphrey
, and
E. R.
Westerberg
,
J. Appl. Phys.
47
,
5248
(
1976
).
77.
R. V.
Latham
,
High-Voltage Vacuum Insulation: Basic Concepts and Technological Practice
(
Academic
,
London
,
1995
).
78.
H.
Andrews
,
K.
Nichols
,
D.
Kim
,
E. I.
Simakov
,
S.
Antipov
,
G.
Chen
,
M.
Conde
,
D.
Doran
,
G.
Ha
,
W.
Liu
,
J.
Power
,
J.
Shao
, and
E.
Wisniewski
,
IEEE Trans. Plasma Sci.
48
,
2671
(
2020
).
79.
T.
Cao
,
L.
Luo
,
Y.
Huang
,
B.
Ye
,
J.
She
,
S.
Deng
,
J.
Chen
, and
N.
Xu
,
Sci. Rep.
6
,
33983
(
2016
).
80.
V. I.
Kleshch
,
V. A.
Eremina
,
P.
Serbun
,
A. S.
Orekhov
,
D.
Lützenkirchen-Hecht
,
E. D.
Obraztsova
, and
A. N.
Obraztsov
,
Phys. Status Solidi B
255
,
1700268
(
2018
).
81.
W.
Lei
,
C.
Li
,
M. T.
Cole
,
K.
Qu
,
S.
Ding
,
Y.
Zhang
,
J. H.
Warner
,
X.
Zhang
,
B.
Wang
, and
W. I.
Milne
,
Carbon
56
,
255
(
2013
).
82.
Z.
Zhang
,
M.
Giesselmann
,
J.
Mankowski
,
J.
Dickens
,
A.
Neuber
, and
R. P.
Joshi
,
J. Phys. D
50
,
185202
(
2017
).
83.
Y. B.
Zhu
,
P.
Zhang
,
A.
Valfells
,
L. K.
Ang
, and
Y. Y.
Lau
,
Phys. Rev. Lett.
110
,
265007
(
2013
).
84.
Y. B.
Zhu
and
L. K.
Ang
,
Phys. Plasmas
22
,
052106
(
2015
).
85.
W. D.
Kesling
and
C. E.
Hunt
,
IEEE Trans. Electron Devices
42
,
340
(
1995
).
86.
J. J.
Petillo
,
E. M.
Nelson
,
J. F.
Deford
,
N. J.
Dionne
, and
B.
Levush
,
IEEE Trans. Electron Devices
52
,
742
(
2005
).
87.
K. L.
Jensen
,
D. A.
Shiffler
,
M.
Peckerar
,
J. R.
Harris
, and
J. J.
Petillo
,
J. Appl. Phys.
122
,
064501
(
2017
).
88.
D.
Biswas
and
R.
Ramachandran
,
J. Vac. Sci. Technol. B
37
,
021801
(
2019
).
89.
V. V. A.
Camargo
,
A. C. J.
Rossetto
,
D.
Vasileska
, and
G. I.
Wirth
,
J. Comput. Electron
19
,
668
(
2020
).
90.
S.
Bosi
and
C.
Jacoboni
,
J. Phys. C.: Solid-State Phys.
9
,
315
(
1976
).
91.
C.
Jacoboni
and
L.
Reggiani
,
Rev. Mod. Phys.
55
,
645
(
1983
).
92.
A. V.
Lugovskoy
and
I.
Bray
,
J. Phys. D
31
,
L78
(
1998
).
93.
H.
Ibach
and
H.
Luth
,
Solid-State Physics: An Introduction to Principles of Materials Science
(
Springer
,
Berlin
,
2003
).
94.
E.
Kartheuser
and
S.
Rodriguez
,
Phys. Rev. B
33
,
772
(
1986
).
95.
J. M.
Ziman
,
Electrons and Phonons
(
Clarendon Press
,
Oxford
,
1962
).
96.
K. L.
Jensen
,
E. J.
Montgomery
,
D. W.
Feldman
,
P. G.
O’Shea
,
J. R.
Harris
,
J. W.
Lewellen
, and
N.
Moody
,
J. Appl. Phys.
110
,
034504
(
2011
).
97.
K. L.
Jensen
,
N. A.
Moody
,
D. W.
Feldman
,
E. J.
Montgomery
, and
P. G.
O’Shea
,
J. Appl. Phys.
102
,
074902
(
2007
).
98.
M.
Deutsch
,
J. Phys. A: Math. Gen.
20
,
L811
(
1987
).
99.
K.
Motizuki
and
M.
Sparks
,
J. Phys. Soc. Jpn.
19
,
486
(
1954
).
100.
J. J.
Quinn
,
Appl. Phys. Lett.
9
,
167
(
1963
).
101.
V. V.
Kabanov
and
A. S.
Alexandrov
,
Phys. Rev. B
78
,
174514
(
2008
).
102.
W. E.
Lawrence
and
J. W.
Wilkins
,
Phys. Rev. B
7
,
2317
(
1973
).
103.
P.
Lugli
and
D. K.
Ferry
,
Physica B
117 & 118
,
251
(
1983
).
104.
M. A.
Osman
and
D. K.
Ferry
,
Phys. Rev. B
36
,
6018
(
1987
).
105.
M.
Mosko
and
A.
Moskova
,
Phys. Rev. B
44
,
10794
(
1991
).
106.
R. P.
Joshi
and
R. F.
Wood
,
J. Appl. Phys.
84
,
3197
(
1998
).
107.
R. P.
Joshi
and
D. K.
Ferry
,
Phys. Rev. B
43
,
9734
(
1991
).
108.
S. N.
Sami
,
R.
Islam
,
R.
Khare
, and
R. P.
Joshi
,
J. Appl. Phys.
129
,
213303
(
2021
).
109.
P.
Lugli
and
D. K.
Ferry
,
IEEE Trans. Electr. Dev.
32
,
2431
(
1985
).
110.
W.
Milestone
,
S.
Nikishin
, and
R. P.
Joshi
,
MDPI Electron.
11
,
2997
(
2022
).
111.
W.
Milestone
,
D.
Guo
,
M.
Sanati
,
K. M.
Dowling
,
S.
Hue-Riege
,
L. F.
Voss
,
A.
Conway
, and
R. P.
Joshi
,
J. Appl. Phys.
129
,
195703
(
2021
).
112.
P.
Lugli
and
D. K.
Ferry
,
IEEE Electron. Dev. Lett.
6
,
25
(
1985
).
113.
J. T.
Manassah
,
M. A.
Mustafa
,
R. R.
Alfano
, and
P. P.
Ho
,
IEEE J. Quantum Electron.
QE-22
,
197
(
1986
).
115.
116.
R. H.
Landau
,
M. J.
Paez
, and
C. C.
Bordeianu
,
A Survey of Computational Physics
(
Princeton University Press
,
Princeton
,
NJ
,
2008
).
117.
D. J.
BenDaniel
and
C. B.
Duke
,
Phys. Rev
152
,
683
(
1966
).
118.
Q. G.
Zhu
and
H.
Kroemer
,
Phys. Rev. B
27
,
3519
(
1983
).
119.
W. A.
Harrison
,
J. Appl. Phys.
110
,
113715
(
2011
).
120.
V.
Ariel
and
A.
Natan
, “
Electron effective mass in graphene
,” in
2013 International Conference on Electromagnetics in Advanced Applications (ICEAA)
,
Turin, Italy
(IEEE,
2013
), pp.
696
698
.
121.
See https://www.vasp.at/ for more information about the VASP simulator.
122.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
123.
G.
Kresse
and
J.
Furthmuller
,
Phys. Rev. B
54
,
11169
(
1996
).
124.
125.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
126.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
127.
J.
Neugebauer
and
M.
Scheffler
,
Phys. Rev. B
46
,
16067
(
1992
).
128.
G.
Makov
and
M. C.
Payne
,
Phys. Rev. B
51
,
4014
(
1995
).
129.
H.
Toijala
,
K.
Eimre
,
A.
Kyritsakis
,
V.
Zadin
, and
F.
Djurabekova
,
Phys. Rev. B
100
,
165421
(
2019
).
130.
131.
132.
Y.
Zhou
and
P.
Zhang
,
J. Appl. Phys.
127
,
164903
(
2020
).
133.
I.
Martini
,
E.
Chevallay
,
V.
Fedosseev
,
C.
Hessler
, and
M.
Martyanov
, “
Photoemission studies of copper
,” in
1st Topical Workshop on Laser Based Particle Source
(CERN,
Geneva
, Switzerland,
2013
).
134.
M.
Reiser
,
Theory and Design of Charged Particle Beams, Wiley Series in Beam Physics and Accelerator Technology
(
Wiley
,
New York
,
1994
), p.
xix
.
135.
C.
Brau
, in
Physics and Applications of High Brightness Electron Beams: Proceedings of the ICFA Workshop Chia Laguna
,
Sardinia, Italy
,
2002
, edited by
J.
Rosenzweig
,
G.
Travish
, and
L.
Serafini
(
World Scientific
,
Singapore
,
2004
), p.
20
.
136.
137.
K. L.
Jensen
,
D. A.
Shiffler
,
J. J.
Petiloo
,
Z.
Pan
, and
J. W.
Luginsland
,
Phys. Rev. Special Topics Accelerators Beams
17
,
043402
(
2014
).
139.
M.
Nedjalkov
,
D.
Querlioz
,
P.
Dollfus
, and
H.
Kosina
,
Wigner Function Approach
(
Springer
,
New York
,
2011
), p.
289
.
140.
M.
Nedjalkov
,
H.
Kosina
,
S.
Selberherr
,
C.
Ringhofer
, and
D. K.
Ferry
,
Phys. Rev. B
70
,
115319
(
2004
).
141.
O.
Muscato
and
V.
Di Stefano
,
Commun. Appl. Ind. Math.
10
,
20
(
2019
).