{'fraction_non_alphanumeric': 0.021386074931661038, 'fraction_numerical': 0.012006217505493917, 'mean_word_length': 4.762198888202594, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anisotropic work function of elemental crystals

The work function is a fundamental electronic property of a solid that varies with the facets of a crystalline surface. It is a crucial parameter in spectroscopy as well as materials design, especially for technologies such as thermionic electron guns and Schottky barriers. In this work, we present the largest database of calculated work functions for elemental crystals to date. This database contains the anisotropic work functions of more than 100 polymorphs of about 72 elements and up to a maximum Miller index of two and three for non-cubic and cubic crystals, respectively. The database has been rigorously validated against previous experimental and computational data where available. We also propose a weighted work function based on the Wulff shape that can be compared to measurements from polycrystalline specimens, and show that this weighted work function can be modeled empirically using simple atomic parameters. Furthermore, for the first time, we were able to analyze simple bond breaking rules for metallic systems beyond a maximum Miller index of one, allowing for a more generalized investigation of work function anisotropy.


Introduction
The work function (Φ) is an electronic surface property of crystalline solids and is crucial to the understanding and design of materials in many applications. It can be directly applied to the engineering of device specifications such as the Schottky barrier of semiconductor junctions or the thermionic currents of electron guns. Furthermore, it has been used to guide the engineering of interfacial interactions between metals and monolayer structures for nanoscale self-assembly. 1 The work function is also an important parameter in characterization techniques where it can influence the tip tunneling current of scanning tunneling microscopes or correct the binding energy in photo-electron spectroscopy (PES).
The work function has also been explored as a parameter for materials design. For example, previous experimental and computational investigations of Ni-alloys by Lu et al. 2 , 3 have established a correlation between the work function and various mechanical properties such as toughness, hardness, ductility and bulk modulus. A more recent study using firstprinciple calculations found similar correlations for elemental crystalline solids. 4 The work function has also been proposed as a possible parameter for the desorption rate of surface adsorbates. 5 Calculated work functions of hcp materials have also been used to screen for more effective metallic photocathodes. 6 Much effort has also been devoted to modelling Φ itself. Michaelson 7 and Miedema et al. 8 , for example, were previously successful in modelling the polycrystalline work function as a linear function of electronegativity. The modeling of the anisotropic work function (Φ hkl ) as a function of surface morphology and chemical environment has also garnered much attention. Smoluchowski smoothing is one such model which describes the contributions to the work function of metals as a result of isotropic electron spreading and anisotropic electron smoothing. 9 The spreading of negative charges increases the work function while the anisotropic smoothing of negative charges at the surface decreases the work function.
Smoothing increases with surface roughness (defined here as the reciprocal of the surface packing fraction 10 ) which decreases the work function. This model is supported by previous observations that the anisotropic surface energy (γ hkl ) is inversely proportional to Φ hkl via the broken bond surface density. 11,12 Similarly, the Brodie model attempts to explain Φ hkl for transition metals as a function of (bulk) electron effective mass, surface atomic radius and inter-planar distance. 13,14 A more recent model using a dielectric formalism has been proposed by Fazylov 15 that describes Φ hkl using surface roughness and surface plasmon dispersion.
An extensive database for Φ hkl would be invaluable for validating and further expanding upon these models. However, experimentally measured work functions are usually for polycrystalline specimens (Φ expt poly ) instead of single crystals. An example of this is the extensive collection of experimentally measured Φ expt poly for 66 polycrystalline elemental solids compiled by Michaelson 16 . Though measurements for anisotropic Φ hkl are not uncommon, values often vary due to the many techniques used or non-standardized methods of implementing the same technique (e.g., PES). 17,18 The sparsity of Φ hkl and the lack of a comprehensive compilation with a single standardized technique makes it difficult to develop and gain insights into work function anisotropy using experimental measurements.
Here, density functional theory (DFT) has the advantage of calculating Φ hkl for a model of any specific solid surface under a controllable set of parameters, making it possible to create a standardized collection of values. Many authors have attempted such compilations for Φ hkl , which are often times accompanied by the corresponding surface energy γ hkl . 11 Here, we report the development of a comprehensive, validated database of work functions for elemental crystalline solids using DFT calculations that addresses all the above limitations in the following ways: 1. Coverage of 142 polymorphs of 72 elements, including rare earth metals.
2. Facets up to an MMI of three and two for cubic and non-cubic crystals, respectively, are considered.
3. Common reconstruction schemes, such as the missing-row (110) fcc and the diamondtype reconstructions, 26 have been taken into account.
We validate our computed work functions with past experimental and computational data for both Φ hkl and Φ expt poly . Finally, we will discuss trends in the work function of the elements, and develop a predictive empirical model for Φ expt poly . The work function is defined as the energy barrier required to move an electron from the surface of a solid material into free space, as given by the following expression:

Methods
where V vac is the electrostatic potential of the vacuum region near the surface and E f is the Fermi energy of the slab. The energy barrier can be visualized in Figure 1 where V vac is obtained when the electron is far enough away from the surface, that the potential remains constant over a small distance in the vacuum. This method has been widely used in previous studies for calculating the work function 11,12,21,22 and has been shown to converge quickly with respect to slab thickness. 19

Modeling non-uniform work functions
For comparison to work functions obtained from polycrystalline specimens, one approach is to calculate the work function of a "patchy" surface by weighting each Φ hkl by the area fraction of its corresponding facet 5,27 as follows: where A hkl and f A hkl are the total area and the area fraction of all facets in the {hkl} family, respectively. A polycrystal is an extreme case of a patchy surface, and as such the same technique can be applied to Φ expt poly . The PES signal of the patch with a lower work function will tend to eclipse those with higher work functions leading to an underestimated measurement of Φ expt poly . Thus, experimental values of the lowest anisotropic work function (Φ lowest hkl ) are only ∼ 90 meV lower than Φ expt poly . Because of this, it has also been suggested that Φ lowest hkl is a good estimate of Φ expt poly . 5,18,21 In this study, we use the facets present in the Wulff shape previously calculated by the current authors 28 as an estimate of the orientation and area fraction present in a polycrystalline sample to obtainΦ and Φ lowest hkl as estimates for Φ expt poly . As mentioned earlier, Smoluchowski smoothing describes the anisotropic work function of metals as being inversely correlated with the broken bonds per surface area. As such, we model our values for Φ hkl normalized byΦ using the ratio of broken bonds-to-bulk coordination number ( N BB CN bulk ) in a slab normalized by the surface area-to-atomic area ratio ( in order to compare across all systems: It is known that for bcc and hcp materials, the first nearest neighbors (1NN) and second nearest neighbors (2NN) are in close proximity to each other, leading to contributions from the latter to the anisotropy of surface energy. 29 Hence, when defining BB for a material, we will limit the maximum coordination number to the 1NN for fcc materials and explore both 1NN and 2NN for hcp and bcc materials. For hcp structures, we omitted Φ 0001 when investigating the effect of the 1NN as the number of broken bonds will always be 0 which is unphysical. The inverse correlation between BB and Φ hkl for each element can be quantified by the Pearson correlation coefficient (r). We will define systems with a negative linear trend between Φ hkl and the normalized broken bonds with r < −0.75 as having a strong correlation, −0.75 < r < −0.5 as having a moderate correlation and r > −0.5 as having a weak correlation. Only ground state metallic fcc, bcc and hcp systems were explored under this context. We used the high-throughput workflow proposed by Sun and Ceder 30 and implemented by Tran et al. 28 and Montoya and Persson 31 to obtain all required data. The workflow was

Computational details and workflow
implemented using the open-source software packages Python Materials Genomics, 32 Fire-Works 33 and Atomate. 34 The work function is extracted from the calculations and inserted into the same database. To handle errors that may arise during calculations, the custodian software package was used as a wrapper around VASP together with a set of robust error handling rules. The database will be continuously improved and will continue growing as more structural data becomes available on the Materials Project. 35

Data Availability
The data can be accessed from the elemental-surface-data-focused Crystalium 36 website, as well as from the Materials Project website 37 on its detail pages for specific crystals.

Results
Due to the vast number of data points for Φ hkl when comparing to literature values, we have adopted a consistent marker shape and color scheme for ease of reference (see Figure 2) for all subsequent plots. The shape and color represents the row and group of the element in the periodic        110) and (210)). An exception to this is  41 reported that Φ Al 100 > Φ Al 110 ∼ Φ Al 111 , which is consistent with our results (see Table 2 for values).
Our values for Φ hkl are in excellent agreement with those calculated using PBE, LDA and RPBE (as shown in Figure 4  ). As found in our previous work, only Pt, Au and Ir have significantly lower surface energies for the (110) missing-row reconstruction compared to the unreconstructed surface, which is in agreement with experimental observations. In general, we find that reconstruction leads to a relatively small increase in the work functions, though the three fcc metals exhibiting a thermodynamic driving force to reconstruct also have the largest work functions.

Discussion
Periodic trends in the work function   (Pt group). When plottingΦ against the weighted surface energy (γ) in Figure 6(b), we observe a split between elements above and below group 8. A similar parabolic trend when plottingγ against group number results in a maxima observed at group 7 rather than 10.
The position of these parabolic peaks are related to the increasing cohesive energy resulting from the increasing number of half-filled d-orbitals as well as the width of the electronic s, p and d bands. For a more in-depth discussion, the interested reader is referred to the references herein. 20, 42 We further note that cohesive energy and thus surface energy are strongly correlated with mechanical properties which suggests a cohesive energy origin to previously observed trends betweenΦ and mechanical properties of transition metals. 2,4,43 To our knowledge, this work represents the first time DFT has been used to calculate the work functions of the lanthanides. Figure 6(c) plotsΦ against the group number. A gradual decrease is observed for the half-filled lanthanides from Ce to Sm with a sharp decrease for Eu. The latter half of the lanthanides has a relatively constant value from Gd to Tm with a sharp decrease for Yb. Afterwards, a sharp increase is observed for Lu. These trends are consistent with trends of the cohesive energies of the lanthanides (see Figure 10), which in turn may be attributed to the gradual filling of the 4f orbitals. The two lowestΦ are observed when the 4f orbitals are half filled (Eu) and completely filled (Yb), and these two elements also have the lowest cohesive energy and melting point among the lanthanides.
With the exception of Ce, the computedΦ underestimates the experimental values of Φ expt poly for the lanthanides with a standard deviation of 0.136 eV from Ce to Yb (see Table 1).

Discrepancies in the comparisons
In general, our computed work functions are consistent with previous computational studies. 11,21,23 It is well known that the GGA(PBE) functional underestimates the intermediate range van der Waals (vdW) forces and Fermi energy while having no long-range vdW forces, which generally leads to an underestimated work function. Although LDA generally yields values closer to experiment than GGA(PBE), this agreement is due to the various errors inherent in LDA that work in tandem to provide an error cancellation. 23 LDA/GGA are also known to have errors associated with overbinding/underbinding leading to smaller/larger cell volumes, lattice parameters and atomic distances, which can in turn influence surface properties. Larger atomic distances will decrease 2NN contributions to γ hkl and Φ hkl . This effect is especially prominent in the refractory metals Mo, Ta, Nb and W where the work function is shown to be dependent only on 1NN, but not 2NN (see Table 3). Without 2NN contributions, the work function will be severely underestimated in these metals when compared to experimental values, thus explaining the increasing deviation for refractory metals. Furthermore, the surface energy scales well with 2NN for Nb and Ta (which have the greatest deviation in Φ hkl of the four refractory metals), indicating that underbinding in GGA(PBE) is more consequential for work function than it is for surface energy.
Experimental error is also a potential source of discrepancies in our comparisons.Φ As is significantly higher whileΦ La andΦ Se are significantly lower than that of their corresponding values for Φ expt poly . Experimental values for these particular elements were taken from Michaelson 16 where surface contamination could lead to inaccuracies of up to 0.5 eV 7 in the reported measurements. Φ As poly is also known to range from 3.75 to 5.4 eV which our value for Φ As lies between. Furthermore, the value of Φ Se poly was also determined using a photoelectric method which is known to yield erroneous values of work function for semiconductors. In Although the calculated Φ lowest hkl has previously been suggested as a good approximation of experimentally measured Φ expt poly , we have shown that the Wulff-area-weightedΦ provides a much closer estimate. This suggests that the eclipsing effect of lower work functions in PES signals of patchy surfaces may not be as prominent as once thought. Despite this, the values of Φ expt poly are higher thanΦ by an average value of 0.18 eV. Kawano 5 has argued that the weighted work function is more likely to follow a Boltzmann distribution with a higher value than that provided by Equation 2. In this context, temperature becomes an important factor in determining the work function. 43,44 However, it is unaccounted for in our calculations which are assumed to take place at 0K.

Effect of reconstruction on work function
In general, the work function of the reconstructed facets are slightly larger than that of the unreconstructed facets. This can be due to the exposure of the     The comprehensive data set presented in this work affords us the ability to develop more robust models for the work function of the elemental metals. It has previously been wellestablished by Michaelson 7 and Miedema et al. 8 that the work functions of metals have a positive linear relationship with the electronegativity χ of the metal. This may be explained by the fact that χ is a measure of how strongly electrons are bounded to the atom, and hence, the higher the χ, the greater the energy needed to bring an electron from the bulk to the free vacuum (Φ). Nevertheless, as can be seen from Figure 8(a), it is clear that χ only explains R 2 = 85.5% of the variation inΦ across the metals.
We carried an investigation of the relationship betweenΦ and various atomic properties.
As can be seen from Figure 8(a), a strong, albeit non-linear, negative relationship is observed betweenΦ and the metallic radius R. From Gauss' law, the potential inside an infinite charged plate is proportional to the bulk charge density times the square of the thickness of the plate, i.e., it scales charge per unit length of material. We postulate that the average work function is proportional to the electron density per unit length, similar in spirit to the traditional jellum work function model for metals. 48 We performed a linear regression ofΦ against χ ratom , where χ is related to the electron charge contributed per atom (in line with previous models) and r atom = 3 √ V atom , where V atom is the unit cell volume per atom. As shown in Figure 8(b), this optimized modelΦ = 1.55 + 3.76 χ ratom exhibits a much improved prediction accuracy forΦ, with a very high R 2 of 0.933 and a small SEE of 0.226 eV.

Conclusion
In conclusion we have constructed the largest database of anisotropic work functions to date. We have validated our database by comparing to both experimental and computational results from the literature and by confirming previously observed trends. In addition, we have also developed a technique for estimating the work function of a polycrystalline specimen using the Wulff shape and showed that it is a significantly more accurate estimate for experimental polycrystalline values than the lowest anisotropic work function. Using this large dataset, we have also extensively probed well-known empirical relationships for the work function, such as the Smoluchowski rule, and developed a substantially-improved prediction model for the work function of the metals from atomic properties such as the electronegativity and metallic radius.     1stN N 2ndN N 1stN N 2ndN

Declarations of interest:
None