Charge density and screened potential near a (5,5) capped carbon nanotube from density functional computation

This paper continues an analysis described earlier¹ (and here identified as MEB), which used the order-N program onetep² to solve by density functional theory (DFT) the distributions of potential and charge density in a capped (5,5) carbon nanotube (CNT). Further inspection of the results from that work, in particular, the axial potential, has found details consistent with the molecular orbitals that were previously deduced and shows how the exchange potential of the emitter produces its work function.

The analysis of field emission by Fowler and Nordheim³ (FN) assumed a simple approximation for the potential barrier outside an emitting surface and deduced a relation between the applied electric field and the emitted current. This theory (with later corrections, refinements, and approximations) is the basis for an experimental method that is widely used to test whether an observed current is due to field emission. Many experimental results are consistent with FN theory, which predicts that transmission is by tunneling through the barrier (rather than by surmounting it).

However, FN also modeled the potential inside the emitter with an even simpler approximation of a uniform value. During the century after the publication by FN, the theory of waves in crystals has evolved to show that in an infinitely periodic material, the electric potential is also periodic; further, time-independent orbitals are possible only for electrons of specific energies. Even so, the calculation of a wave function to define emission from many states, in an arbitrary emitter material and geometry, and in a potential that is consistent with all occupied states, would clearly need very extensive computation. Thus, instead of solving for a wave function, it is more convenient to solve by first principles for the distribution of particle densities.

DFT^4,5 defines the distribution of “noninteracting electrons” that move only in an effective potential (here called V_T) that is the sum of (a) the classical electrostatic potential (V_H) due to electrons and nuclei, (b) a contribution (V_X) produced by the quantum condition known as “exchange,” and (c) a correction (V_C) known as “correlation.” The last two are often identified together as (V_XC). The total V_T is used in a version of Schrödinger’s equation, whose eigenfunctions define the permitted orbitals for the system. Solution by DFT can thus identify many states (each with charge density consistent with the total potential), together with their occupations. The total charge density is the sum of the intensities of the occupied orbitals.

The exchange condition, an application of Pauli’s principle⁶ from quantum physics, implies that electronic charge in the emitter is expanded relative to the purely electrostatic prediction. The exchange potential is a defined function of the self-consistent charge distribution and that charge is unlikely to have a surface that is simply shaped or clearly defined. Hence, an exchange calculation in DFT is likely to be more accurate than an approximation by use of image theory.

Here, we present more details from the time-independent cluster (or single-cell) calculation with the zero applied field, as reported by MEB. The electronic charge in the end cap is shown to diffuse around the atomic cores. The axial electric field shows small discontinuities; these correspond with local increases in charge density that form small perturbations on the exponential fall at the entry to the barrier. The locations of these perturbations are consistent with reflection of the orbitals reported by MEB. We also present the axial distribution of components of potential, showing that the work function for our structure is due almost entirely to the exchange condition in the emitter (as found by Bardeen⁷ for sodium). This implies that a strong retarding field, pummeled by the incident electrons, exists at the Fermi equipotential just outside the cores in the CNT cap, even in the absence of the applied field.

As usual in DFT, the term “external field” is used here for the field of a source such as a positive nuclear charge in the emitter that is distinct from the electron distribution. Where a background electric field is added by external electrodes or boundary conditions, it is described here as an “applied field.” A “state” refers to a molecular orbital of specific energy, with two electrons per state. Spin effects were not calculated.

The notation for variables follows that defined earlier.⁸ The axial coordinate z increases from the emitter into the barrier. The total electrostatic potential (energy) is denoted by V_T(r, θ, ϕ), or by V_T(z) where a 1D model is sufficient. An electric field F is treated as having units of voltage/distance. Where F_z(z) accelerates electrons in the positive z direction, F_z(z) is defined to be positive. At the start of the potential barrier, electrons are retarded but ∂V/∂z is positive, so in general e F = −grad V. Any difference between the chemical potential and the Fermi level of energy for the system is ignored in this paper. The Fermi level for the system is denoted by E_F. Potentials are quoted relative to E_F except where otherwise specified. Here, the equipotential surface at E_F in physical space is called the Fermi equipotential; it lies slightly outside the surface defined by the atoms in the nanotube and cap. The location at which the Fermi equipotential meets the CNT axis is denoted by z_F; this is the start of the potential barrier as seen by electrons with the energy of the Fermi level. On the axial scale used for this computation, z_F was located at about z = 0.08 nm. The notation used in Fig. 6, quoted from the thesis of S. M. Masur,⁹ is described in Sec. IV B 3.

A brief outline of DFT was given in the paper by MEB. In DFT, the force on a “non-interacting” electron is − (∇V_T), where V_T includes both the electrostatic potential and the contribution of V_XC from the exchange and correlation relations. As V_T is computed to be consistent with the charge density and the exchange condition, it cannot be chosen arbitrarily but is assembled from the components mentioned in Sec. I. In practice, V_XC is estimated with a functional such as the local density or generalized gradient approximation.

Computations by onetep for a capped (5,5) (CNT) have been described earlier.^1,6–8 The basis set used is defined by the orbitals themselves, which are generated by iteration, together with their occupations.

In the computation for our CNT cap with a generally applied field, MEB used a single (non-periodic) cubic “DFT cell” of side 5 nm and containing a tip section of the CNT. The standard method of DFT solution does not need to specify boundary conditions (other than matching) in any direction in which the system is periodic, but our single-cell calculation needs either the potential or its normal derivative to be specified over the whole boundary surface, together with a value for the net charge enclosed in the cell. With these conditions, the solution is expected to be numerically stable, while if only the axial potential were specified, it would be unstable.

The method used by MEB for fixing the boundary conditions around the DFT cell in the applied field was, as far as we know, novel. A separate macroscopic system was defined to contain cathode and anode planes plus a conducting cylinder and hemisphere, in contact with the cathode plane, to simulate the full length and shape of the capped CNT. The diameter of the simulated CNT was equal to that of the Fermi equipotential surface as found (iteratively) in the DFT computation. Some of the mesh points for the macroscopic system were defined at positions corresponding to the surface of the DFT cell. The macro system was solved by a classical solver with classical conditions at the surface of the simulated emitter. Interpolation between mesh points on the cell boundary then provided boundary potentials at the same points for the DFT solution.

MEB found that the greatest contribution to the axial potential from exchange occurs where electron density is highest: near the end pentagon of carbon atoms in the cap, as in Fig. 6. For an ambient temperature that is not zero, onetep initially assigns the occupation of states in the emitter according to Fermi–Dirac statistics and uses ensemble DFT to minimize the Helmholtz free energy.¹⁰ 

The rest of this paper describes results found with the zero applied field. It is a benefit of DFT that results are produced directly in physical coordinates. The previous work¹ produced files of potential and charge density at 416 × 416 × 416 points (correcting a typo in MEB) on a mesh in the DFT cell. These files were combed for the results shown here.

We then tested other functional relations. We estimated the slope dV_T/dz for z > z_F from the differences between computed values of V_T(z) over intervals of 2 pm. On plotting dV_T/dz against (V_T− E_F), a piecewise linear graph appeared as in Fig. 3, which also includes a plot of (z− z_F) against (V_T− E_F). The separate straight sections of dV_T/dz have slopes and intercepts (found from trendlines) as shown in Table I. The discontinuities in dV_T/dz appear at equal intervals in z of 0.050 nm. Each of the two discontinuities visible in Fig. 3, separating sections (i) from (ii) and (ii) from (iii), occurs at the same z as a small local increase in the electronic charge density (see Sec. IV B 5). The calculated dV_T/dz has a minimum at z ≈ 0.46 nm.

The charge density is expected to fall off exponentially with the emitter radius where the potential reaches the vacuum level, as can be seen in Fig. 1(b). The same results are plotted on expanded scales in Fig. 4 and show more clearly that some of the charge penetrates into the barrier region. The presence of charge at locations out to z = 0.2 nm (or 0.4 nm as shown in Fig. 1(b) at r ∼ 0.9 nm) implies that this charge, which may be a steady state of outgoing and reflected current, has an energy up to 3 eV or more above E_F. This point is discussed in Sec. IV B 5.

When the computed values of V_T are used to plot the left side of (3) against (z − z_F) as in Fig. 5, the changes in slope corresponding to the discontinuities in Fig. 3 are not obvious. Those discontinuities are discussed in Sec. IV B 5. The work function for Fig. 5 was set as 4.581 eV to include a small axial dipole from inside to outside the nanotube wall (visible in Fig. 6).

Equation (4) shows that (V_vac − V_T(z)) is computed to fall exponentially as z increases. At first sight, this is typical of screening, as discussed for a free electron gas (in a metal) within Chap. 17 of Ashcroft and Mermin.¹⁴ However, Eguiluz et al. report¹⁵ a more relevant calculation by the GW method¹⁸ (believed to be accurate for correlation at a long range) for a surface with “r_s = 3.93” showing that “V_XC becomes imagelike outside the surface.” The question of what dependence on z in the barrier is realistic is not discussed further here, because onetep itself imposed certain limits which may have restricted its computation of V_T(z).

According to the description of onetep,¹⁰ the aim is that the code achieves “linear-scaling [computational] cost with near-complete basis set accuracy.” The program subjects the orbitals and the density kernel (both constructed at run time) to “localization constraints, so as to achieve the variational freedom of a large basis set while retaining the … small matrix sizes of a minimal basis” [of atomic orbitals]. For the computation reported here,

• the maximum radius of orbitals was set at 12 Bohr radii (∼0.64 nm);

• the kinetic energy cutoff was set at 1000 eV;

• the local density approximation was used for most calculations.

Further tests are planned to vary these limits, for example, by using a method of known long-range accuracy or by varying the settings for computation.

Figure 6 shows the computed total potential V_T(z) on the axis near z = z_F for the zero applied field, together with contributions to V_T from its components. These contributions were plotted relative to the vacuum level, while the potential scale is shown relative to the Fermi level. The total V_T consists of components named in Fig. 6 as Hartree (electrostatic), Ion-ion, and XC, and identified here for brevity as V_H, V_c−c (for core-core), and V_XC. The plot labeled “Ion-ion potential” shows (V_vac − V_c−c), to ensure that the negative of the slope correctly suggests the direction of the force on an electron (as described in the  Appendix). The short-dashed line identified here as V_H is the sum of the (negative) potential V_e−c induced by electron-core attraction and the (positive) V_e−e due to electron-electron repulsion. The distributed electrons provide partial screening, so the combined axial potential V_H = V_e−c + V_e−e is negative relative to the vacuum level. The minimum in the valley in V_H is at z = −0.09 nm (close to the end pentagon of C cores in the cap) and this valley is almost completely confined to z < z_F. Outside z_F, V_H rises locally above E_F by about 0.2 eV, possibly caused by the curvature of the end cap of the CNT.

Independently of the Coulombic forces, the exchange condition keeps electrons further apart on average than if the condition did not apply. Thus, exchange lowers the total potential below V_H, which shows the “direct” electrostatic potential alone. The effect of exchange is represented here by the potential V_XC, which in this case goes the furthest negative where the charge density is greatest. The reduction in electron density caused by exchange around a single electron is called an “exchange-correlation hole.” It may alternatively be regarded as a localized volume of positive charge superimposed on the negative background density, thereby contributing to the potential V_XC.

Figure 6 shows clearly that V_XC rises from E_F to V_vac at z > z_F, while V_H rises to V_vac at z ≈ z_F. Thus, for this structure, exchange causes the transition of V_T from E_F to V_vac to move outward from z_F by about 0.4 nm. In the range 0 < (z − z_F) < 0.42 nm, the total energy V_T(z) resembles the exchange contribution V_XC. This parallels the effect reported by Bardeen⁷ for sodium that the potential barrier outside the emitter is chiefly due to the exchange constraint for electrons within the emitter material. The V_H shown here was calculated with exchange interaction present.

In Fig. 3, a discontinuity in axial dV_T/dz together with a slight change in its slope can be seen at (V_T − E_F) = 1.63 eV and (z − z_F) = 0.04 nm. According to macroscopic electrostatics, this behavior should be accompanied by a local concentration of negative charge at z = 0.12 nm. An earlier plot by S. M. Masur⁹ of axial electron density in zero field as a function of z for the same structure is reproduced here as Fig. 4. The inset to Fig. 4 shows that at z = 0.12 nm, a small increase in electron density is superimposed on the exponentially falling background. Also in Fig. 4, an even smaller increase in density can just be seen at z ∼ 0.175 nm. This location is close to that of the second discontinuity visible in Fig. 3 and listed in Table I. Further small discontinuities in dV_T/dz are deduced from V_T(z) and listed in Table I but are not visible in Fig. 4.

The steady-state DFT solution of Fig. 4 shows that a constant small proportion of the available charge density enters the barrier. This constant charge density is believed to be caused by a stream of electrons, each electron being reflected at the location of its own (normal) energy in the barrier. From Figs. 3 and 4, the small density (of order 0.1% of the maximum or less) that penetrates as far as (z − z_F) = 0.175 nm includes energies up to at least 3 eV above E_F. The Local Density of States plot for the end cap⁹ shows that a range of states (or orbitals) with energy 1–3 eV above E_F is permitted; but from the Fermi–Dirac distribution one would not expect these to be occupied in the emitter at room temperature. However, the substantial review of field emission by Gadzuk and Plummer¹⁶ reported experimental observation of a low-density tail of up to 4 eV above E_F, in current emitted from tungsten at 78 K. The same authors reported that the current density of these higher-energy electrons was proportional to the square of the total current, suggesting that they could be produced by two-particle scattering.

If our calculated discontinuities in axial dV_T/dz are indeed due to reflected electrons in allowed states above E_F, it is not obvious why they appear at regular intervals in z of 0.05 nm as shown in Table I. Friedel-like screening from all the Kohn–Sham states in the strong retarding field may be relevant but is not pursued here.

Further results from the earlier computation with onetep for the ground state of a capped (5,5) CNT in zero applied field show details of charge density and axial potential in the cap that have been calculated self-consistently using adaptive orbitals. The electronic charge diffuses around the atomic framework on a scale of the order of a bond length.

The potential is built from the first principles of electrostatic interaction plus quantum contributions for exchange and correlation. A plot of these components on the axis shows that in this structure, the exchange condition creates a potential change between the Fermi level and the vacuum, in the barrier region outside the Fermi equipotential. Thus, it creates the work function of the emitter and provides the initially retarding field in the potential barrier. Subject to checking the long-range calculation of correlation, the DFT computation provides a physical basis showing that the widely used Schottky–Nordheim potential¹⁷ is a good approximation for the start of the potential barrier.

The axial electric field shows small discontinuities near the start of the potential barrier. These appear at potentials that align with those calculated previously for specific molecular orbitals in the emitter. The multistate DFT calculation, thus, suggests that a small fraction of electrons incident on the barrier have energies between the Fermi level and vacuum level. This energy spread of a minority is greater than is predicted by the Fermi–Dirac function and recalls an earlier experimental observation¹⁶ with a different emitter material. These hot electrons may contribute to the current at the onset of field emission, but whether they can provide a current-voltage relation comparable to the experimental evidence of tunneling needs further consideration.

C.J.E. thanks A. Kyritsakis for stimulating discussions. The plots of Fig. 1 were produced by J.H. as part of an unpublished report¹¹ for Part III of the Natural Sciences Tripos, Cambridge, 2023. Figure 3 was generated from a file computed by S. M. Masur with onetep in the course of her research⁹ as a Marie Curie Fellow for the Simdalee2 Project. That project received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme FP7/2007-2013 under REA Grant Agreement No. 606988. Figures 4 and 6 are adapted from figures published by S. M. Masur in her Ph.D. thesis.⁹ 

The authors have no conflicts to disclose.

C. J. Edgcombe: Conceptualization (equal); Data curation (supporting); Formal analysis (lead); Investigation (equal); Project administration (lead); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (lead). J. Huns: Data curation (lead); Investigation (equal); Validation (equal); Visualization (equal).

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

In Fig. 6, the dotted potential plot identified as “Ion-ion potential” or V_c−c shows the potential on the axis produced by the interaction between the carbon cores. In the cylindrical tube of the CNT containing hexagons of atoms, each carbon atom with sp² bonding has given up a mobile π* electron.

For this capped single-wall CNT, the magnitude of V_c−c on the axis is greatest near the end pentagon of C atoms, at z ≈ −0.09 nm. Since like charges repel, any pair of adjacent cores will raise the potential nearby and repel another positive charge. If V_c−c is shown as a hill in potential, then the usual rule, that the negative of the slope of the plot shows the direction of force, will show the direction of force on another positive charge. To find the direction of force on an electron, we choose to modify the usual rule. In Fig. 5, the V_c−c plot has been inverted, so that application of the usual rule (whether consciously or not) will correctly show the direction of force on an electron due to ion-ion interaction. Hence, with V_c−c as shown in Fig. 5, the maximum deviations from the vacuum level satisfy

$V T= V H+ V XC− V c − c at z≈−0.09 nm.$