Atlas for the properties of elemental 2D metals

Common two-dimensional (2D) materials have a layered 3D structure with covalently bonded, atomically thin layers held together by weak van der Waals forces. However, in a recent transmission electron microscopy experiment, atomically thin 2D patches of iron were discovered inside a graphene nanopore. Motivated by this discovery, we perform a systematic density-functional study on atomically thin elemental 2D metal films, using 45 metals in three lattice structures. Cohesive energies, equilibrium distances, and bulk moduli in 2D are found to be linearly correlated to the corresponding 3D bulk properties, enabling the quick estimation of these values for a given 2D metal and lattice structure. In-plane elastic constants show that most 2D metals are stable in hexagonal and honeycomb, but unstable in square 2D structures. Many 2D metals are surprisingly stable against bending.

Still, 2D materials with metallic bonding have numerous potential applications [25], including catalysis [26] and gas sensing [27], which makes them an inviting research subject. Also, recent literature encourages studying 2D metals further. For example, in a recent experiment, an atomically thin iron membrane was grown inside a graphene nanopore [28]. The membrane appeared to be iron atoms in a 2D square lattice structure. Later theoretical works suggested that the iron patches are more stable as carbides, containing also carbon [29,30]. Similarly, free-standing monolayer-thick zinc oxide [31] and copper oxide [32] membranes have been observed inside graphene nanopores. In addition to these 2D structures, * pekka.koskinen@iki.fi theoretical works have predicted the existence of stable 2D gold [33], silver [34], and copper [35] membranes. Further, a computational study found that Au membranes may form in graphene nanopores by way of only small energy barriers for Au diffusion [36]. Previous works also show that Au [37] and Pt [38] readily diffuse on graphene. While these examples are about solid 2D structures, simulations have predicted even the existence of a 2D liquid [39]. A gold membrane in a graphene nanopore was investigated with molecular dynamics simulations and a solid-to-liquid phase transition was observed. Later theoretical calculations extended the prediction of 2D liquid phase also to Cu, Ag, and Pt [40]. This phase transition is facilitated by the flexible metallic bonding, which increases the mobility of the atoms compared to materials with directional and rigid covalent bonds. Therefore, materials with flexible metallic bonding have potential to establish a new class of 2D materials with novel properties. Alas, only a few elemental 2D metals have been studied; a systematic investigation with many elements is missing.
In this paper, we provide this missing systematic investigation. We present a density-functional study of 45 elemental 2D metals (Fig. 1). Each element is considered in hexagonal, square, and honeycomb structures. These structures are chosen to obtain atoms with different coordination numbers. We compare the cohesive energies, bond lengths, and elastic parameters between the 2D structures. We find that many of the properties of the 2D metals are inherited from the corresponding properties of 3D bulk. Therefore, by using the well-known 3D bulk properties, the cohesive energies, bond distances, and bulk moduli can be quickly estimated for a given 2D structure and metal. The same properties correlate linearly also among the studied 2D lattice structures.

II. COMPUTATIONAL METHODS
The calculations, based on density-functional theory (DFT), were done in the plane-wave mode of the GPAW code [41][42][43], using an 800-eV plane-wave cutoff energy and default setups. Since the calculated systems were hypothetical, a well-known nonempirical exchange and correlation functional was preferred. The Perdew-Burke-Ernzerhof (PBE) exchange and correlation functional [44] is nonempirical, computationally inexpensive, and reproduces the bulk cohesive energies well [45]. Although the PBE functional is known to reproduce bond lengths less accurately, we expect that the general trends are not that sensitive to possible errors in individual systems. Therefore, the PBE functional was used for all calculations. Three different 2D structures were studied: the hexagonal (hex), square (sq), and honeycomb (hc) lattices (Fig. 2). The hexagonal and square lattices were modeled using two atoms and the honeycomb lattice using four atoms in the unit cell. All calculations had 5-Å vacuum regions for nonperiodic directions. The flat 2D structures used a 12×12×1 and the 3D bulk used a 12×12×12 Monkhorst-Pack k-point sampling [46,47]. Bent structures used a 1×12×1 sampling. All calculations were spin polarized and convergence was checked with respect to vacuum layers, k-point grids, and plane-wave cutoff.

A. 2D cohesive energies correlate between structures
We begin by considering the 3D bulk and the three ideal 2D structures for all 45 metals. Most structures are nonmagnetic so we focus on other properties, beginning with cohesive energy. The cohesion is determined by calculating the energies of systems with ideal bonding angles and minimum-energy bond lengths. These are then subtracted from the energies of free atoms. The resulting values of 2D structures are compared relative to one another and to the 3D cohesion. The cohesive energy is defined as where E atom is the energy of a free atom and E/N is the energy of the structure divided by the number of atoms in the unit cell. With this definition, the larger the cohesion energy, the more   energy is required to sublimate the system. In order to validate the computational method, the bulk cohesive energies are calculated for all 45 metals and compared to the corresponding experimental values [48]. As expected, the calculated 3D bulk cohesive energies follow the experimental values well, with a mean absolute error of 0.28 eV. This indicates that the used numerical method is sufficiently accurate to identify trends between 2D geometries. The mean absolute error for 3D bulk cohesion is also in line with a previously reported value of 0.24 eV [45].
The calculated cohesive energies for all considered structures and metals are shown in Fig. 3. They range from nearly zero (Hg) to almost 9 eV (5d metals W and Os). The weak cohesion of Hg is in part caused by its notoriously strong electron correlations, difficult to capture by any functional [49]. Nevertheless, as discussed before, the PBE functional is otherwise expected to capture the correct trends in cohesion energies. The cohesion is, in general, the highest in the middle of the transition-metal series. This can be rationalized by a simple model [50]: Assuming that the energy of the metal can be approximated as a sum of single-particle energies, the cohesive energy can be written as where f is the Fermi energy, 0 the energy of the d state for the free atom, and ρ the density of d states. If one further assumes that the density of d states is constant ρ 0 in an energy range w located symmetrically around 0 , Eq. (2) gives with zero energy chosen so that 0 = 0. The number of d electrons N d is obtained by integrating the density of d states to the Fermi energy and the number of electrons in a full d band by a similar integral over the entire width of the d band By combining Eqs. (3), (4), and (5), we obtain a parabolic relation for E coh as a function of N d : with the maximum cohesion in the middle of the d series. While this simple model can be used to rationalize the qualitative trend in the cohesive energies, it fails to produce quantitative agreement. For example, the model holds better for the 4d and 5d series than for the 3d series. The deviation in the 3d metals can be attributed to strong Coulomb correlations, which cause band splitting [51]. In general, the simple s metals have lower cohesion than the d metals due to the lack of d-electron contribution to bonding. Let us next look more closely how the cohesion energy is affected by the 2D structure. The energies of the calculated 2D structures correlate linearly. The calculated cohesive energy of the square lattice E sq is an approximate function of the calculated cohesive energy of the hexagonal lattice E hex , with a function E sq = α×E hex and a fit parameter α. Values for the fit parameters between the different 2D structures are as follows: E sq = 0.932×E hex , E hc = 0.773×E hex , and E hc = 0.831×E sq . Most important, the 2D cohesive energies can be linearly correlated to the experimental 3D bulk cohesive energies [ Fig. 3(b)]. While the emerging energy correlations are not perfect, they well indicate the general trends, which are as follows. First, the cohesion is the strongest for the 3D bulk structures, as expected. Second, generally the hexagonal geometry is the most stable and the honeycomb lattice the least stable. The square-lattice cohesions are somewhere in-between hexagonal and honeycomb cohesions. Previous studies on 2D iron [52], silver [34], and gold [33] membranes obtained the same relative ordering. A more recent study also found the close-packed hexagonal Au lattice energetically more stable than the honeycomb lattice [53]. Third, despite the dramatically reduced coordinations, the most stable 2D structure has about 70% of the 3D bulk cohesion. After all, most of the metals have 12 nearest neighbors in 3D bulk, but only six nearest neighbors in the hexagonal 2D lattice. Similarly, the square lattice has about 65% of the bulk cohesion with four nearest neighbors and the honeycomb lattice about 54% with only three nearest neighbors. Therefore, the average bond in a 2D structure is stronger than in the 3D bulk.

B. Equilibrium bond lengths correlate between structures
Let us continue to study the bond lengths of the metals in the three ideal 2D structures. The reported bond lengths correspond to the lowest energies and range from 2 to 5.5Å (Fig. 4). Comparison to the previous section shows that, while near the middle of the d series the cohesive energies have maxima, the bond lengths have minima. This is in agreement with the maxim of stronger bonds being shorter [54]. Also, the changes in the interatomic distances have a parabolic shape for all d series. Just as the cohesive energies, the bond lengths correlate linearly between the 2D structures. Performing fits as in the previous section, we obtain the values d sq = 0.950×d hex , d hc = 0.943×d hex , and d hc = 0.992×d sq for the bond lengths. For most metals, the bonds are the longest in the hexagonal lattice, intermediate in the square lattice, and shortest in the honeycomb lattice.
The bond lengths again correlate to the experimental 3D bulk values [ Fig. 4(b)]. Hg deviates from the general trend, again perhaps due to inaccurate description of electronic correlation. As a trend, the bond lengths are the longest when the number of nearest neighbors is the largest. In 2D, the number of nearest neighbors is six for hexagonal, four for square, and three for honeycomb 2D lattices. In 3D, the number of nearest neighbors for most metals is 12. While the cohesion is reduced by about 30% from bulk values, the bond lengths are only about 1% shorter in the hexagonal 2D geometry. It seems that, while metals do not prefer layered structures, the geometries of the individual metal layers do not change much as the thickness of the metal reduces to a monolayer. For the square and honeycomb lattices, however, the bond lengths do shrink considerably compared to the 3D bulk. Yet, the bond lengths in square lattice are only 1% longer than in honeycomb lattice. This could reflect the smaller change in the number of nearest neighbors: When the structure is changed from the hexagonal lattice to the square lattice, the number of nearest neighbors is reduced by one third, but when the structure is changed from the square lattice to the honeycomb lattice, the number of nearest neighbors is reduced only by one fourth. In the previous section, we observed that the 2D bonds are stronger than the 3D bonds. The shorter bond lengths for the 2D structures are in line with this observation.

C. Bulk moduli correlate between structures
In this section, we compare the 2D bulk moduli between different structures and find similar correlations to 3D bulk 035411-3 moduli as in the two previous sections for cohesive energies and bond lengths. The bulk moduli are determined from a fit to isotropically deformed structures. Like cohesive energies and bond lengths, the bulk moduli B correlate between different 2D structures and fittings give B sq = 0.860×B hex , B hc = 0.576×B hex , and B hc = 0.666×B sq . The 2D bulk moduli display roughly similar behavior as a function of the proton number as the cohesion energies (Fig. 5). The bulk moduli are the largest near the middle of the d series and increases from 3d metals to 5d metals. The bulk moduli for the 2D structures correlate linearly with the experimental 3D bulk moduli. Cr is an exception for this correlation, with experimental 3D bulk modulus of 190 GPa.
Since the cohesions are the strongest and bond lengths the shortest near the middle of the d series, it is reasonable that Calculated 2D bulk moduli for hexagonal (hex), square (sq), and honeycomb (hc) structures, as a function of the proton number (a) and experimental bulk moduli (b). Lines display the linear fits between calculated and experimental values [48]. also the 2D bulk moduli are the highest near the middle of d series. Shorter bonds are harder to contract and elongate, which leads to higher bulk modulus. Note that while the correlations between calculated and experimental cohesive energies and bond lengths have the same units, the bulk moduli have different units in 2D. The 3D bulk moduli are in units GPa, while the 2D bulk moduli are in units GPa nm. The results for hexagonal 2D structures are summarized in Fig. 6.

D. Elastic constants
So far, all the considered structures have had ideal geometries, with bond angles 60 • (hex), 90 • (sq), and 120 • (hc), but now we consider also anisotropic strains. We determine the 2D where U is the energy density, x the strain in x direction, and y is the strain in y direction around a critical point of the potential energy surface. The elastic constants are obtained by calculating multiple strained structures surrounding the ideal geometry and fitting a second-order polynomial to the potential energy surface. The fit has the form where x = x − x 0 , y = y − y 0 , x 0 , y 0 , and U 0 are fit parameters, and x and y are strains in x and y directions with respect to the ideal structure (Fig. 2). The fit parameters x 0 and y 0 are included to allow for deviations from the ideal geometry in order to center on a critical point of the potential energy surface. Therefore, nonzero values for parameters x 0 and y 0 indicate that the critical point of the potential energy surface is located on a geometry with nonideal bonding angles. The energy density U used in the fit is obtained by dividing the energy of the strained structure by the area of the computational cell of the ideal equilibrium structure. The values of the resulting elastic constants are given in the Appendix and shown in Fig. 7. The general behavior of the elastic constants follows the behavior of the bulk moduli. The values of elastic constants are the highest near the middle of d series. While the elastic constants c ij are hard to determine experimentally, they can be used to calculate the Young moduli and Poisson ratios for the structures, which may be more accessible to experiments.
For the square structures the values of c 12 are large compared to c 11 and c 22 . This indicates structural instabilities against in-plane strains. To quantify this observation and to estimate the stability of a given geometry with respect to in-plane deformations, we calculate the determinant of the Hessian matrix of U : If |H| is negative, the critical point close to the ideal structure is a saddle point and the geometry is unstable. The values for |H| for studied metals and structures are shown in Fig. 8. For most metals in the hexagonal structure, |H| is positive, and structures thus stable. Many metals are stable also in the honeycomb structure. However, despite few exceptions, most metals in the square structures are unstable.

E. Bending modulus
Thus far, all deformations have been restricted to the atomic plane, but now the atoms are displaced also out of plane. We introduce deformations to calculate the bending moduli for the hexagonal and square lattices. To regain periodicity, we construct bent structures consisting of two connected cylindrical sections (Fig. 9). All the atoms are located on cylindrical surfaces, yielding constant radii of curvature throughout the entire structure, with the bond lengths fixed to the flat equilibrium values. The method ignores possible in-plane strains, caused by expansion or contraction due to bending. However, for most metals and small curvatures this contribution is small, as evidenced by fair agreement with energy density fits that are based on pure bending alone. The construction of a bent structure for a given radius of curvature R and bond length d is as follows: the angle of rotation between neighboring atoms is φ = 2 arcsin d 2R (12) and the length of the required computational cell is where N x is the number of atoms in the x direction required to make the full bend. Here, we used 14 atoms for the bent square lattices. The bending amplitude is given by  Fig. 9 is obtained. To get the bent hexagonal structures, a second row of atoms is added, with atoms rotated by φ/2 at a y distance The directions are indicated by the axes in Fig. 9. In the hexagonal structures, one bond at the connection point between the cylinders is slightly strained, but the effect is small enough to be neglected. The bending modulus κ is then obtained by fitting the energy density as moduli in 3d series have behavior qualitatively similar to the bulk moduli: the bending moduli are the highest near the middle of the series, corresponding to short bond lengths. For the 4d and 5d series, however, the behavior is different. This time, negative values for bending moduli appear near the middle of the series. Most important, the bending moduli are surprisingly large, some of them even around ∼1 eV, and thus well comparable to bending modulus in graphene (κ = 1 eV) and other covalently bound 2D materials [55,56].

F. Electron density changes mostly in plane
Finally, we consider how the 2D structure affects the electron density. In this section, we again consider ideal equilibrium structures. We begin by investigating the electron density integrated perpendicular across the atomic plane. Surprisingly, the electron density perpendicular to atomic plane remains nearly unaffected by the choice of 2D geometry. To quantify this observation, we computed the second moment of electron density n( r) for free atoms and the different 2D structures where the z direction is perpendicular to the atomic plane. These results are summarized in Fig. 11. With the exception of the lightest metals Li and Be, the second moment is almost independent of the 2D structure. However, mostly the second moment is reduced from the value of the free atoms, indicating the relocation of charge from the more diffuse free atom states to the in-plane bonding regions between the atoms. For the simple metals, the second moment reduces the most while for the transition metals it can even increase, as compared to free atoms. This is probably due to spin structures, which are more complicated in free transition-metal atoms. Let us proceed to investigate the electron density in the atomic planes. Figure 12 shows the electron density of Au for the three different structures. In all three cases, the electron density has a minimum at the corner of the Wigner-Seitz cell. In addition to the large maxima at the positions of the nuclei, there are also saddle points between the atoms in all structures. The arrangement of these critical points are rather similar in all studied metals. Further, these critical points of the electron density are located in the atomic plane.
To quantify the changes in the electron density between different metals and structures, we calculate the ratios between the densities in the local minima at the corners of the Wigner-Seitz cells and at the saddle points halfway between the atoms. This measure has been used to characterize the metallicity of a bond: the higher the ratio, the flatter the electron density in the bonding region and the more metallic the bond [57].
The resulting density ratios are larger for the simple metals than for the transition metals (Fig. 13). This is probably due to the more localized nature of the d electrons: While the s and p electrons delocalize and produce smoother electron density, the d electrons remain partly localized and produce a larger difference in the electron density in the bonding region between the atoms. An even clearer trend emerges between the structures. In general, the ratio decreases when the number of nearest neighbors decreases. The ratio is largest for hexagonal lattice, intermediate for square lattice, and smallest for the honeycomb lattice. The behavior in density ratio is related to the geometries of the Wigner-Seitz cells of the different structures. While the hexagonal and square lattices are rather densely packed, the honeycomb lattice is sparser and has large distance between nuclei and the corners of the Wigner-Seitz cell. Therefore, the electron density ratio for honeycomb is low. Large density ratios can be associated with metallic bonding and small density ratio values with covalent bonding. We note that the density ratios of square structures are intermediate between metallic and covalent bonding, which could therefore partly explain the instability of the square structures. While changing the geometry from hexagonal to square lattice decreases the density ratio and weakens the metallic bonds, it still fails to fully lead to covalent bonding as in the honeycomb structure.
Finally, we consider the density of states (DOS) of the 2D structures, using Au as an example (Fig. 14). We use Gaussian smearing of 0.1 eV and set Fermi energy to zero for each geometry. Spin-up densities are shown above zero and spin-down densities are shown below zero. Some qualitative differences are observed between the different geometries. The hexagonal structure has a large peak near −2 eV, a feature not present in the square structure. In the square geometry, the DOS is flatter and does not have clear peaks. The honeycomb structure has DOS more featured and peaked than those of the hexagonal and square geometries. In order to compare the electronic structures between different metals and structures, we calculated the DOS at Fermi energy (Fig. 15). Most metals have nonzero DOS at Fermi energy, as expected. No general trends, however, appear in these DOS values.

IV. SUMMARY AND CONCLUSIONS
We studied semi-infinite free-standing monoatomic metal layers with density-functional calculations. We found that the cohesive energies, bond lengths, and bulk moduli in the hexagonal, square, and honeycomb 2D geometries correlate linearly with the corresponding 3D bulk values. Therefore, the same properties correlate linearly also between the 2D structures. The hexagonal structures, in general, have the highest cohesion, longest bond lengths, and highest bulk moduli. Moreover, the cohesive energies are the highest, bond lengths the shortest, and 2D bulk moduli the highest near the middle of the d series. The trend in cohesion was associated with the d-band filling. The half-filled d band in the middle of the d series corresponds to filled bonding orbitals and empty antibonding orbitals. Simple metals lack the d-electron contribution to cohesion and therefore have smaller cohesive energies. Further, stronger bonds lead to shorter bond lengths and higher bulk moduli. Since cohesive energies correlate between 2D and 3D bulk structures, choosing a suitable metal for a stable 2D structure faces a tradeoff between sufficiently high 2D cohesion and sufficiently low 3D bulk cohesion. Too low 2D cohesion renders the 2D structure unstable, while too high 3D bulk cohesion leads to growth of a 3D structure.
The calculated in-plane elastic constants indicate that most metals are stable in the hexagonal and honeycomb geometries, but unstable in the square-lattice geometry. Analysis of the bending moduli showed that many metals are stable against bending, even in the square geometry. Since all the structures considered in this paper are semi-infinite, determining whether a certain geometry is stable when placed inside a 2D support structure, such as a graphene nanopore, is not directly accessible using only the values and perspectives presented here. For example, bonding to a support structure might change the stability of a given geometry due to charge transfer between support and metal. Still, the effects of strain to the total energy of a combined 2D metal in a support structure system can be estimated using the calculated elastic moduli. For instance, since graphene is extremely stiff against compression in the atomic plane, most of the deformations due to lattice mismatch will occur in the 2D metal, not in graphene. The change in energy due to this strain can be approximated using our values for semi-infinite 2D metals.
Electron density analysis showed that changing the 2D structure leaves the perpendicular density profile nearly unchanged. This was quantified by calculating the second moment of electron density in the direction perpendicular to the atomic plane. While the second moment of the electron density changes during the formation of a 2D structure from free atoms, it remains nearly constant when the 2D geometry is changed between hexagonal, square, and honeycomb lattices. Yet, in the atomic plane the electron density does change. This change was quantified by studying the critical points of the density in the plane of the atoms. The electron density is smoother for the structures with more nearest neighbors, as indicated by the larger ratios between electron densities at the minima far from the atomic nuclei and the saddle points halfway between the atoms. Decrease in this density ratio is related to the increase of covalency of the bonds. All these energetic, geometric, elastic, and electronic trends, tabulated also in the Appendix for completeness and future reference, make up an atlas that provides essential guidance in the research involving supported and free-standing 2D metal nanostructures.

ACKNOWLEDGMENT
We acknowledge the Academy of Finland for funding (Project No. 297115). Table I shows energetic, geometric, and elastic properties of hexagonal, square, and honeycomb structures.