Beyond ideal two-dimensional metals: Edges, vacancies, and polarizabilities

Recent experimental discoveries of graphene-stabilized patches of two-dimensional (2D) metals have motivated also their computational studies. However, so far the studies have been restricted to ideal and infinite 2D metallic monolayers, which is insufficient because in reality the properties of such metallic patches are governed by microstructures pervaded by edges, defects, and several types of perturbations. Here we use density-functional theory to calculate edge and vacancy formation energies of hexagonal and square lattices of 45 elemental 2D metals. We find that the edge and vacancy formation energies are strongly correlated and decrease with increasing Wigner-Seitz radii, analogously to surface energies. Despite a radical reduction in atomic coordination numbers, the 2D and 3D vacancy formation energies and work functions are nearly the same for each metal. Finally, static polarizabilities reveal a clear cubic dependence on bond length. These trends provide useful insights when moving towards reality with elemental 2D metals.

The experimental evidence for 2D structures composed of metal atoms has motivated much computational research. Multiple elements and 2D lattices have been studied including Au, Ag and Cu monolayers [26][27][28], transition metal monolayers [29], and our recent study of elemental monolayers from 45 metals in three 2D lattices [30]. However, so far most of the computationally studied free-standing elemental monolayers have been ideal and periodic in the plane of the atoms. These calculations are relevant for sufficiently large, ordered systems. Yet under realistic conditions 2D systems will have defects, such as edges and vacancies [31]. While the edges can be stabilized by supporting materials, edge formation energies provide information about the stabilities of finite systems compared to the periodic ones. While the vacancy formation energy is related to the stability of a * pekka.koskinen@iki.fi 2D structure, it is also connected to the atom mobility and is useful in identifying promising elements for 2D liquids [32,33]. Realistic systems can also have perpendicular perturbations. For example, 2D systems can host adsorbates [34], can be grown on surfaces [35][36][37][38][39][40][41][42], or can be a part of a layered heterostructure [43]. In these cases the interactions perpendicular to the monolayers changes the properties of the ideal free-standing 2D lattices.
In this paper we aim to address how the finite size, vacancies, and simple perpendicular perturbations affect the properties of ideal free-standing monolayers of 2D metals. Using a density-functional theory (DFT) approach, we calculate the edge and vacancy formation energies for 45 metals (Fig. 1) in hexagonal and square geometries (Fig. 2). We correlate these properties with the ones of conventional 3D bulk structures and find that the edge formation energies are related to the surface energies, both decreasing with increasing Wigner-Seitz radii. Despite the drastic changes in coordination numbers, the 2D vacancy formation energies for many metals are close to 3D vacancy formation energies. Further, since vacancy can be considered to consist of a round edge encircling the missing atom, the vacancy formation energies are related to the edge energies. As a measure of sensitivity to perpendicular perturbations, we also consider the static polarizability of 2D monolayers and find that the polarizability per atom increases with increasing 3D bond length. The work functions of hexagonal 2D films are relatively close to those measured for polycrystalline 3D

II. COMPUTATIONAL METHODS
The edge and vacancy formation energies were obtained from total energies calculated with the densityfunctional approach as implemented in the GPAWcode [44,45]. For consistency with respect to earlier work, the exchange and correlation energies were approximated with the Perdew-Burke-Ernzerhof (PBE) functional [46]. Also previously converged computational parameters and lattice constants were used [30]. All structures were calculated without relaxation using ideal bond lengths. The plane-wave cut-off was 800 eV and 5Å vacuum region separated atoms from the non-periodic unit-cell edges. The atomic ribbons were modeled with 1×12×1 Monkhorst-Pack k-point sampling [47,48] and vacancy formation energies were calculated with a constant k-point density in the atomic plane. The polarizabilities of monolayers were calculated with dipolelayer corrections. Jellium calculations were done as spincompensated and rest as spin-polarized.

A. Edge energies decrease with increasing
Wigner-Seitz radii We begin by considering the edge energy of all 45 metals in hexagonal and square geometries. Spin polarization is taken to account but since most systems are nonmagnetic we focus on other properties, starting with the edge energy [49] defined as where 2D is the energy per atom in periodic 2D structure, E r (N ) the energy of a ribbon with N atoms and L edge the total length of the edge i.e. twice the length of the computational cell in the periodic direction. To get the edge energy one could choose a ribbon of specific width, calculate its energy, and subtract the corresponding 2D lattice energy. In this case the result would depend on the width of the chosen ribbon. Fortunately, this dependence can be removed with the ansatz from which the edge energy is obtained by calculating ribbons of varying widths and fitting the properties edge and 2D simultaneously. Further, comparing the fitted 2D lattice energy 2D to one from a periodic calculation gives a convergence test for the edge energy. This method is a 2D analog of a similar approach for determining surface energies [50]. The energy of a ribbon as a function of its width is linear already for very narrow ribbons, indicating that the ansatz (2) holds already for small N .
To rationalize this ansatz we consider a simple model of non-interacting electrons in a 1D box. To model a ribbon with varying width L x we set the number of particles N proportional to the width of the well N = λL x . The energy as a function of N in atomic units is which displays the observed linear behavior for large N . Further, second term on the rightmost side is independent of N and corresponds to the edge energy. We expect that similar calculation with a 3D box introduces some complications, but leaves the general trend unaffected [51]. The edge energies obtained this way are in general high near the middle of the d-series and particularly high for 5d-metals in hexagonal structures (Fig. 3). As discussed in our previous work [30], metals near the middle of the dseries have occupied bonding orbitals and unoccupied antibonding orbitals [52]. This makes their bonds stronger and edge energies higher. The trend is qualitatively similar for both hexagonal and square lattices. The previously reported value of 0.2 eV/Å for Au agrees with our result [34].
Next, we consider edge energies as a function of the Wigner-Seitz radius r s defined by the equation where N is the number of valence electrons in volume V . This new viewpoint emphasizes how the edge energies span almost two orders of magnitude and display roughly monotonic decrease with increasing r s (Fig. 4). The trend holds especially well for simple metals. Most important, similar behavior has been reported for surface energies [54]. We conclude that metals with high surface energies will have high edge formation energies, a trend that calls for closer inspection. This trend can be rationalized using jellium ribbons. A jellium ribbon has a finite width L x and thickness L z but length L y that approaches infinity. The r s for a ribbon is defined as where N is the number of electrons in the unit cell and L y the length of the unit cell in the periodic direction. We show that DFT energies for ribbons with varying L x , L z and r s are reasonably well described by a liquid drop model [56]. This will imply a simple connection between edge and surface energies. The liquid drop model gives the total energy of the jellium ribbon as (6) where u is the energy density of bulk jellium and σ the surface energy. Using equation (5), the energy per elec- where b is the energy per electron for bulk jellium. We calculate ribbons with thickness of single fermi wave length λ f and change the width L x from λ f to 5λ f . and fitted surface energies (Fig. 5a) agree with previously reported values [55]. Therefore the energy of a ribbon is reasonably well described by the liquid drop model. Since the energy contribution from the jellium edges is 2L y L z σ, it is natural that the edge energies have similar r s trend  (7) (circles). The dashed line shows the bulk energy given by equation (8) as the surface energies. Nevertheless, the correspondence between edge and surface energies is somewhat surprising considering that the ribbons are only monolayer thick.

B. 2D and 3D vacancy formation energies correlate well
Next we study the vacancy formation energies. In practice they depend on the vacancy density, but this dependence can be removed using a method analogous to the one used to calculate the edge energy. We calculate five monolayers of different size, each with a single vacancy. The vacancy formation energy v is then obtained from a fit where E v (N ) is the energy of the monolayer with N atoms and a single vacancy. As in the previous section, we can compare the 2D from the fit to the value from periodic calculation to confirm convergence. The resulting vacancy formation energies range from nearly zero to almost 4 eV (Fig. 6a). Comparison to calculated 3D bulk vacancy formation energies shows that for many metals the vacancy formation energies in 2D and 3D have similar values (Fig. 6b). This supports our previous observation of 2D bonds being stronger than the 3D bonds [30].
Moreover, the vacancy formation and edge energies dis-  play qualitatively similar trends. The vacancy formation energies are generally high near the middle of d-series and again the highest for the metals near the center of 5d-series. This can be understood by considering the formation of a vacancy as a formation of finite-length edges to the monolayer. To confirm this interpretation, we plot the vacancy formation energy as a function of the edge energy times the length of the formed edge. We simply assume that the vacancy is circular with edge length πd, where d is the nearest neighbor distance. While the vacancy formation energies from this simple approximation are generally overestimated they are fairly close to the directly calculated ones (Fig. 7).

C. 2D and 3D work functions are nearly identical
Next, we consider the work functions of hexagonal 2D monolayers of all 45 metals. The work function is a global reactivity descriptor related to the cost of electron removal [59,60]. As a result, we find that the calculated work functions of 2D monolayers are close to the work functions measured for polycrystalline 3D bulk samples (Fig. 8). This is somewhat surprising since the polycrystalline samples contain crystals with different sizes and surface structures. For nanosized systems the variation in size is known to lead to oscillations in the work function. These changes are due to quantum size effects that arise when some sample dimensions are close to the electron wave lengths [61][62][63]. While the quantum size effects are important for accurate calculations of small samples, the correlation between monolayer and polycrystalline work functions indicates that the quantum size effects do not drastically alter the work function. Therefore, the work function of a 3D system is a reasonable first approximation for the work function of a monolayer.

D. 2D polarizability increases with 3D bond length
To keep our approach as generic as possible, we also consider a perturbation of a constant electric field perpendicular to the atomic plane. We quantify the results in terms of the static polarizabilities calculated for hexagonal 2D metals only. Calculation of polarizabilities is motivated because the polarizability of a molecule is related to its reactivity. For example, larger polarizability of a molecule correlates with stronger physisorption [64]. Similarly, larger surface polarizabilities increase physisorption energies. We calculated the polarizability α by applying an electric field perpendicular to the atomic plane and varying the field strength E from 0.05 to 0.5 V/Å. The polarizability was then given by a fit where p is the dipole moment. The resulting polarizabilities per atom for the monolayers range from 1 to 10 A 3 , which are significantly lower than the corresponding polarizabilities of free atoms that range from 5 to 60Å 3 [65] (Fig. 10). In order to visualize the difference between free atom and a monolayer, we define the local dipole moment where n is the electron density and z is the direction perpendicular to the atomic plane. While free Na atom has large dipole moment near the nucleus, p loc is almost constant for a hexagonal Na monolayer (Fig. 9). To understand the difference between free and bound atoms, we consider a simple model. The static polarizability per atom is given by the E → 0 limit of the equation [66] where ML is the energy per atom for a monolayer. For a periodic system the energy per atom in terms of cohesive energy coh is where free is the energy of a free atom. Therefore, if coh (E) and free (E) are known, the polarizability of an extended system can be calculated. In practice, exact expressions for coh (E) and free (E) are not known. However, by approximating coh (E) the polarizability of many-atom system α ML can be calculated in terms of polarizability of free atom α free . A simple approximation for coh (E) can be obtained by considering interaction between aligned dipoles. Assume that the application of an electric field gives rise to an equal dipole moment for each atom. The dipole-dipole interaction energy is dip (r) = p 2 /r 3 , where p is the dipole moment and r the distance between dipoles. Since p = α ML E, the energy as a function of the applied field is dip (r, E) = α 2 M L E 2 /r 3 . Approximating the cohesion energy as a sum over all dipoles gives where the summation is over (n, m) ∈ Z 2 \ (0, 0) and a and b are the lattice vectors. Taking the second derivative with respect to E gives where d is the bond length and S is the lattice sum S = (n,m) =(0,0) as discussed in ref [67]. Solving for α ML gives According to equation (17), the polarizability of an extended system is always smaller than that of a free atom because ML (E) varies slower than free (E) due to the cost from from dipole interactions. Also α ML → α free when S → 0 or d → ∞, as expected. If we approximate the polarizability of a free atom to be proportional to the size of the atom [68,69], which in turn is proportional to the bond length, the polarizability of a monolayer becomes proportional to the cube of the bond length. This cubic dependence is indeed what we observe (Fig. 10).

IV. SUMMARY AND CONCLUSIONS
We studied the properties of edges and defects in mono-atomic free-standing 2D structures composed of metal atoms by density-functional theory. We considered 45 metals in hexagonal and square geometries and calculated their edge and vacancy formation energies. To keep our results general, we removed ribbon size and vacancy density dependence by utilizing a linear ansatz with a correct asymptotic behavior. We rationalized the ansatz for edge energies with a simple model of particles in 1D box. The edge energies ranged from almost zero to 0.6 eV/Å and had the highest values near the middle of d-series. Further, they decreased almost monotonically with increasing Wigner-Seitz radius, especially for the simple metals. A similar trend has been observed for surface energies. We explained this connection and the dependence between edge and surface energies by jellium ribbons and the liquid drop model.
The 2D vacancy formation energies ranged from nearly zero to almost 4 eV and were highest near the middle of the d-series. For many metals the 2D vacancy formation energies were unexpectedly close to the 3D vacancy formation energies. This was in line with our earlier observation of 2D bonds being stronger than 3D ones. The vacancy formation energies were approximated well by the edge energies after considering the formation of a va-cancy as the formation of a hole with a round edge. This connection implies that the edge energies can be used to estimate the formation energies also for vacancies of multiple atoms. Further, the energies of flat clusters can be quickly estimated using the edge energies.
Last, we perturbed the hexagonal monolayers by a constant electric field and calculated the dipole polarizabilities. We found that the polarizability per atom is significantly lower for monolayers compared to free atoms. We gave a simple model based on dipole interactions and obtained a cubic dependence between bond length and monolayer polarizability. The polarizability is relevant for physisorption of the monolayer on substrate and for adsorbates on the monolayer since larger polarizabilities lead to stronger physisorption. The work functions of hexagonal monolayers were found to be close to work functions measured for polycrystalline samples, indicating that the changes in surfaces structures and quantum size effects do not drastically alter them when going from 3D to 2D. The work function is a global reactivity descriptor that is important in charge transfer processes. These results, which have been collected to a single table (Table I) for readers' benefit, contribute to the growing field of 2D metals and especially to advancing from idealized semi-infinite systems towards more realistic finite systems with defects and interactions between the environment.