Ab initio surface free energies of tungsten with full account of thermal excitations

The surface free energies of seven different facets of tungsten (W) are obtained up to the melting point with full account of all the relevant thermal excitations; in particular, thermal atomic vibrations, electronic excitations, and their mutual coupling. The latter is done using ab initio molecular dynamics simulations coupled with the thermodynamic integration technique. In this way, the calculations contain almost no error but the one related to the used exchange-correlation functional, which makes the results truly ﬁrst principles. The obtained results are compared with previous quasiharmonic calculations for the surface free energies of W and experimental data. The anharmonic contribution is, as expected, important for open surfaces at high temperatures, which leads to a temperature dependence of the surface energy anisotropy. The calculated Wulff shapes and surface energies are in excellent agreement with experimental data close to the melting point, where the crystalline structure of the surface layers is destroyed by a dramatic mobility of the atoms there.


I. INTRODUCTION
Experimentally, surface energies of solids can only be measured at temperatures close to the melting point, and thus high temperature data from calculations are the only way to connect calculations to experiments [1]. Such calculations are in general a nontrivial task due to the necessity to account for all the relevant thermal excitations and their coupling.
The most important contribution from the thermal lattice vibrations can be accounted for within the quasiharmonic approximation (QHA), which has been previously used in a number of ab initio surface energy calculations at finite temperature [2][3][4][5][6]. However, while the QHA is reasonably accurate for most defect-free bulk systems, it misses anharmonic contributions. These contributions can be quite substantial, especially at high temperatures in systems with open defects. For example, defect formation energies of vacancies [7,8] and stacking faults [9] show strong temperature dependencies, which are not properly captured by the QHA. A strong anharmonic contribution has also been obtained for Al surfaces in modeling with classical potentials [10].
Recently, a methodology for accurate account of the fully anharmonic contribution to the surface free energy from ab initio molecular dynamics (AIMD) has been developed and applied to the TiN(001) surface [11]. The methodology is based on the two-staged upsampled thermodynamic integration using Langevin dynamics (TU-TILD) method [8,12].
In the specific case of TiN, the commonly used QHA breaks down already for the bulk, making it necessary to use methods fully accounting for the presence of thermal lattice displacements, such as the TU-TILD method, to obtain any realistic properties, including surface free energies, at finite temperatures above 900 K. For the TiN(001) surface, however, the contribution from the surface anharmonic vibrations to the surface free energy turned out to mostly cancel that of the bulk, so that the net dependence on temperature of the TiN(001) surface free energy was relatively weak.
The latter is also partly due to the specific geometry and composition of this surface, which is nonpolar with strong bonding between the Ti and N atoms. However, it is obvious that anharmonic effects are going to be different for different surfaces, leading to different temperature dependencies of their free energies, affecting also the temperature dependence of their anisotropies.
In this paper, we study such temperature dependencies of the surface free energies and their anisotropies for W using the previously developed thermodynamic integration method for surfaces based on AIMD simulations [11]. W is a technologically important system, for which there exists experimental information for surface energies and their anisotropies at high temperature [13][14][15].
The surface free energies of W and their anisotropies have also recently been determined in ab initio QHA calculations by Scheiber et al. [6]. Although their obtained Wulff shapes are in qualitative agreement with experiments, some features are missing. The method to obtain the QHA results in Ref. [6] was to "rescale" the atomic positions with the bulk thermal expansion, without reoptimization of surface relaxations, as motivated in Ref. [3].
Such a rescaling is not needed in AIMD, which eliminates this assumption in our study in addition to the accurate account of the fully anharmonic contribution. Further, AIMD simulations allow us to take into consideration the coupling of vibrational and electronic excitations in W, which is already quite strong at temperatures half of the melting point (about 1600 K) and increases rapidly with temperature as has been demonstrated by Zhang et al. [16] in calculations for the bulk.

A. Surface free energy
The methodology used here in essence follows that used in Ref. [11]. Adjustments to the method are stated in this section, and for convenience, a short description of the method is also given in the following.
We use the slab technique to determine the surface free energy of the (hkl ) facet, γ hkl (T ), according to where F hkl (a T , T ) and F bulk (a T , T ) refer to the Helmholtz free energies of the "(hkl )" slab and bulk with the same number of atoms, for the lattice constant a T and temperature T , and A T is the surface area of the slab. The factor of 1 2 accounts for the two alike surfaces in the slab calculation. The lattice constant a T is the equilibrium one for the bulk, a eq (T ), at the given temperature and at zero pressure. For the slab geometry, a T fixes the two lateral dimensions within the surface plane. In the remaining dimension perpendicular to the surface plane, the interplane distances are in practice usually shifted for a few surface layers, but then converge to their bulk value in the middle of the slab as the slab thickness is increased. Here, a thickness of around 20 Å was estimated to be enough.
The Helmholtz free energy for bulk or slab can be adiabatically decomposed into the following contributions: where E denotes the conventional 0 K total energy of the system, F el the thermal-electronic-excitation contribution for the static lattice, and F vib the vibrational free energy of the lattice, obtained here either in the QHA or in the fully anharmonic form using the TU-TILD method. When calculated with the latter method, the vibrational free energy F vib contains the (adiabatic) coupling between one-electron excitations and lattice vibrations.
To determine an accurate thermal expansion at zero pressure required to fix a T , the full free energy surface, F bulk (V, T ), as a function of volume, V , and temperature, T , is calculated for the bulk including the same contributions as above, The thermal expansion is determined through the minima with respect to volume on this surface. In our implementation, we perform a Legendre transformation of F bulk (V, T ) to the Gibbs energy surface G bulk (P, T ), where P is the pressure, from which the equilibrium volume is obtained through the derivative with respect to pressure. Other equilibrium thermodynamic properties of the bulk, e.g., the heat capacity at constant pressure, are accessible from G bulk (P, T ) as well. As already stated, the surface free energy is given at P = 0, and all the results for the bulk are also given for P = 0.
To calculate the full vibrational free energy including the anharmonic contribution, a modified version of the original TU-TILD method [8,12] is used: where Further, F Einst is the free energy of an optimized Einstein crystal; E Einst , E MTP , and E DFT low are the energies of a particular atomic configuration calculated for the Einstein crystal, with a moment tensor potential (MTP) [17], and with low-converged density functional theory (DFT) parameters, respectively; · · · λ denotes a thermodynamic average for a particular coupling constant λ and at a certain temperature and volume/lattice constant; finally, the term E UP is obtained within free energy perturbation theory and accounts for the difference in free energy between the low-and high-converged DFT calculations.

A. Bulk
To get the bulk properties from AIMD for bulk W, we used the projector augmented-wave (PAW) method [18] as implemented in the Vienna ab initio Simulation Package (VASP) [19,20]. For the exchange-correlation energy functional we used the generalized gradient approximation in the PBEsol parametrization [21].
The choice of this functional is dictated by two facts: (1) it accurately reproduces the ground state properties of W; in particular, its equilibrium lattice constant [22] and bulk modulus, which is important for account of the vibrational contribution, and (2) it is supposed to be accurate for surfaces [21]. Thus, if the used simulation methodology does not introduce additional errors, the simulation results are expected to be a benchmark for surface energies and its anisotropies in W.
We used the W PAW potential as provided with VASP version 5.4.4, including only the 6s and 5d electrons in the valence band. (Note that several previous theoretical investigations [6,16] used different potentials with the semicore 5p electrons or 5p and 5s electrons included in the valence band.) Our tests have shown that it yields bulk properties close to the PAW potential with both the 5p and the 5s electrons in the valence band, and almost the same surface free energies due to error cancellation (at 3695 K the effect is less than 5 meV/Å 2 ).
Throughout this study, we used a plane wave energy cutoff of 300 eV for the low-converged calculations, and 500 eV for the high-converged calculations. For bulk, the integration over the Brillouin zone was performed using only the point (low converged) or using a 5 × 5 × 5 k-point grid determined according to the Monkhorst-Pack scheme [23] (high converged) for a supercell based on a 4 × 4 × 4 expansion of the cubic bcc unit cell with two atoms. This supercell contains 128 atoms. If another supecell size was used, the k-point grid was scaled correspondingly.
The self-consistent field approach [16] within finitetemperature DFT [24] was used for thermal electronic excitations, also including adiabatic coupling to vibrations during molecular dynamics (MD) runs. For obtaining the 0 K equilibrium volume, energies were fitted to a Vinet equation of state [25].
The MLIP software [17] was used to train a bulk MTP to low-converged DFT MD for a larger 6 × 6 × 6 supercell with 432 atoms for bulk at 3695 K and five different volumes spanning the thermal expansion. In these training MD runs, thermal electronic excitations were excluded and the Methfessel-Paxton scheme [26] of order one was used for integraton over the k points, with a smearing of 0.2 eV. The maximum radius used in the MTPs was 5 Å. For the bulk MTP a level 16 potential was used, with a radial basis size of eight. The level is a degreelike hyperparameter controlling the number of MTP parameters [17], which in the present work was equal to 125 for the bulk MTP. The energy root-meansquare error (RMSE) and force RMSE for the bulk MTP were 1.7 meV/atom and 0.10 eV/Å, respectively.
For the thermodynamic integration from Einstein solid to the bulk MTP, F Einst→MTP bulk , also the larger 6 × 6 × 6 supercell was used. In this way, the most important dynamic contributions are captured at the level of the bulk MTP, while the 128-atom supercell used for integration from bulk MTP to DFT, F MTP→DFT bulk , is significantly less demanding computationally, only providing a maximum error due to the smaller size of 2 meV/atom at 3695 K and a = 3.225. At lower temperatures/smaller volumes, this error is significantly smaller.
For the MD in the thermodynamic integration, Langevin dynamics [27] was used with a dampening parameter of 0.01 fs −1 . The van Gunsteren-Berendsen algorithm [28] was used for integration of the equations of motion. A time step of 2 fs was used in the thermodynamic integration, in which each MD run lasted for 50 000 time steps. In the thermodynamic integration of F Einst→MTP bulk , a dense set of λ 1 values (typically 26) was used. One thermodynamic integration was performed for each point on a mesh of five volume times 13 temperature points.
For the computationally heavier AIMD to get F MTP→DFT bulk , a less dense mesh with five temperature points was used instead. For this term, also the λ 2 values were less dense, using five λ 2 = {0, 0.15, 0.5, 0.85, 1}. To optimize statistics, two MD runs initiated with different random seeds were used at each λ 2 value, keeping the number of MD steps to a minimum (500 steps, excluding equilibration of 3000 steps with MTP).
Finally, the resulting potential energy surface was fitted to [11]  where ω(V ) is expanded in a polynomial in V up to second order and ω i (V ) polynomials in V up to third order. n was set to two.

B. Slab
We investigate five surface facets previously shown to be stable [6,29]: (100), (110), (111), (211), and (310). Among these is the (100) facet which has a known surface reconstruction [30,31] that was also considered in this work. In addition, two more facets were considered, the (210) and the (320) facets, due to their potential to be stable at higher temperatures. The slab sizes were chosen to be close to 20 Å in all three dimensions. This corresponds approximately to the dimensions of the bulk supercell (six times a). The dimensions chosen are presented in Table I.
As shown previously [11], the integration of Eq. (6) converges earlier with respect to supercell size. A test calculation showed that at 2000 K the difference was below 1 meV/atom. Therefore, for the (110) slab, Eq. (6) was performed with the smaller cell, while Eq. (5)  In principle, F hkl (a T , T ) needs to be calculated only along the thermal expansion of the bulk [a T = a eq (T )]. Previously, however, to increase numerical stability and to be able to easily account for changes in the thermal expansion (induced, e.g., by accounting for thermal contributions from the electronic system), we extended the mesh of investigated lattice constants at each temperature [11].
In this study, on the other hand, we had to optimize computer time due to the large number of facets calculated. Therefore, we calculated Eq. (6) [F hkl (a T , T ) MTP→DFT ] only along a T , as this was the most computationally expensive part. However, Eq. (5), F hkl (a T , T ) Einst→MTP was still calculated on a mesh similar to what was done for the bulk.
The latter calculations were computationally cheap, but also necessary to have on a mesh of in-plane lattice dimensions. The point is that because of the increasingly large mobility of the surface atoms above 2000 K, the integration in Eq. (5) from the Einstein solid could not be converged. Therefore, an integration over temperature had to be performed from 2000 up to 3695 K at the level of the auxiliary MTP for each slab/facet, Note that the reason for this is that the Einstein model uses reference 045403-3 3.0 0.14 positions, while the integration of Eq. (6), on the other hand, does not suffer from this problem. In these MD runs for integration over temperature the time step was 0.5 fs, and they ran for 200 000 steps, for temperatures from 2000 to 3700 K in steps of 25 K, for the same five in-plane lattice constants as used in the integration from the Einstein solid.
Further computational details specific to the contributions to F hkl (T ) are provided in the following; however, only if they differ from the respective bulk calculations; otherwise we refer to the parameter settings given for bulk in Sec. III A.
For each slab/facet, a separate MTP ["(hkl ) MTP"] was fitted to optimize the TU-TILD slab calculations. To train the (hkl ) MTPs we utilized active learning [17,32] through additional MD "training simulations." The training simulations required initial potentials, which were trained to DFT MD simulations for each (hkl ) slab, in the same way as the bulk MTP was trained. The training simulations were done for all used temperatures 2000 K and in-plane lattice constants (using different initial seeds) and over a time span of 15 ps to ensure that the trained MTPs were stable and accurate. (In total, more than 6 ns per slab/facet were run in the training simulations.) The levels of the (hkl ) MTP parameters were adjusted such as to obtain a reasonable error (lower than 5 meV/atom) for the trained MTPs. This meant that potentials of level 20 were needed for the (100) and the (110) facets, while a level of 16 was enough for the rest. With the given basis size of eight this corresponds to 329 (level 20) and 125 (level 16) parameters. The errors for MTP potentials in terms of energy and force are presented in Table II. Note, however, that the errors do not affect the accuracy of the AIMD results; rather, a smaller error indicates that the thermodynamic integration in Eq. (6) converges faster, and the computational effort is decreased.
A mesh of five in-plane lattice constants times six temperature points from 500 to 2000 K was used for the calculation of F Einst→MTP hkl . The other parameters were set as for the respective bulk calculations. The reference positions around which the Einstein contribution was calculated were chosen as the unrelaxed positions corresponding to those of ideal bulk.
In calculating Eq. (6) for the slabs, F MTP→DFT hkl , we used a 5 × 5 × 1k-point mesh for the high-converged parameters. To optimize the computational load, we used a parametrization of the temperature dependence obtained for five temperatures for Relative thermal expansion for bulk W from fully anharmonic calculations (black), the QHA (yellow dashed) and QHA plus static thermal electronic excitations (yellow). We also show results from earlier calculations with the QHA by Scheiber et al. [6], as well as experimental data from White and Minges [33] and Dubrovinsky and Saxena [34].

100
. The obtained parametrization form was then used for all the other surfaces, but adjusted for each surface using the results for three different temperatures. The adjustments had simple dependencies on temperature, and Eq. (6) could in this way be obtained with fewer calculations.
Finally, all the contributions to the slab free energies were combined at a T up to the melting temperature, T m = 3695 K, for each facet. Then, the surface free energies were calculated including all relevant contributions according to Eq. (1).

A. Bulk thermal properties
The 0 K lattice parameter obtained here was 3.145 Å. This is slightly smaller than the value obtained with the PBEsol functional by Haas et al. [22], which is practically the same as the experimental one 3.160 Å. This difference is due to the W PAW potential used here, with only six electrons in the valence shell. Nonetheless, the calculated thermal expansion is in excellent agreement with experimental data as is demonstrated in Fig. 1 where the relative thermal expansion is shown obtained in fully anharmonic AIMD simulations together with QHA results calculated here and previously [6]. For the QHA results, we also show the contribution from electronic excitations.
There is a small difference between our QHA results and those of Scheiber et al. [6], which can be due to several reasons, such as, for example, the supercell size (which the QHA can be sensitive to), the difference in the used PAW potential, or other computational parameters. We note, however, that there is also a slight difference between the effect of thermal electronic excitations ("+el") in these two calculations, The calculated temperature dependence of the adiabatic bulk modulus and isobaric heat capacity. The results are compared with experimental data for the bulk modulus from Lowrie and Gonas [35], Featherston and Neighbours [36], Qi et al. [37], and Bolef and De Klerk [38], and for the heat capacity from Cezairliyan and McClure [39] as well as from a CALPHAD assessment by Gustafson [40].
probably caused by whether its volume dependence is considered or not.
In any case, the results from AIMD are in an excellent agreement with the experimental data from Refs. [33,34] shown with red and purple crosses and pluses, respectively. While the QHA results deviate above around 2500 K, the AIMD results agree all the way up to close to the melting temperature.
The calculated heat capacity also shows an excellent agreement with experimental data as can be seen in Fig. 2. The bulk modulus is slightly higher than the experimental values, as expected from a slightly smaller lattice parameter. However, as discussed in Sec. III A, it practically does not affect the surface free energies, which deviate less than 5 meV/Å 2 at the melting point. Figure 3 shows the surface free energy of the seven different facets considered (100, 110, 111, 210, 211, 310, and 320) as functions of temperature, in the range of 500-3695 K, which is above the Debye temperature of 383 K [41,42]. The solid lines are the results of the fully anharmonic DFT calculations for each facet, respectively, including the thermal electronic contribution coupled with thermal atomic movements at the corresponding temperature. These results can be compared with the corresponding results from Ref. [6] obtained in the QHA, except for the (210) and (320) surfaces which were not considered there.

B. Surface free energies
The agreement between the fully anharmonic and QHA results is excellent at temperatures up to 2000 K, despite the difference in methodology. Above 2000 K, the results start to deviate especially for the more open surfaces, such as (111) and (310), and this deviation increases with temperature. So, around 2500 K, the surface free energies in the fully anharmonic calculations decrease significantly faster with temperature than their QHA counterparts.
It is interesting that the initial spread of the surface energies at 0 K becomes much smaller close to the melting point. The reason for that will be discussed below. The important point here, however, is that our results are in very good agreement with experiments for the liquid state assessed and performed by Paradis [14]. In particular, the surface free energy of the close packed (110) facet is practically on top of the range of experimental surface energies. Let us note that there also exists other experimental data assessed by Tyson [1]; however, these are for temperatures between 1300 and 2030 K, and the spread of this data is very large and we therefore do not include it here.

Thermal excitations' contribution
Returning to the difference between fully anharmonic and QHA results, it is clear from Fig. 3 that the difference depends on the facet. For example, for the (110) facet, the anharmonic effects and thermal-electronic-excitation coupling to the vibrations cancel those of the bulk, so that their contribution to γ 110 is very small even up above 2000 K, and it essentially follows that obtained with only the QHA and thermal electronic excitations for the static ideal lattice (from Ref. [6]). In contrast, these contributions are not small for the (100) facet, and the difference γ 100 − γ 110 is lowered by 13% even at 2000 K, as compared to the QHA from Ref. [6]. In terms of anisotropy, γ 100 /γ 110 is lowered from 1.145 for the QHA, to 1.124 in this study. Likewise, γ 310 is lowered in a similar way.
To demonstrate the origin of this different behavior, we extract adiabatically decomposed contributions to the surface free energies corresponding to Eq. (2), further dividing γ vib into representing contributions from the QHA, anharmonicities, and thermal-electronic-excitation coupling to vibrations, respectively, for different facets. With the QHA results from Ref. [6], and our calculations of γ el , and γ vib with and without thermal electronic excitations accounted for, all contributions can be determined. The results for one more open facet, (100), and for the most close packed (110) facet are shown in Table III, at two different temperatures: one intermediate, 2000 K, and at T m = 3695 K. Note, however, that the static thermal electronic contribution, just as the force constants in the QHA, depends on the reference positions used for F (a T , T ) QHA/el . Therefore, whether one uses 0 K relaxed reference positions or, for example, positions scaled with the thermal expansion, can make a considerable difference (see also, e.g., Ref. [10]). This is especially important at large a T (corresponding to high temperature) where the relaxation at 0 K is pronounced.
The γ qh values in Table III are for slabs with some layers relaxed at a T =0 K and then rescaled to a T and will be referred to as semirelaxed. For consistency, we present γ vib contributions based on semirelaxed reference positions, as well as results where γ el was based on relaxed positions for 045403-5 comparison at 3695 K, as the difference is especially large there. We will, however, in the rest of the analysis below, consider only the semirelaxed results.
From the values in Table III for (100) it can be seen that at 2000 K, γ el 100 is −3.4 meV/Å 2 . By comparing AIMD results with and without thermal electronic excitations, the contribution to γ 100 in that case is around −3.0 meV/Å 2 . The difference between these two values, i.e., γ el-vib , is negligible, and we can thereby account our lower γ 100 , compared to the QHA results, to explicit anharmonic vibrations. Likewise, γ el-vib 110 is small. So, even though the thermal-electronic-excitation coupling to vibrations have previously proven to be strong in bulk W from about half the melting temperature [16], we can see that at 2000 K the coupling effect is essentially the same at the surface, and the net γ el-vib is not large enough to influence the surface free energy.
In terms of anharmonic vibrations, γ ah 100 (2000 K) = −4.9 meV/Å 2 , which is noticeable. For the (110) facet at the same temperature, the difference from the previous QHA results is close to zero and the effect of explicit anharmonic vibrations is small in this surface, just as in the bulk.
At higher temperature, there is a significant increase in the electronic-vibrational coupling effect for the (110) surface, that increases γ 110 . The effect on γ 100 , however, is not at all that large. So, while γ ah at 3695 K is of similar magnitude for the (110) facet as for the (100) facet, γ el-vib differs and contributes to the stronger temperature dependence of γ 100 . Thus, both γ ah and γ el-vib contribute to different extents at 045403-6  different temperatures, and both are important for a correct description of the surface free energies at high temperatures.

Surface thermal disorder
In order to investigate this further and get a more clear picture of what happens with the surface structure at elevated temperatures, projected atomic trajectories from MTP MD runs (capturing the major parts of the anharmonic vibrations) at 2000 and 3000 K are obtained and plotted in Fig. 4. Although at 2000 K the overall crystal surface structure is well preserved, the thermal displacements of the atoms in the surface layers are noticeably larger than those in the bulk, especially for the (100) and (111)   above the surface layer, which also means creating vacancies in the surface layer. It is also clear that the more open surfaces become more disordered compared to the close packed (110). Obviously, at this temperature and above, the QHA should break down, and this is the origin of a strong deviation from a linear decrease of all the surface free energies at around 3000 K.
At the melting temperature, the trajectories of atoms of not only surface but several subsurface layers become diffuse. This can be clearly seen in Fig. 5, where we show atomic trajectories for the (111) and (310) facets at 3695 K. For both surfaces, the atomic displacements in several surface layers are larger than the nearest neighbor distance, which indicates an extremely high atomic mobility leading to an appearance of new "adlayers"-even two in the case of the (310) facet-and surface premelting (we do no study this effect in this work, since it requires the use of much larger supercells for a proper quantitative modeling).
Such a high degree of atomic disorder at the surface clearly shows that the crystalline character at the surface is either greatly reduced, or already broken, for the more open surfaces at high temperatures. The picture of atoms vibrating around mean positions breaks down, and a correct modeling of the surface free energy at this temperature requires going beyond the QHA. This increasingly large disorder also explains why the surface energy anisotropy decreases with temperature and that the surface free energies approach that of the structurally disordered liquid state.

C. Equilibrium crystal shape
To be able to visualize the impact of the change in surface free energy as well as compare to earlier calculations, we construct Wulff shapes based on the surface free energy at a given temperature. The resulting shapes, constructed and visualized using the WULFFPACK software [43], are presented in Fig. 6 for 500, 1000, 2000, and 3000 K. One can notice that the (210) and the (320) facets whose surface free energy was not considered earlier are present in the resulting shapes from 1000 K and up. Further, the (111) and (211) surfaces are more stable than previously suggested. This is in close agreement with the particle shapes that can be seen in Ref. [6].
At 500 K, as expected from its low surface free energy (see Fig. 3), the (110) facet occupies the most surface area. The (211) facet comes next, and we see also (310), (111), and (100). At higher T the (100) surface occupies more and more area, at the expense of (310), while (111) and (211) remains relatively constant. Between 500 and 1000 K the (320) facet appears between the (110) and (310) facets, but is at even higher temperatures, between 2000 and 3000 K, exchanged for the (210) facet. There is an interval around 2600 K where both the (210) and (320) surfaces are present.

V. SUMMARY
We have used AIMD coupled with the thermodynamic integration technique to determine surface free energies of 045403-8 1000 K seven different facets which include all the relevant thermal excitations with a high degree of accuracy, so that the remaining errors can mostly be attributed to the PBEsol exchange-correlation energy functional used in this work. The high-temperature results for the surface free energies and surface energy anisotropies are in good agreement with experimental data. The analysis of thermal-excitation contributions shows that anharmonic vibrational effects become substantial close to about 3000 K, i.e., 80% of the melting point, and they are significantly more pronounced for more open surfaces, leading to a decreasing surface energy anisotropy with temperature. At the same time, the surface energy of the most close packed (110) facet, shows a surprisingly similar dependence in our anharmonic calculations compared to earlier QHA calculations by Scheiber et al. [6], deviating only from around 2500 K and up to the melting point.
The analysis of atomic trajectories between 2000 K and the melting point indicates a dramatic increase of the mobility of atoms in the surface layers starting from 3000 K. The surface atoms can spend substantial time in adatom positions, creating corresponding vacant sites in the surface layer. This effect becomes very strong at the melting temperature, where the crystalline surface structure is almost destroyed in several surface layers.