Step free energies at faceted solid surfaces: Theory and atomistic calculations for steps on the Cu(111) surface

TL;DR

This paper develops a thermodynamic formalism for steps on faceted crystalline surfaces, applying atomistic simulations to calculate step free energies and their temperature dependence on Cu(111) surfaces.

Contribution

It introduces a new theoretical framework linking step free energy to excess quantities, enabling atomistic calculation of temperature effects on surface steps.

Findings

Step free energy remains finite up to melting point.

Weak temperature dependence of excess energy and stress below 0.6 Tm.

Significant decrease in step free energy at high temperatures.

Abstract

A theory for the thermodynamic properties of steps on faceted crystalline surfaces is presented. The formalism leads to the definition of step excess quantities, including an excess step stress that is the step analogy of surface stress. The approach is used to develop a relationship between the temperature dependence of the step free energy ($\gamma^\mathrm{st}$) and step excess quantities for energy and stress that can be readily calculated by atomistic simulations. We demonstrate the application of this formalism in thermodynamic-integration (TI) calculations of the step free energy, based on molecular-dynamics simulations, considering $\left< 110 \right>$ steps on the $\{111\}$ surface of a classical potential model for elemental Cu. In this application we employ the Frenkel-Ladd approach to compute the reference value of $\gamma^\mathrm{st}$ for the TI calculations. Calculated results for excess energy and stress show relatively weak temperature dependencies up to a homologous temperature of approximately 0.6, above which these quantities increase strongly and the step stress becomes more isotropic. From the calculated excess quantities we compute $\gamma^\mathrm{st}$ over the temperature range from zero up to the melting point ($T_\mathrm{m}$). We find that $\gamma^\mathrm{st}$ remains finite up to $T_\mathrm{m}$, indicating the absence of a roughening temperature for this $\{111\}$ surface facet, but decreases by roughly fifty percent from the zero-temperature value. The strongest temperature dependence occurs above homologous temperatures of approximately 0.6, where the step becomes configurationally disordered due to the formation of point defects and appreciable capillary fluctuations.

Equations

$$E^{\text{st}}\equiv E^{\text{st}}(S^{\text{st}},N^{\text{st}},A^{\text{st}},L,{\epsilon }_{ij}).$$

$$E^{\text{st}}=T{S}^{\text{st}}+\mu N^{\text{st}}+\gamma A^{\text{st}}+{\gamma }^{\text{st}}L,$$

$$E^{\text{t}}$$

$$E^{\text{b}}$$

$${\gamma }^{\text{st}}L=[E-TS-\mu N-\gamma A{]}_{XY}$$

$$[Z{]}_{XY}=\frac{\left|\begin{array}{ccc}{Z}^{\text{st}}& X^{\text{st}}& Y^{\text{st}}\\ Z^{\text{t}}& X^{\text{t}}& Y^{\text{t}}\\ Z^{\text{b}}& X^{\text{b}}& Y^{\text{b}}\end{array}\right|}{\left|\begin{array}{cc}{X}^{\text{t}}& Y^{\text{t}}\\ X^{\text{b}}& Y^{\text{b}}\end{array}\right|}.$$

$$[X{]}_{XY}=[Y{]}_{XY}=0.$$

$${\gamma }^{\text{st}}L=[E{]}_{AN}-T[S{]}_{AN}.$$

$$d{E}^{\text{st}}=Td{S}^{\text{st}}+\mu d{N}^{\text{st}}+\sum _{i,j}^{x,y}{\sigma }_{ij}^{\text{st}}{V}^{\text{st}}d{\epsilon }_{ij},$$

$$d(E^{\text{st}}-T{S}^{\text{st}}-\mu N^{\text{st}})=-S^{\text{st}}dT$$

$$+\sum _{i,j}^{x,y}{\sigma }_{ij}^{\text{st}}{V}^{\text{st}}d{\epsilon }_{ij}.$$

$$d({\gamma }^{\text{st}}L)=-S^{\text{st}}dT$$

$$-A^{\text{st}}d\gamma +\sum _{i,j}^{x,y}({\sigma }_{ij}^{\text{st}}{V}^{\text{st}}-{\delta }_{ij}\gamma A^{\text{st}})d{\epsilon }_{ij}.$$

$$0=-S^{\text{t}}dT$$

$$-A^{\text{t}}d\gamma +\sum _{i,j}^{x,y}({\sigma }_{ij}^{\text{t}}{V}^{\text{t}}-{\delta }_{ij}\gamma A^{\text{t}})d{\epsilon }_{ij},$$

$$0=-S^{\text{b}}dT-N^{\text{b}}d\mu +\sum _{i,j}^{x,y}{\sigma }_{ij}^{\text{b}}{V}^{\text{b}}d{\epsilon }_{ij}.$$

$$d({\gamma }^{\text{st}}L)=-[S{]}_{XY}dT$$

$$-[A{]}_{XY}d\gamma$$

$$[{\tau }_{ij}{]}_{XY}\equiv \frac{1}{L}\frac{\partial \left({\gamma }^{\text{st}}L\right)}{\partial {\epsilon }_{ij}}=\frac{1}{L}{\left[{\sigma }_{ij}V-{\delta }_{ij}\gamma A\right]}_{XY},$$

$$d{\gamma }^{\text{st}}=$$

$$-\frac{[A{]}_{XY}}{L}d\gamma +\sum _{i,j}^{x,y}([{\tau }_{ij}{]}_{XY}-{\delta }_{ix}{\delta }_{jx}{\gamma }^{\text{st}})d{\epsilon }_{ij},$$

$$[{\tau }_{ij}{]}_{XY}={\delta }_{ix}{\delta }_{jx}{\gamma }^{\text{st}}+\frac{\partial {\gamma }^{\text{st}}}{\partial {\epsilon }_{ij}}.$$

$$d{\gamma }^{\text{st}}=-\frac{[S{]}_{AN}}{L}dT+\sum _{i,j}^{x,y}([{\tau }_{ij}{]}_{AN}-{\delta }_{ix}{\delta }_{jx}{\gamma }^{\text{st}})d{\epsilon }_{ij}.$$

$$[{\tau }_{ij}{]}_{AN}=\frac{1}{L}{\left[{\sigma }_{ij}V\right]}_{AN}.$$

$$d\left(\frac{{\gamma }^{\text{st}}L}{T}\right)=$$

$$-\frac{[A{]}_{XY}}{T}d\gamma +\sum _{i,j}^{x,y}\frac{[{\tau }_{ij}{]}_{XY}L}{T}d{\epsilon }_{ij},$$

$$\operatorname{d}\!{\left(\frac{\gamma^{\text{st}}L}{T}\right)}=-\Bigg{(}\frac{[E]_{AN}}{T^{2}}-\sum_{i,j}^{x,y}\frac{[\tau_{ij}]_{AN}L}{T}\frac{\operatorname{d}\!\varepsilon_{ij}}{\operatorname{d}\!T}\Bigg{)}\operatorname{d}\!{T},$$

$$\alpha_{ij}\equiv\left(\frac{\partial\varepsilon_{ij}}{\partial T}\right)_{\mbox{\boldmath$\sigma$}^{\text{b}}=0}$$

$$\operatorname{d}\!{\left(\frac{\gamma^{\text{st}}L}{T}\right)}=-\Bigg{(}\frac{[E]_{AN}}{T^{2}}-\frac{2\alpha\,[\tau_{\text{avg}}]_{AN}L}{T}\Bigg{)}\operatorname{d}\!{T}$$

$$[{\tau }_{\text{avg}}{]}_{AN}=\frac{[{\tau }_{xx}{]}_{AN}+[{\tau }_{yy}{]}_{AN}}{2}$$

$${\gamma }^{\text{st}}(T)=$$

$$\displaystyle-\frac{T}{L(T)}\int_{T_{0}}^{T}\Bigg{(}\frac{[E]_{AN}}{T^{\prime 2}}-\frac{2\alpha\,[\tau_{\text{avg}}]_{AN}L}{T^{\prime}}\Bigg{)}\operatorname{d}\!{T^{\prime}}.$$

$$H_0=\sum _{i=1}^N\frac{{\mathbf{p}}_i^2}{2m}+U({\mathbf{r}}_1,{\mathbf{r}}_2,...,{\mathbf{r}}_N)$$

$$H_{\text{E}}=\sum _{i=1}^N\frac{{\mathbf{p}}_i^2}{2m}+\sum _{i=1}^N\frac{1}{2}k{\left({\mathbf{r}}_i-{\mathbf{r}}_i^0\right)}^2$$

Full text

I Introduction

In theories of crystal morphologies and growth kinetics, a property of fundamental importance is the step free energy, $\gamma^{\text{st}}$, i.e., the excess free energy of a step line defect on a faceted solid-liquid interface or crystal surface Williams (1994); Jeong and Williams (1999). The magnitude of $\gamma^{\text{st}}$ controls the depth of the cusp in the interfacial free energy versus orientation plot for faceted interfaces, and this property is thus fundamental in determining the equilibrium crystal shape Einstein (2014). The step free energy also plays an important role in governing crystallization kinetics from the melt or vapor Teng et al. (1998), by controlling the magnitude of the barrier to island nucleation in the growth of a faceted surface or interface. The step free energy can depend strongly on temperature, and this dependence ultimately leads to the vanishing of $\gamma^{\text{st}}$ above the thermodynamic roughening temperature Williams et al. (1993); *[Erratum:]faceting_erratum. Below the roughening transition, steps with free energies that are low relative to the thermal energy will display pronounced capillary fluctuations, which have important consequences for their kinetic properties, bunching instabilities Tersoff et al. (1995), and morphologies Bartelt et al. (1992, 1990).

Despite the importance of $\gamma^{\text{st}}$ described above, measurements of this quantity remain relatively rare. Further, reported values are often available only for a fixed value of the temperature Schlößer et al. (1999); Giesen et al. (2001); Steimer et al. (2001); Gabrisch et al. (2001) and measurements over a wide temperature range have been undertaken in few systems Schulze Icking-Konert et al. (1999); Zandvliet (2003). As a consequence, knowledge of the nature of the temperature dependence of step free energies and understanding of the microscopic factors that underlie it remain incomplete. This situation presents a challenge for the development and application of quantitative mesoscale theories in studies of faceted crystal growth phenomena in real systems, and robust methods for the direct calculation of temperature-dependent step free energies from atomic-scale simulations are thus of fundamental interest. In the present paper we present a thermodynamic formalism that relates the temperature dependence of $\gamma^{\text{st}}$ on faceted crystal surfaces to excess quantities that can be computed directly by atomistic simulations. This formalism provides a framework for the calculation of $\gamma^{\text{st}}$ as a function of temperature through the thermodynamic integration of an appropriate adsorption equation.

The remainder of this paper is organized as follows. In Sec. II the thermodynamic formalism is introduced, and the relevant step excess quantities and other fundamental thermodynamic equations are defined. In Sec. III we demonstrate how the thermodynamic formalism can be combined with the calculation of a reference step free energy at low temperatures by the Frenkel-Ladd method Frenkel and Ladd (1984); Freitas et al. (2016), to compute $\gamma^{\text{st}}$ up to high temperatures, accounting for contributions arising from the formation of surface point defects and capillary fluctuations. In Secs. IV and V we present simulation details and results, respectively, of an application of the equations derived to the calculation of the free energy of $\langle 110\rangle$ steps present on the $\{111\}$ surface of face-centered-cubic copper using molecular dynamics simulations. In Sec. V we also compare the step free energy obtained here to experimentally measured Steimer et al. (2001) and first-principles-calculated Feibelman (1999); Da Silva et al. (2006); Stasevich et al. (2006) results available in the literature. Finally, in Sec. VI we summarize the main findings, and discuss applications of the formalism presented in this work more generally.

II Thermodynamic theory of surface steps

II.1 Step excess quantities

Consider a thermodynamic system that consists of a homogeneous solid with a stepped surface, where the step separates two flat surface terraces as shown in Fig. 1. Both terraces have the same structure and thermodynamic properties, while the surface around the step has properties that are different from those of the flat terraces. The step region, terraces, and bulk are in thermodynamic equilibrium with each other. We assume that atoms can migrate by diffusion between the bulk and the surface regions, allowing the concentration of point defects to vary everywhere in the system, in a way required to maintain equilibrium. We also assume that atoms can attach and detach from the step, so the system is in equilibrium with an infinite source and sink of atoms.

The properties of the bulk far away from the surface region and the properties of the flat terraces far away from the step region are well described by standard bulk and interfacial thermodynamic relations Gibbs (1906); Cahn (1979). In this section we will address the thermodynamic properties of the steps on the crystalline surface. Consider an imaginary region that contains a finite segment of the step, as shown in Fig. 1. The lower boundary of this region is located inside the homogeneous part of the bulk crystal, while the side boundaries parallel and perpendicular to the step line cross the system surface normal to the terraces. The latter condition is important because it defines the total surface area inside the region.

The extensive thermodynamic properties of the step region depend on its dimensions, since the enclosed system is not homogeneous. We postulate that the total energy of the region is a function of the following extensive and intensive variables:

[TABLE]

The superscript “st” refers to variables of the step region shown in Fig. 1: $S^{\text{st}}$ is the entropy, $N^{\text{st}}$ is the number of atoms, $A^{\text{st}}$ is the surface area enclosed by the region, $L$ is the step length, and $\varepsilon_{ij}$ are the lateral components of strain, with $i=x,y$ and $j=x,y$. The lateral components of strain in Eq. (1) correspond to the macroscopic strain in the homogeneous bulk lattice far away from the step, not to be confused with the local inhomogeneous strain around the step. As illustrated in Fig. 1, the coordinate system is chosen such that $\mathbf{\hat{z}}$ is normal to the terraces, while $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ are parallel and normal to the step line, respectively.

Consider a variation when the physical state of the system is fixed and we extend the boundaries of the region from zero to some finite values. Assuming that $E^{\text{st}}$ is a homogeneous function of degree one with respect to $S^{\text{st}}$, $N^{\text{st}}$, $A^{\text{st}}$, and $L$ we obtain

[TABLE]

where $T$ is the temperature, $\mu$ is the chemical potential, $\gamma$ is the surface free energy per unit area, and $\gamma^{\text{st}}$ is the step free energy per unit step length. Note that by definition $\gamma$ is the property of the terrace uninfluenced by the surface step. We assume that the external pressure is zero since the solid is in contact with vacuum.

At this point a comment should be made about the meaning of the quantities introduced above. The current thermodynamic treatment is focused on steps on solid surfaces, and it is well known that such steps produce long-range elastic fields Marchenko and Parshin (1980); Shilkrot and Srolovitz (1996); Müller and Saúl (2004). As a result, both terraces and the bulk crystal are strictly speaking inhomogeneous in the entire system. Equation (2) can still be used to describe the system if $\gamma$ and $\gamma^{\text{st}}$ are understood as the properties of the terraces and the step in the limit when the system size goes to infinity. In other words, even though the inhomogeneity due to the strain fields induced by the step can extend far away from the step line, its total contribution to the energy of the system is finite. This will be demonstrated using atomistic simulations in Sec. IV of this study. This property of surface steps should be contrasted with the case of lattice dislocations, which are also line defects. Different from steps, the elastic contribution to the total energy of a dislocation diverges with the system size Hirth and Lothe (1982) and Eq. (2) would not apply.

The amount of bulk and terrace inside the step region at this point is arbitrary, and hence quantities in Eq. (2) depend on the choice of the step region. In order to define the step excess quantities we need to subtract the bulk and terrace contributions from the quantities of the step region in Fig. 1. To this end we write equations analogous to Eq. (2) for the terrace and bulk regions shown in Fig. 1:

[TABLE]

where superscripts “t” and “b” refer to terrace and bulk respectively. These two regions are located sufficiently far away from the step that their extensive properties are not affected by it. The terrace region includes the surface as well as a portion of the homogeneous bulk phase, while the bulk region is unaffected by the surface. Solving the system of equations given by Eqs. (2), (3), and (4) using Cramer’s rule, we obtain an expression for step free energy $\gamma^{\text{st}}$:

[TABLE]

where $X$ and $Y$ are any of the extensive quantities $S$, $N$ or $A$. Terms $[Z]_{XY}$ are Cahn’s determinants and are calculated as the ratio of two determinants Cahn (1979)

[TABLE]

The first row of the numerator contains extensive thermodynamic quantities of the region containing the step, while the second and the third rows contain properties of regions enclosing the terrace and the bulk, respectively. According to properties of determinants if any two columns are equal, the determinant is zero:

[TABLE]

Thus, two terms in Eq. (5) automatically vanish.

The quantity $[Z]_{XY}$ has the meaning of the excess property of a step when the region with the step has the same amount of $X$ and $Y$ as the terrace and the bulk regions combined. The excess quantities generally depend on the choice of the extensive variables $X$ and $Y$. On the other hand, Eq. (5) shows that all different choices of $X$ and $Y$ result in the same excess amount of $\gamma^{\text{st}}L$.

Considering a particular example when $X=A$ and $Y=N$, we obtain

[TABLE]

The step excess quantities $[E]_{AN}$ and $[S]_{AN}$ are the excess energy and entropy when the step region has the same surface area and number of atoms as terrace and bulk regions combined. The choice of $A$ as one of the extensive variables means that the excess area of a step is zero. The step is represented as a dividing line on the surface and the properties of terraces are extended all the way to this line. Using this representation, step excess quantities can be formulated in a manner similar to the Gibbs dividing surface construction for interfaces. On the other hand, the derivation that uses Cahn’s determinants provides expressions for excess quantities that are more general. The ability to choose different definitions can be useful in applications because some excess quantities are more accessible than others to measurements or calculations. In Sec. IV we describe how several step excess quantities can be calculated directly from atomistic simulations by making use of the flexibility provided by Cahn’s determinants.

II.2 Adsorption equation

In the previous section we derived an expression for the step free energy and other excess quantities. We are now in a position to derive an equation that describes how $\gamma^{\text{st}}$ changes with temperature and mechanical deformation, namely the adsorption equation. Consider a variation of state when the system exchanges heat and does mechanical work. For the region containing the step the change in total energy is given by

[TABLE]

where $\sigma^{\text{st}}_{ij}$ is the stress tensor and $V^{\text{st}}$ is the volume. The product $\sigma^{\text{st}}_{ij}V^{\text{st}}$ is defined as the derivative of $E^{\text{st}}$ with respect to elastic deformation $\varepsilon_{ij}$. Equation (9) assumes that the surface area changes due to elastic deformation of the lattice and not by incorporation of new lattice units. At the same time the number of atoms in the region and the relative areas of the terraces can change by diffusion and attachment of atoms to the step. The conditions for mechanical equilibrium between the system and the vacuum Landau and Lifshitz (1986); Müller and Saúl (2004) require that $\sigma^{\text{st}}_{iz}=0$ for $i=x,\,y,\,\text{or}\,z$. Thus, all summations involving the stress tensor are over the $x$ and $y$ indices only.

Performing a Legendre transformation on terms containing entropy and number of particles we obtain from Eq. (9)

[TABLE]

Combining Eqs. (2) and (10) we obtain

[TABLE]

The intensive variables on the right-hand side in Eq. (11) are not independent since equations similar to Eq. (11) for the terrace and bulk regions impose additional constrains. For the terrace we have Frolov and Mishin (2009a)

[TABLE]

while the Gibbs-Duhem equation for the bulk reads

[TABLE]

Solving Eqs. (11), (12), and (13) using Cramer’s rule Cahn (1979), we obtain the adsorption equation for steps

[TABLE]

where $X$ and $Y$ are any of the extensive quantities $S$, $N$, $A$, or $(\sigma_{ij}V-\delta_{ij}\gamma A)$. Notice that the coefficients of the differentials in Eq. (14) are the step excess quantities introduced earlier in Eq. (6), and are independent of the particular choice of the regions illustrated in Fig. 1. Due to the property of determinants in Eq. (7), two terms in the adsorption equation can be eliminated by specifying $X$ and $Y$, leaving only independent variables. The number of variables should coincide with the number of degrees of freedom available to the system. Consider the same example given in Sec. II.1, where we choose $X$ and $Y$ equal to $A$ and $N$. In this case the four possible variations are changes in temperature and deformation described by strains $\varepsilon_{xx}$, $\varepsilon_{yy}$, and $\varepsilon_{xy}$. It is natural to have the step free energy be a function of these variables. The differential of surface free energy $\gamma$ that appears in Eq. (14) is an unusual variable to describe the changes in the thermodynamic state of the step. While such an exotic form of the adsorption equation can be formulated and is consistent with the Gibbs phase rule, in most practical cases it is more convenient to eliminate this term by specifying $X=A$.

II.3 Step stress

Equation (14) introduces a new excess property in addition to the quantities that appeared in Eq. (5). The last term in Eq. (14) describes changes in $\gamma^{\text{st}}L$ due to elastic deformations and defines the step excess stress as

[TABLE]

where $i=x,y$ and $j=x,y$. $[\tau_{ij}]_{XY}$ is a quantity with units of energy per length that represents the additional force exerted on the perimeter of the stepped surface due to the presence of the step. Different from $\gamma^{\text{st}}$, the step excess stress $[\tau_{ij}]_{XY}$ is not a unique quantity: it is a direct consequence of the derived adsorption equation that one can introduce several valid step excess stresses by specifying different extensive properties $X$ and $Y$. Notice that by the derivation above [\mbox{\boldmath\tau}] is a second rank tensor, not a scalar like step free energy; hence it has nonzero components parallel and normal to the step line Li et al. (2011).

Differentiating the product $\gamma^{\text{st}}L$ in Eq. (14) and using $\operatorname{d}\!{L}=\sum_{i,j}^{x,y}\delta_{ix}\delta_{jx}L\operatorname{d}\!{\varepsilon_{ij}}$ we obtain the intensive form of the adsorption equation:

[TABLE]

where the differential coefficients are the step excess quantities per unit step length. From Eq. (16) we can now obtain the relation between $[\tau_{ij}]_{XY}$ and $\gamma^{\text{st}}$:

[TABLE]

Equations (15) and (17) are the step analogs of the stress equations for solid surfaces Shuttleworth (1950); Frolov and Mishin (2009a, 2012a). They are a direct consequence of the derived adsorption equation, Eq. (14), and give a recipe for how $[\tau_{ij}]_{XY}$ can be calculated as an excess property using the determinant formalism.

Consider the example discussed earlier (Secs. II.1 and II.2) when $X=A$ and $Y=N$. This choice of extensive variables eliminates differentials of chemical potential $\mu$ and surface free energy $\gamma$, leaving only independent variations with temperature and deformation:

[TABLE]

The second term in Eq. (18) describes how the step free energy changes when the surface is deformed at constant temperature. Notice that during such a process the chemical potential $\mu$ and free energy of the terraces $\gamma$ are not constant. Equation (18) defines a particular step excess stress given by

[TABLE]

The components of this stress tensor have been calculated in the present work from atomistic simulations, and the magnitudes of this quantity will be presented in Sec. V below.

III Thermodynamic integration formalism

In this section we describe how the equations derived in Sec. II provide a framework for a thermodynamic-integration approach to computing the temperature dependence of the step free energy $\gamma^{\text{st}}$ by atomistic simulations. We also demonstrate how the absolute free energy of the step can be derived at low temperatures (i.e., where the concentration of kinks and surface adatoms are sufficiently low that we can neglect their contribution to the free energy) using the Frenkel-Ladd Frenkel and Ladd (1984) method, to provide a reference value in the thermodynamic integration approach. The combination of these two methods provides a general framework for the calculation of step free energies over a wide temperature range, accounting naturally for vibrational and configurational disorder.

III.1 Gibbs-Helmholtz relation for step free-energy integration

The temperature dependence of the step free energy can be obtained by directly integrating $\operatorname{d}\!{(\gamma^{\text{st}}L)}$, given in Eq. (14), along a reversible thermodynamic trajectory. However, in many applications the calculation of the excess entropy $[S]_{XY}$ can be challenging. Fortunately, we can avoid the explicit calculation of $[S]_{XY}$ by integrating $\operatorname{d}\!{(\gamma^{\text{st}}L/T)}$ instead of Eq. (14). We can compute $\operatorname{d}\!{(\gamma^{\text{st}}L/T)}$ explicitly by combining Eq. (5) and (14):

[TABLE]

where $[\tau_{ij}]_{XY}$ depends on the choice of the $X$ and $Y$ variables. Equation (20) is the surface step analog of the Gibbs-Helmholtz equation from bulk thermodynamics. A similar equation for interfaces was derived previously Frolov and Mishin (2009b), and was demonstrated to be efficient for calculating the temperature dependence of interface free energies Frolov and Mishin (2009a, b, 2010a, 2012b); Laird et al. (2009); Laird and Davidchack (2010).

Before integrating Eq. (20) we need to choose $X$ and $Y$ since the selection of these variables determines which quantities need to be calculated to perform the thermodynamic integration. A convenient choice for the applications considered here is $X=A$ and $Y=N$. In this case Eq. (20) becomes

[TABLE]

where $[\tau_{ij}]_{AN}$ is given by Eq. (19). Equation (21) can be integrated along a reversible thermodynamic path, where the temperature is increased from $T_{0}$ to $T$ while the solid is expanded to accommodate the thermal expansion, effectively maintaining zero bulk stress, i.e., \mbox{\boldmath\sigma}^{\text{b}}=0. Notice that this thermodynamic path couples the a priori independent variables $T$ and $\varepsilon$:

[TABLE]

where $\alpha_{ij}$ is the linear thermal-expansion factor and $i$ and $j$ are equal to $x,y,\text{ or }z$. We will assume here that the crystal lattice has cubic symmetry, allowing us to define our coordinate system in a way that eliminates the dependence of $\alpha_{ij}$ on the indexes $i$ and $j$. One further implication of following this thermodynamic path is that the system is not subject to shear strain during the thermal expansion; hence $[\tau_{ij}]_{AN}$ for $i\neq j$ performs no mechanical work. With these considerations Eq. (21) becomes

[TABLE]

where

[TABLE]

is the average step stress. Upon integration of Eq. (22) following the thermodynamic path described above we obtain

[TABLE]

Note that all quantities inside the integral depend on the temperature $T^{\prime}$.

Equation (24) allows for the calculation of the temperature dependence of $\gamma^{\text{st}}$ if we know how to calculate all quantities on its right-hand side. The excess quantities inside the integral on the right-hand side of Eq. (24), $[E]_{AN}$ and $[\tau_{\text{avg}}]_{AN}$, can be computed readily from atomistic simulations since they only involve the calculation of the system energy and stress tensor. Thus, the only remaining term on the right-hand side of Eq. (24) is the step excess free energy at a reference temperature $\gamma^{\text{st}}(T_{0})$. This type of term, present in all thermodynamic integration methods, cannot be trivially computed using atomistic simulations, since it involves the calculation of the absolute free energy of the system. In the next section we present a method due to Frenkel and Ladd Frenkel and Ladd (1984); Freitas et al. (2016) which enables the calculation of the absolute free energy of solid systems. In the present context, this method enables the calculation of $\gamma^{\text{st}}(T_{0})$ provided the temperature $T_{0}$ is chosen low enough such that the steps are structurally ordered (i.e., without an appreciable concentration of kinks, adatoms, or vacancies). Once we know the free energy of the step at this reference temperature, we can use Eq. (24) to compute the absolute free energy of the step at any other temperature $T$ from values of $[E]_{AN}$ and $[\tau_{\text{avg}}]_{AN}$ at temperatures between $T_{0}$ and $T$.

III.2 Application of Frenkel-Ladd approach for calculation of step free energies

The Frenkel-Ladd Frenkel and Ladd (1984) (FL) method is a type of thermodynamic integration approach that allows calculation of the absolute free energy of crystalline solids from atomistic simulations. Consider a system composed of $N$ identical particles with the Hamiltonian

[TABLE]

where $m$ is the mass of the particles and $U(\mathbf{r}_{1},\mathbf{r}_{2},...,\mathbf{r}_{N})$ is a many-body interatomic potential. We assume that, at the temperature and pressure of interest, the system’s stable phase is a solid with a known crystalline lattice structure. Considering this lattice structure we will construct a second Hamiltonian for a reference Einstein crystal, which consists of particles of the same mass $m$ attached to the equilibrium lattice sites by harmonic springs with spring constant $k$:

[TABLE]

where $\mathbf{r}_{i}^{0}$ is the equilibrium lattice position of particle $i$ in the system described by $H_{0}$.

In the FL method we use a Hamiltonian which is a linear interpolation of the Hamiltonians given by Eqs. (25) and (26):

[TABLE]

where $\lambda$ is a parameter of this Hamiltonian. The free energy of the system $H(\lambda)$ is

[TABLE]

where $k_{\text{B}}$ is the Boltzmann constant, $h$ is the Planck constant, $\mathbf{x}=\{\mathbf{r}_{1},\mathbf{r}_{2},...,\mathbf{r}_{N},\mathbf{p}_{1},\mathbf{p}_{2},...,\mathbf{p}_{N}\}$ is a point in the phase space of the particles of this system, and $\beta=1/k_{\text{B}}T$. It can be easily shown, by computing the derivative of Eq. (28), that

[TABLE]

where $\langle\ldots\rangle_{\lambda}$ is the canonical ensemble average for a specific value of the parameter $\lambda$. From direct integration of the equation above from $\lambda=0$ to $\lambda=1$ we obtain

[TABLE]

where $F_{0}(N,V,T)\equiv F(N,V,T;\lambda=0)$ is the free energy of the solid described by $H_{0}$, $F_{\text{E}}(N,V,T)\equiv F(N,V,T;\lambda=1)$ is the free energy of the Einstein crystal, and $U_{\mathrm{E}}$ is the potential energy of the harmonic springs in the Einstein crystal. Since $H_{\text{E}}$ is composed of independent harmonic oscillators we can calculate its free energy analytically:

[TABLE]

where $\omega=\sqrt{k/m}$ is the natural frequency of the harmonic oscillators.

Equations (29) and (30) allow calculation of the absolute free energy of the solid $H_{0}$ from atomistic simulations. The only unknown in Eq. (29) is the integrand on the right-hand side, which is an equilibrium ensemble average and, therefore, can be calculated directly using atomistic simulation techniques Frenkel and Smit (2001) such as molecular dynamics or Monte Carlo with the Hamiltonian given by Eq. (27). The evaluation of Eq. (29) can be performed in a straightforward manner using equilibrium simulations to obtain averages necessary to calculate the integral on the right-hand side numerically. However, this is an inefficient way to perform this thermodynamic integration. State-of-the-art methods Freitas et al. (2016) for evaluating Eq. (29) based on nonequilibrium simulations have been developed and are now implemented in high-performance atomistic simulation software such as LAMMPS Plimpton (1995) (Large-scale Atomic/Molecular Massively Parallel Simulator). These methods drastically reduce the computational cost of the thermodynamic-integration calculation and provide robust error-control criteria. An in-depth description of these techniques and detailed account of how they can be implemented in practice is given in Ref. Freitas et al., 2016.

We have shown in Sec. II, Eqs. (5) and (8), that the step free energy $\gamma^{\text{st}}(T_{0})$ is a quantity that can be computed from the free energies of the three different regions shown in Fig. 1. Hence, our approach in the work presented below is to obtain the $\gamma^{\text{st}}(T_{0})$ using the FL method to compute the absolute free energies of the relevant required systems. In so doing we have followed closely the methodology described in Ref. Freitas et al., 2016 to perform the FL calculations. Note, however, that the FL method has its applicability limited to low-temperature surfaces (flat and stepped), since at high temperatures the presence of surface vacancies, adatoms, and kinks on the steps breaks the FL method assumption that the atomic motion occurs around the equilibrium lattice positions, Eq. (26). Thus, the free energy computed with the FL method is only used as an initial integration point for the thermodynamic integral approach of Sec. III.1, more specifically in Eq. (24).

IV Atomistic simulations

IV.1 Methodology

To demonstrate the application of the methodology described in the previous section, for computing step free energies by atomistic simulations, we focus on the $(111)$ surface of face-centered-cubic Cu, modeled with the embedded-atom-method (EAM) interatomic potential due to Mishin et al. (2001). In previous simulations it has been found that this surface remains faceted at all temperatures up to the melting point of the EAM model ($T_{\text{m}}=1327\,\text{K}$ for the potential model considered Frolov and Mishin (2010b)). No evidence for surface premelting was observed in these previous simulations, such that the surface maintains the layered crystalline structure up to $T_{\text{m}}$. Since the surface remains faceted, the step free energies are expected to remain finite up to this temperature.

We have chosen molecular dynamics (MD) as the atomistic simulation technique to evaluate the step excess quantities necessary for the thermodynamic integration equations. All calculations were performed using LAMMPS Plimpton (1995), an open source implementation of MD. The Langevin thermostat Schneider and Stoll (1978) was employed to sample the phase space, according to the canonical ensemble distribution. The relaxation time used for the thermostat was $\tau_{\text{L}}\equiv m/\gamma=20\text{ps}$, where $\gamma$ is the friction parameter and $m$ is the atomic mass. The timestep was chosen based on the highest-frequency normal mode of the system ($\nu_{\text{max}}=7.8\times 10^{12}\text{ Hz}$); we have taken $\Delta t$ to be approximately $1/60$th of the oscillation period of that normal mode: $\Delta t=2\,\text{fs}$.

IV.2 System geometry and dimensions

In Sec. III Eq. (24) was derived for the temperature dependence of $\gamma^{\text{st}}$, based on the choice $X=A$ and $Y=N$. In this sub-section we elaborate further why this is a convenient choice for the calculation of the excess quantities that appear in Eq. (24) from atomistic simulations. From Eq. (6) we have

[TABLE]

The bulk region in Fig. 1 does not have any surface which means $A^{\text{b}}=0$. Heretofore the regions shown in Fig. 1 had arbitrary dimensions; from now on we choose the dimensions of the step and terrace regions in such a way that they have the same surface area, i.e., $A^{\text{st}}=A^{\text{t}}$. Furthermore, we choose the depth of these regions such that they contain the same number of atoms: $N^{\text{st}}=N^{\text{t}}$. With this particular choice of dimensions, the excess quantities shown in Eq. (31) become $[Z]_{AN}=Z^{\text{st}}-Z^{\text{t}}$. Thus, the need to compute the thermodynamic properties for the bulk ($Z^{\text{b}}$) is minimize and $[Z]_{AN}$ becomes a simple difference between the properties of the step and terrace regions.

To calculate step excess quantities we modeled two different simulation blocks illustrated in Fig. 2. The simulation block shown in Fig. 2a is a solid film with two flat $(111)$ surfaces. Periodic boundary conditions were applied for the directions parallel to the surface. The second simulation block illustrated in Fig. 2b was obtained from the first one by adding half of an atomic plane on the top surface and removing half of the atomic plane from the bottom surface. As a result of this construction, the second block has four surface steps. At the same time the construction ensures that the two simulation blocks have the same number of atoms and the same surface area. Properties $Z^{\text{st}}$ and $Z^{\text{t}}$ were then calculated for the two blocks with and without steps, respectively. The difference between these quantities gives the step excess $[Z]_{AN}$ given by Eq. (31). Indeed, $Z^{\text{st}}-Z^{\text{t}}$ represents the excess of property $Z$ due to steps, when the reference system has the same surface area and the same number of atoms. We remind the reader that the excess quantities inside the integral on the right-hand side of Eq. (24) are $[E]_{AN}=E^{\text{st}}-E^{\text{t}}$ and $[\tau_{ii}]_{AN}=(\sigma_{ii}^{\text{st}}V^{\text{st}}-\sigma_{ii}^{\text{t}}V^{\text{t}})/L$ and can be readily computed from atomistic simulations since they involve only the calculation of the energy and stress tensor of each of the systems in Fig. 2.

The steps considered in the MD simulations of the simulation cells illustrated by Fig. 2b are directed along the close-packed $\langle 110\rangle$ direction. The crystallographic symmetry of the $(111)$ surface is such that the two steps shown in each of the surfaces of Fig. 2b are slightly different; it can be seen in Fig. 3 that they have different nearest-neighbor configurations on the $(111)$ plane immediately below the surface. The step with lowest zero-temperature energy Müller and Saúl (2004) ($U_{0}=103.13\,\text{meV}/\,\mathring{\mathrm{A}}$) is a $\langle 110\rangle$A step, while the step with the slightly larger energy ($U_{0}=104.08\,\text{meV}/\,\mathring{\mathrm{A}}$) is a $\langle 110\rangle$B step. Using the terminology of Ref. Feibelman, 1999 $\langle 110\rangle$A steps have $\langle 100\rangle$ microfacets and $\langle 110\rangle$B steps have $\langle 111\rangle$ microfacets.

We have chosen the simulation box size in such a way that the step-step interaction energy of all four steps in Fig. 2b was negligible compared to the step self-energy (i.e., the energy of an isolated step). Steps are abrupt interruptions of the surface first layer; hence, they deform the atomic structure around them, creating an elastic field Marchenko and Parshin (1980). The interaction energy due to the overlap of the strain fields of the two steps decays as $d^{-2}$ with the step-step separation and exponentially with the bulk depth (see Ref. Shilkrot and Srolovitz, 1996 and references therein). Following the work of Shilkrot and Srolovitz (1996) we have verified this behavior for the step elastic interaction energy 111If these interactions are appreciable we would expect a significant difference between results obtained from supercells with symmetric (as shown in Fig. 2) versus antisymmetric indentations on the top and bottom surfaces. But if the crystal slab is thick enough to minimize step interactions the differences would be negligible. We performed a separate independent test of the effect of system size, by comparing the step self energy derived from supercells with the geometry shown in Fig. 2, with results obtained using non-orthogonal boundary conditions Shilkrot and Srolovitz (1996) where there is only one step on each surface. The values for the step self energy agreed for both geometries to within $0.01\%$. of our model and we have determined the step-step distance ($d$) and bulk depth ($h$) such that $E_{\text{int}}/U_{0}\leq 10^{-4}$, where $E_{\text{int}}$ is the total step interaction energy and $U_{0}$ is the step self-energy at zero temperature. The box dimensions obtained are $d=70.8\,\mathring{\mathrm{A}}$ and $h=53.2\,\mathring{\mathrm{A}}$ at $T=0\,\text{K}$; for finite temperatures we have increased the system dimensions to account for thermal expansion, for zero bulk stress.

The simulation box length along the step line cannot be determined based on static simulations. In order to determine the step length necessary to eliminate finite-size effects along the step direction it is necessary to consider fluctuations of the step line that appear at finite temperatures, known as capillary fluctuations. The accurate evaluation of step excess quantities requires satisfactorily sampling the normal modes of these fluctuations (i.e., the capillary waves) during the simulation. If the step length used is too small the sampling of long-wavelength modes is suppressed. On the other hand, an excessively lengthy step would make the thermodynamic integration calculations prohibitively long due to the need to sample normal modes with very long wavelengths and associated long relaxation times. Thus, to determine the step length required in the simulations we need to analyze the convergence of $[E]_{AN}$ and $[\tau_{\text{avg}}]_{AN}$ with the step length since, according to Eq. (24), the thermodynamic integration equation depends on the computation of these two quantities.

Using the values of $d$ and $h$ determined above we have run simulations at $T=1300\,\text{K}$ for systems with different step lengths ($L$) and calculated the step excess energy and stress. Figure 4 shows the convergence of the step excess quantities with step length for these simulations. Based on these results we have chosen $L=30.7\,\mathring{\mathrm{A}}$ as the step length for the next simulations since the step excess quantities are seen to be well converged for steps of this size.

V Results and discussion

V.1 Step free energies from Frenkel-Ladd simulations

Using the box dimensions specified in Sec. IV.2, we have constructed two systems to be used in the FL simulations, one with a flat surface and the another with stepped surfaces, as shown in Figs. 2a and b. Both systems have $39168$ atoms and the same surface area. In Fig. 5 we show plan-view snapshots Stukowski (2010) of the top layer of atoms from typical configurations for the stepped system at different temperatures. At temperatures of $700\,\text{K}$ or lower, the step line is mostly straight with small fluctuations due to atomic vibrations, while as we raise the temperature closer to the melting point ($T_{\text{m}}=1327\,\text{K}$) configurational disorder due to capillary fluctuations and the formation of vacancies and adatoms becomes pronounced. The presence of appreciable configurational disorder limits the application of the FL method to temperatures below $700\,\text{K}$. The formation of defects above this temperature causes a sharp increase in the dissipation during the switching to the Einstein crystal. This excessive dissipation is characteristic of irreversible processes and violates the assumptions necessary for the derivation of the FL method equations, namely, the reversibility of the integration path and that the atomic motion occurs around average positions given by the equilibrium lattice positions.

The FL method was applied to both systems in Fig. 2 according to the nonequilibrium techniques presented in Ref. Freitas et al., 2016. We have employed a switching time of $t_{\text{s}}=4\,\text{ns}$ and the S-shaped Freitas et al. (2016); de Koning and Antonelli (1996) functional form for the $\lambda(t)$ parameter. The simulations were carried out at temperatures ranging from $100\,\text{K}$ to $700\,\text{K}$ in intervals of $100\,\text{K}$. Estimates for the statistical errors were obtained by performing three independent switching simulations (forward and backward) for each temperature.

Based on the discussion in Sec. IV.2, the step excess free energy was calculated from the difference of the free energy of the two systems in Fig. 2: $\gamma^{\text{st}}=(F^{\text{st}}-F^{\text{t}})/L$, where $L$ is the total length of the four steps in Fig. 2b. Since the surfaces in the system illustrated in Fig. 2b contain both $\langle 110\rangle$A and $\langle 110\rangle$B types of steps, the FL method provides the average of the free energy of both of these step types. The results of the FL simulations are shown as the red and green dots in Fig. 6. Note that the error bars are smaller than the points on the plot. The standard error of the mean of the points in Fig. 6 is $\approx 1\,\text{meV}/\,\mathring{\mathrm{A}}$, which requires the calculation of the free energy per atom for each system, which is achieved with an accuracy of $\approx 3\,\mu\text{eV/atom}$. Such high statistical accuracy is achievable due to the high efficiency and accuracy of the nonequilibrium Frenkel-Ladd method used in this work. Further details about the technique as well as an in-depth analysis of error control and estimation are provided in Ref. Freitas et al., 2016.

V.2 Step excess quantities

The step excess quantities were calculated for systems with the same size and number of atoms as the systems used for the FL calculations. From the MD simulations we obtained the average energy of the systems with the step $E^{\text{st}}$ and the flat terrace $E^{\text{t}}$, and also the components of stress tensor \mbox{\boldmath\sigma}^{\text{st}} and \mbox{\boldmath\sigma}^{\text{t}}. The step excess properties $[E]_{AN}$ and $[\tau_{\text{avg}}]_{AN}$ were then computed by taking the difference between the quantities of the stepped system and the flat-terrace system, as described in Sec. IV.2. The MD simulations were performed for temperatures ranging from $100$ to $1300\,\text{K}$, in intervals of $100\,\text{K}$. Additionally, we also performed one simulation at the melting temperature for the potential $T_{\text{m}}=1327\,\text{K}$. The systems were equilibrated for $6\,\text{ns}$ before calculating the values of $E^{\text{st}}$, $E^{\text{t}}$, \mbox{\boldmath\sigma}^{\text{st}}, and \mbox{\boldmath\sigma}^{\text{t}}. After equilibration, these values were sampled at intervals of $2\,\text{ps}$ for $400\,\text{ns}$ at each temperature. Figure 7 shows the temperature dependence of $[E]_{AN}$ and $[\tau_{\text{avg}}]_{AN}$. The error bars correspond to the standard error of the mean for each data point, obtained through a block average analysis of the data collected for each temperature. Note that the error bars of $[E]_{AN}$ are too small to appear on the plot.

The results for $[E]_{AN}$ in Fig. 7 show that the excess energy increases with temperature from $(103.61\,\pm\,0.02)\,\text{meV/}\,\mathring{\mathrm{A}}$ at $T=0\,\text{K}$ to $(302\,\pm\,4)\,\text{meV}/\,\mathring{\mathrm{A}}$ at $T_{\text{m}}$. Within a large temperature interval, from zero to approximately $800\,\text{K}$ (i.e., a homologous temperature of approximately $0.60$), the value of $[E]_{AN}$ remains essentially constant, and then begins to increase much more rapidly as the melting temperature is approached. The simulations show that $[E]_{AN}$ remains finite, and does not diverge as the melting point is approached.

The results for excess stress in Fig. 7 show that [\mbox{\boldmath\tau}]_{AN} is appreciably anisotropic at low $T$: the step stress component perpendicular to the step, $[\tau_{\perp}]_{AN}$, is compressive at low temperatures while $[\tau_{\parallel}]_{AN}$ is tensile. Although they are similar in magnitude we notice that they have measurably different values at $0\,\text{K}$: $[\tau_{\perp}]_{AN}=-38.3\,\text{meV/}\,\mathring{\mathrm{A}}$ and $[\tau_{\parallel}]_{AN}=34.3\,\text{meV/}\,\mathring{\mathrm{A}}$. Both parallel and perpendicular components increase with temperature, but the perpendicular does so faster, the consequence being that the anisotropy becomes reduced at high temperatures, where both components become compressive. Notice also that the excess average stress, $[\tau_{\text{avg}}]_{AN}$, remains almost constant for low temperatures, before the onset of large capillary fluctuations of the step. As for $[E]_{AN}$, only for temperatures above approximately $800\,\text{K}$ is a significant temperature dependence of the step excess stress observed.

V.3 Step free energies from thermodynamic integration calculations

In this subsection we focus on the temperature-dependent step free energies, obtained by the thermodynamic-integration (TI) approach described in Sec. III.1. The results obtained from this approach are shown as the solid line in Fig. 6, which plots the value of $\gamma^{\text{st}}$ over the entire temperature range from $T=0\,\text{K}$ up to $T_{\text{m}}$. In performing the TI calculations, we have chosen $T_{0}=400\,\text{K}$ as the reference point for the thermodynamic integration, and the integration was performed in both directions, from $T_{0}$ to $T_{\text{m}}$ and from $T_{0}$ to $\approx 0\,\text{K}$. As noted above, the TI values for $\gamma^{\text{st}}$ agree well with those from the FL method that were not used in the integration (red points in Fig. 6), demonstrating the consistency of the predictions for the temperature dependence of $\gamma^{\text{st}}$ at low homologous temperatures obtained from these two independent methods. We present in the Appendix a discussion of the numerical convergence of the TI results, including error calculations and the independence of the final results on the choice of the reference temperature $T_{0}$.

Overall, the TI results in Fig. 6 show that the temperature dependence of $\gamma^{\text{st}}$ is large and highly nonlinear over the full temperature range. Although the magnitude of $\gamma^{\text{st}}$ remains finite at the melting point, indicating that the surface remains faceted up to $T_{\text{m}}$, the net effect of increasing temperature is a sizable decrease of $\gamma^{\text{st}}$. Specifically, increasing the temperature up to melting leads to a decrease in magnitude of $\gamma^{\text{st}}$ by more than half, from a value of $(103.61\pm 0.02)\,\text{meV/}\,\mathring{\mathrm{A}}$ at $T=0\,\text{K}$ to $(45.8\pm 0.4)\,\text{meV/}\,\mathring{\mathrm{A}}$ at $T=T_{\text{m}}$.

Considering the temperature dependence of $\gamma^{\text{st}}$ in further detail, we divide the results into two temperature ranges: low homologous temperature up to $800\,\text{K}$ (i.e., from homologous temperatures of zero to approximately 0.60), and high homologous temperatures from $800\,\text{K}$ up to the melting point. In the first temperature range, the excess quantities presented in the previous section are approximately constant in value, and $\gamma^{\text{st}}$ displays a relatively weak rate of decrease with temperature. Over this temperature range the value of $\gamma^{\text{st}}$ decreases approximately linearly, by roughly $13\%$ percent, from a value of $(103.61\pm 0.02)\,\text{meV/}\,\mathring{\mathrm{A}}$ to $(90.3\pm 0.2)\,\text{meV/}\,\mathring{\mathrm{A}}$. Since the steps are observed to remain straight on the simulation length and time scales (i.e.., no evidence of appreciable kinks, adatoms, or surface vacancies is observed) for temperatures up to $800\,\text{K}$, we interpret the temperature dependence of $\gamma^{\text{st}}$ over this temperature range to arise primarily from atomic vibrational contributions to the step excess thermodynamic quantities.

Above $T=800\,\text{K}$, $\gamma^{\text{st}}$ displays a much more pronounced temperature dependence. From $800\,\text{K}$ up to $T_{\text{m}}$ the value of $\gamma^{\text{st}}$ decreases by roughly $51\%$ percent, from a value of $(90.3\pm 0.2)\,\text{meV/}\,\mathring{\mathrm{A}}$ to $(45.8\pm 0.4)\,\text{meV/}\,\mathring{\mathrm{A}}$. In this temperature range, the concentration of surface adatoms and vacancies increases significantly, and the magnitudes of the step capillary fluctuations become more pronounced. The larger temperature dependence of $\gamma^{\text{st}}$ over this temperature range is thus interpreted to be a manifestation of the effect of such configurational disorder on the step excess thermodynamic quantities.

V.4 Comparison with previous measured and calculated results

Although we are not aware of previous results presenting the temperature dependence of step free energies in Cu all the way up to the melting point, there have been measurements and previously published calculations at low temperatures for this system, to which the present simulation results can be compared.

The step free energy of Cu$(111)$ $[110]$A and B steps at selected temperatures has been obtained experimentally from the analysis of adatom and vacancy islands observed using scanning tunneling microscopy Schlößer et al. (1999); Schulze Icking-Konert et al. (1999); Giesen et al. (2001); Steimer et al. (2001). For a comprehensive comparison of two available methods for computing $\gamma^{\text{st}}$ experimentally, we refer the reader to Ref. Steimer et al., 2001, where Steimer et al. report $\gamma^{\text{st}}=256\,\pm\,22\,\text{meV/a}$ for an average of $A$ and $B$ steps, where $a$ is the atomic distance along the $[110]$ direction and the measurement is for an average temperature of $\overline{T}=360\,\text{K}$ ($T\in[280,440]\,\text{K}$). The present results are remarkably close to this value, as indicated in Fig. 6: we obtain values of $\gamma^{\text{st}}=254.8\,\pm\,0.2\,\text{meV/a}$ for the same temperature. Moreover, the temperature dependence of $\gamma^{\text{st}}$ shown in Fig. 6 is consistent with the analysis in Ref. Steimer et al., 2001, suggesting that the step free energy has a weak temperature dependence for the temperature range at which the experiments were conducted.

This good agreement between the present simulation results and experimental measurements is achieved despite the approximations inherent in the classical description of the interatomic interactions by an EAM potential model. Importantly, a similar level of agreement is also obtained by the EAM model and available ab initio values at zero temperature obtained by density functional theory (DFT). Specifically, the value of the step energy given by the EAM potential considered in this work is $\gamma^{\text{st}}=264.8\,\text{meV/a}$ at $T=0\,\text{K}$, which agrees well with the DFT result of $\gamma^{\text{st}}=270\,\text{meV/a}$ reported in Refs. Feibelman, 1999; Da Silva et al., 2006; Stasevich et al., 2006. The fact that the current results agree well with DFT at zero temperature, and with experiment at finite temperatures, suggests that the latter agreement is not a result of cancellation of errors resulting from inaccurate energetics and temperature dependencies. Rather the EAM model for Cu of Mishin et al. Mishin et al. (2001) employed in this work appears to yield accurate values for $\gamma^{\text{st}}$ and its temperature dependence, at least at low homologous temperatures.

VI Summary

We present a thermodynamic formalism for steps on faceted surfaces of single-component crystalline solids, resulting in the derivation of a general adsorption equation, Eq. (14), relating changes in step free energy ($\gamma^{\text{st}}$) to variations in chemical potential, surface free energy, temperature, and strain. The rate of change of $\gamma^{\text{st}}$ with respect to variations in these variables is related to surface excess quantities of particle number, surface area, entropy, and stress, respectively. Due to the existence of Gibbs-Duhem relations for the bulk and surface, which give rise to constraints on the variations of the intensive variables, Cramer’s rule can be used to express the adsorption equation in terms of a particular choice for the set of independent variables. The approach results in the definition of step excess quantities formulated in terms of determinants, following the formalism first introduced in the context of interfacial thermodynamics by Cahn Cahn (1979). A direct result of the formulation developed in the present work is the definition of a step excess stress, Eq. (17), which is the step analog of the familiar surface stress quantity, and which represents the excess force on the perimeter of a stepped surface due to the presence of a step. Although the formalism presented in this work is developed only for the special case of single-component crystalline surfaces, the underlying approach is more general, and can be extended to multicomponent/multiphase situations, as demonstrated recently by Frolov and Mishin Frolov and Mishin (2015).

The thermodynamic formalism presented in this work is demonstrated to provide a convenient framework for thermodynamic-integration calculations of the temperature dependence of $\gamma^{\text{st}}$ by atomistic simulations. For this purpose, it is natural to employ a particular choice for the set of independent intensive variables that leads to the definition of step excess quantities, in a manner that is similar to choosing a Gibbs Gibbs (1928) dividing surface leading to zero excess volume and particle number. By combining the resulting expression for the adsorption equation with the Gibbs-Helmholtz relation, we derive an expression for the temperature dependence of the step free energy, Eq. (22), in terms of step excess energy and excess stress quantities that can be readily calculated in atomistic simulations. It is straightforward to extend the proposed TI approach to steps at faceted solid-liquid interfaces, grain boundaries and phase boundaries in multicomponent systems Frolov and Mishin (2015). This approach can provide full temperature and composition dependence of step free energy from atomistic simulations, provided that a reference free energy value is known at some temperature and composition. In the present work we have demonstrated how the Frenkel-Ladd method can be employed for this purpose, when the interfaces of interest involve only solid phases. For solid-liquid or solid-vapor interfaces alternative approaches would be needed such as those based on nucleation simulations (e.g., Ref. Frolov and Asta, 2012) or analyses of capillary fluctuations (e.g., Refs. Hoyt et al., 2001; Mishin, 2014).

We demonstrate the application of the thermodynamic integration formalism for the case of $\langle 110\rangle$ steps on faceted $\{111\}$ surfaces of element Cu, employing MD simulations based on a classical EAM potential due to Mishin et al. Mishin et al. (2001). By combining the thermodynamic-integration formalism with the the Frenkel-Ladd method for computing a reference value of $\gamma^{\text{st}}$ at low temperatures, where the step structure remains highly ordered, we present a calculation of the step free energy over the entire temperature range from zero up to the melting point.

In the process of performing the thermodynamic-integration calculations, we compute temperature-dependent values for the step excess energies and stresses, as shown in Fig. 7. The excess energy is found to display a weak temperature dependence up to a homologous temperature of approximately 0.60; beyond this temperature the excess energy increases strongly as the step displays growing configurational disorder due to the formation of surface adatoms and vacancies and appreciable capillary fluctuations. For the step excess stress, we have obtained negative $[\tau_{\perp}]_{AN}$ and positive $[\tau_{\parallel}]_{AN}$ at low homologous temperatures, with both terms having similar magnitudes. With increasing temperature, the low-temperature anisotropy of the step stress is greatly reduced, and at high temperatures $[\tau_{\perp}]_{AN}$ becomes positive. Therefore, thermal effects such as thermal expansion, vibrational fluctuations, and configurational disordering affect each step stress component differently. It is worth noting that this behavior is not unique to step stresses; it has been observed before in atomistic simulations Becker et al. (2007); Frolov and Mishin (2010c) that the surface stress for solid-liquid interfaces also presents positive and negative values, depending on the system properties and thermodynamic conditions.

For the temperature dependence of the calculated step free energy, our findings are shown in Fig. 6 and can be summarized as follows. At low homologous temperatures (i.e., less than approximately 0.6), where the thermal effects are interpreted to be associated primarily with atomic vibrations, $\gamma_{\text{st}}$ is calculated to display a relatively weak temperature dependence. At these low temperatures, the calculated magnitudes of $\gamma^{\text{st}}$ show good agreement with previously reported experimental measurements and DFT calculations, indicating the accuracy of the employed EAM potential for the present application. The calculated temperature dependence of $\gamma^{\text{st}}$ increases strongly at higher homologous temperatures, as the step becomes increasingly configurationally disordered. The net effect is a reduction in the step free energy by more than half as the temperature is increased from zero up to the melting temperature. Such a strong temperature dependence at high homologous temperatures would be expected to have important consequences for kinetic processes such as surface island nucleation and growth kinetics.

We emphasize that the formalism presented in this work provides a general framework for the calculation of step free energies for elemental systems using atomistic simulation methods, and it is applicable beyond the application demonstrated in this work for elemental Cu modeled by an EAM classical potential. In practical applications to other systems, several considerations should be taken into account. First, the Frenkel-Ladd approach provides a methodology to compute step free energies only at temperatures where contributions of configurational disorder due to kinks, adatoms, and vacancies can be ignored, i.e., where vibrational contributions to the temperature dependence of the excess properties dominate; to ensure that this is the case sufficiently long simulations are required to guarantee structural equilibration, or theoretical analyses based on calculated kink and point-defect formation energies should be performed. For temperatures where the steps remain structurally ordered, the Frenkel-Ladd approach converges sufficiently rapidly that it is expected to be applicable to systems with more complex interatomic potentials, or even within the framework of DFT-based MD simulations, provided large enough systems can be considered to account for the strain fields around the steps and adequate sampling of the phonon spectra. Once reference values have been computed by the Frenkel-Ladd approach, the thermodynamic-integration formalism developed in this work can be used to compute step free energies incorporating configurational and vibrational contributions on an equal footing. In general, such calculations require combinations of efficient interatomic potential models and/or advanced sampling methods to enable equilibration of kink and point-defect densities. Additionally, the contributions due to capillary fluctuations can give rise to large size effects (due to the long-wavelength modes) particularly near the roughening temperature (e.g., Refs. Gelfand and Fisher, 1990 and Block et al., 2014), and to account for these effects calculations with different system sizes and/or analysis of the capillary wave spectra may be necessary. Nevertheless, provided these various considerations are taken into account, the formalism presented in this work provides a framework for computing benchmark results against which theories for vibrational (e.g., Refs. Rahman et al., 2003 and Durukanoğlu et al., 2003) and configurational (e.g., Refs. Nelson et al., 1993 and Gelfand and Fisher, 1990) contributions to the step free energies can be compared. We thus anticipate the approach to be useful for furthering understanding of the thermodynamic properties of steps on crystalline surfaces well beyond the application demonstrated in this paper.

Acknowledgements.

The authors would like to thank Professor J.J. Hoyt for valuable discussions. The research of R.F. and M.A. at UC Berkeley was supported by the US National Science Foundation (Grants No. DMR-1105409 and No. DMR-1507033). R.F. acknowledges additional support from the Livermore Graduate Scholar Program. T.F. acknowledges partial support through a postdoctoral fellowship from the Miller Institute for Basic Research in Science at the University of California, Berkeley. Additional support for T.F. was provided under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344.

Appendix A Error calculation and numerical convergence analysis for the thermodynamic integration results

In this appendix we present analyses of the statistical sampling errors and numerical convergence for the thermodynamic integration results presented in Sec. V.3.

The numerical integration in Eq. (24) was performed considering a linear interpolation of the excess quantities shown in Fig. 7. We have tested interpolation schemes using polynomials of different orders and the difference compared to the linear interpolation was negligible. The reason is that the integrand of Eq. (24) is already smooth for the linear interpolation due to the renormalization of the excess quantities by $T^{2}$ or $T$, as shown in Fig. 8. The second term in the integrand of Eq. (24) (involving the excess stress) was found to be at least 50 times smaller than the first term (involving the excess energy) and therefore it is numerically negligible for the result of the integral.

We have chosen $T_{0}=400\,\text{K}$ as the initial point to perform the thermodynamic integrations to compute the temperature dependence of $\gamma^{\text{st}}$. The integration was performed in both directions, from $T_{0}$ to $T_{\text{m}}$ and from $T_{0}$ to $\approx 0\,\text{K}$. The choice of $T_{0}=400\,\text{K}$ as the initial integration point in Fig. 7 was arbitrary and, within the statistical accuracy of the calculations, it should not influence the final results for $\gamma^{\text{st}}$. The free energy calculated with the FL method at any of the other temperatures (red points in Fig. 6) should all be equally valid as an initial integration point. Thus, to verify the accuracy of the calculations, we have performed the integration in Eq. (24) starting from all the different $T_{0}$ values for which we have available FL simulations. The result is shown in Fig. 9 where we plot the value of $\gamma^{\text{st}}$ at $T_{\text{m}}=1327\,\text{K}$ obtained from the integration of Eq. (24) using different initial points. The error bar of each point corresponds to the error of the mean for the particular value of $T_{0}$. The error in the mean was obtained by a resampling process of the excess quantities and the initial value of $\gamma^{\text{st}}$ used in the integration: each of the data points involved in the integration was picked randomly from a normal distribution with a mean value corresponding to the calculated average value of that quantity, and the standard deviation corresponding the calculated standard error of the mean value. The linear interpolation and numerical integration of the excess quantities for the given choice of $T_{0}$ was performed and the resulting step free energy at $T_{\text{m}}$ was averaged over $2000$ of these resampled data sets. For completeness we also show in Fig. 10 the $\gamma^{\text{st}}(T)$ curve obtained from the integration starting from the different $T_{0}$ values. From Figs. 9 and 10 it is clear that the choice of the initial integration point $T_{0}$ does not influence the final result of the thermodynamic integration significantly.