Capillary fluctuations of surface steps: An atomistic simulation study for the model Cu(111) system

TL;DR

This study uses molecular dynamics simulations to analyze capillary fluctuations of surface steps on copper, revealing isotropic step stiffness and diffusion-governed fluctuation lifetimes, with results compared to thermodynamic-integration data.

Contribution

First atomistic simulation study of capillary fluctuations on Cu(111) steps, providing detailed fluctuation spectra and stiffness measurements near melting temperature.

Findings

Step stiffness is isotropic within statistical error.

Fluctuation lifetimes vary over four orders of magnitude with wave number.

Simulation results differ significantly from thermodynamic-integration calculations.

Abstract

Molecular dynamics (MD) simulations are employed to investigate the capillary fluctuations of steps on the surface of a model metal system. The fluctuation spectrum, characterized by the wave number ($k$) dependence of the mean squared capillary-wave amplitudes and associated relaxation times, is calculated for $\left<110\right>$ and $\left<112\right>$ steps on the $\{111\}$ surface of elemental copper near the melting temperature of the classical potential model considered. Step stiffnesses are derived from the MD results, yielding values from the largest system sizes of $(37\pm1) \, \mathrm{meV}/\mathring{\mathrm{A}}$ for the different line orientations, implying that the stiffness is isotropic within the statistical precision of the calculations. The fluctuation lifetimes are found to vary by approximately four orders of magnitude over the range of wave numbers investigated, displaying a $k$ dependence consistent with kinetics governed by step-edge mediated diffusion. The values for step stiffness derived from these simulations are compared to step free energies for the same system and temperature obtained in a recent MD-based thermodynamic-integration (TI) study [Freitas, Frolov, and Asta, Phys. Rev. B 95, 155444 (2017)]. Results from the capillary-fluctuation analysis and TI calculations yield statistically significant differences that are discussed within the framework of statistical-mechanical theories for configurational contributions to step free energies.

Tables

Table 1. Table 1: Simulation box size used for each step orientation. d 𝑑 d is the step-step distance, h ℎ h the bulk depth, L 𝐿 L the step length, and U 0 subscript 𝑈 0 U_{0} is the step self-energy. The box geometry is illustrated in Fig.  1 .

	$d\left(\mathring{A}\right)$	$h\left(\mathring{A}\right)$	$L\left(\mathring{A}\right)$	$U_0\left(\text{meV}/\mathring{A}\right)$
$\left[110\right]$A	$60.5$	$42.8$	$104.8$	$103.1$
$\left[110\right]$B	$66.6$	$42.8$	$104.8$	$104.1$
$\left[211\right]$	$62.9$	$42.8$	$108.9$	$120.5$

Table 2. Table 2: Step stiffness for different step orientations and convergence with box size. d 𝑑 d is the step-step separation distance, h ℎ h is the bulk depth, L 𝐿 L is the step length, and U 0 subscript 𝑈 0 U_{0} is the step self energy. The dimensions used are with respect to the box sizes presented in Table 1 . The box geometry is illustrated in Fig.  1 .

Step stiffness $\tilde{\sigma }$111All columns have error bar $\pm 1\text{meV/}\mathring{A}$, except for the $U_0$ column where the error bar is $\pm 0.02\text{meV/}\mathring{A}$. $\left(\text{meV/}\mathring{A}\right)$
	$(d,h,L)$222The dimensions are given relative to the values shown in Table 1	$(d,2h,L)$	$(2d,1.7h,L)$	$U_0$
$\left[110\right]$A	$37$	$36$	$37$	$103.13$
$\left[110\right]$B	$38$	$39$	$37$	$104.08$
$\left[211\right]$	$36$	$36$	$37$	$120.53$

Equations

$$[F{]}_{AN}\equiv F^{\text{st}}-F^{\text{t}}={\gamma }^{\text{st}}L,$$

$$[F{]}_{AN}=-k_{\text{B}}T\ln Q,$$

$$Q^{\text{st}}=\int\operatorname{d}\!{\mathbf{R}}\,\exp\bigg{[}-\beta U_{\text{cg}}(\mathbf{R})\bigg{]},$$

$$U_{\text{cg}}(\mathbf{R})=-k_{\text{B}}T\ln\left\{\int\operatorname{d}\!{\mathbf{r}}\,\exp\bigg{[}-\beta U^{\text{st}}(\mathbf{R},\mathbf{r})\bigg{]}\right\},$$

$$Q=\frac{Q^{\text{st}}}{Q^{\text{t}}}=\int\operatorname{d}\!{\mathbf{R}}\,\exp\bigg{[}-\beta\mathcal{H}(\mathbf{R})\bigg{]},$$

$$\gamma^{\text{st}}L=-k_{\text{B}}T\ln\bigg{\{}\int\operatorname{d}\!{\mathbf{R}}\,\exp\big{[}-\beta\mathcal{H}(\mathbf{R})\big{]}\bigg{\}}.$$

$$\mathcal{H}\left[y(x)\right]={\int }_y\sigma (\theta )d\ell ={\int }_0^L\sigma (\theta )\sqrt{1+y^{\prime }(x{)}^2}dx,$$

$$\mathcal{H}[y(x)]=\sigma L+\frac{1}{2}\tilde{\sigma }{\int }_0^L{y}^{\prime }(x{)}^{2}dx,$$

$$y(x)=\sum _{n=-(M-1)/2}^{(M-1)/2}{A}_n\exp \left(i{k}_{n}x\right),$$

$$\mathcal{H}(\{A_n\})=\sigma L+\tilde{\sigma }L\sum _{n=1}^{(M-1)/2}{k}_n^2{\left|A_n\right|}^2,$$

$$\tilde{\sigma }L{k}_n^2⟨{\left|A_n\right|}^2⟩=k_{\text{B}}T,$$

$$⟨{\left|A_n\right|}^2⟩=\left(\frac{k_{\text{B}}T}{\tilde{\sigma }L}\right)\frac{1}{k_n^2}.$$

$$C_n(t)=\frac{⟨f_n(t)f_n(0)⟩}{⟨f_n(0{)}^2⟩}=\exp (-t/{\tau }_n),$$

Full text

I Introduction

Capillary fluctuations are a ubiquitous phenomenon at fluid interfaces, line defects, and crystalline interfaces that are atomically rough Aarts et al. (2004); Hoyt et al. (2001); Mishin (2014); Bartelt et al. (1994); Braslau et al. (1988); Ocko et al. (1994); MacDowell et al. (2014). These equilibrium fluctuations, which lead to variations in the line length of a linear defect, or area of a rough interface at finite temperature, have been widely studied by advanced experimental characterization techniques and computer simulations, as they provide insights into the thermodynamic and kinetic properties of the interfaces on which they form. While a detailed overview of such studies is beyond the scope of the present manuscript, we refer the reader to comprehensive reviews and representative experimental and computational studies Jeong and Williams (1999); Einstein (2014); Müller and Saúl (2004); Bartelt et al. (1990); Constantin et al. (2007); Williams (1994) in the context of steps at faceted crystalline interfaces, which provide the focus of the present work. The properties of such steps play a critical role in governing the kinetics of crystal growth from melt, solution, or vapor phases, due to their influence on the thermodynamics of island nucleation and the kinetics of interface migration (e.g., Ref. Libbrecht (2017)).

Over the last decade analyses of capillary fluctuations in molecular-scale computer simulations, based on molecular dynamics (MD) or Monte Carlo (MC) methods, have been employed extensively within the so-called capillary-fluctuation method (CFM) approach to computing interfacial free energies and their associated crystalline anisotropies for crystal-melt interfaces, grain boundaries, and solid-solid heterophase interfaces (e.g., Refs. Hoyt et al. (2001); Sun et al. (2004, 2006); Karma (1993); Buff et al. (1965); Mishin (2014); Foiles and Hoyt (2006); Hoyt et al. (2010)). Recently, the CFM approach has been employed also for steps at faceted crystal-melt interfaces Saidi et al. (2017) to derive temperature-dependent step stiffnesses, which are relevant in the context of modeling solidification rates and associated crystal growth morphologies. For crystal-melt interfaces, liquid surfaces, and fluid-fluid interfaces, detailed comparisons of CFM results with those obtained using alternative thermodynamic-integration and nucleation based MD methods have been undertaken to understand the range of applicability and associated accuracies of these alternative approaches (e.g., Refs. Sides et al. (1999); Ismail et al. (2006); Lacasse et al. (1998); Davidchack et al. (2006); Davidchack and Laird (2000, 2003); Morris and Song (2003)). At the present time we are unaware of such comparisons for the applications of the CFM for step properties.

In the present work we consider the application of the CFM approach for studying thermodynamic and kinetic properties of steps on crystalline surfaces, focusing on Cu$(111)$ as a representative model metal system. The results of equilibrium MD simulations near the melting temperature of the potential model considered are analyzed to compute step fluctuation spectra, characterized by the wave number ($k$) dependence of the mean-square amplitudes $\langle\left|A(k)\right|^{2}\rangle$, as well as the fluctuation relaxation times $\tau(k)$. From the dependence of $\langle\left|A(k)\right|^{2}\rangle$ on $k$ we derive step stiffnesses for $\langle 110\rangle$ and $\langle 112\rangle$ step orientations, obtaining values that are isotropic (independent of orientation) within the statistical precision of the simulations. Further, we obtain values of fluctuation lifetimes that are consistent with the $k^{-4}$ scaling associated with dynamics that are governed by step-edge diffusion.

The focus on the Cu$(111)$ system in the present study enables a comparison of CFM results with step free energies obtained in a recent study published by the authors Freitas et al. (2017) using an alternative thermodynamic-integration (TI) approach. The values of the step stiffnesses derived by the CFM are lower by approximately $25$% compared with the step free energies calculated by the TI approach for $\langle 110\rangle$ oriented steps at the same temperature. The discrepancy is discussed within the framework of statistical-mechanical theories of the configurational contributions to step free energies (e.g., Refs. Gelfand and Fisher (1990); Kayser (1986)) associated with capillary fluctuations.

The remainder of this paper is organized as follows. In Sec. II we present a brief derivation of the main results from capillary-wave theory that are used in the remainder of the paper; although similar derivations appear already in many places in the literature, the overview is included to emphasize key concepts and equations required for the analysis and interpretation of the present MD results. In Sec. III we describe the details of the MD simulations of step capillary fluctuations and in Sec. IV we present the simulation results. In Sec. V a discussion is presented focusing on the comparison of the present CFM results to step free energies obtained previously by thermodynamic integration Freitas et al. (2017). Finally, the results and conclusions are summarized in Sec. VI.

II Capillary-wave model

In this section we summarize the main equations required for CFM analysis of surface step fluctuations. We describe how the capillary-wave Hamiltonian results from the coarse-graining of the atomic partition function, and also highlight several nuances of the CFM that will be discussed in the context of the analysis of the MD simulation results in Sec. III.

II.1 Step effective Hamiltonian

Following the notation from Ref. Freitas et al. (2017), the excess free energy of a step can be defined thermodynamically through the relation:

[TABLE]

where $\gamma^{\text{st}}$ is the step free energy per unit length and $L$ is the system dimension along the average step direction, as illustrated in Fig. 1. $F^{\text{st}}$ and $F^{\text{t}}$ are the absolute free energies of systems with the same surface area ($A$), number of atoms ($N$), and temperature ($T$). These systems can be considered to be identical except that the system corresponding to $F^{\text{t}}$ has a flat surface, while the one corresponding to $F^{\text{st}}$ contains a surface step of length $L$. The free energy of the system with a flat surface can also be written as $F^{\text{t}}=-k_{\text{B}}T\ln Q^{\text{t}}+3Nk_{\text{B}}T\ln\Lambda$, where $k_{\text{B}}$ is the Boltzmann constant, $Q^{\text{t}}$ is the configurational part of the system’s partition function, and $\Lambda=(h^{2}/2\pi mk_{\text{B}}T)^{1/2}$ is the thermal de Broglie wavelength. Similarly, the free energy of the system with a step is $F^{\text{st}}=-k_{\text{B}}T\ln Q^{\text{st}}+3Nk_{\text{B}}T\ln\Lambda$. Hence, we can rewrite the step free energy in Eq. (1) as follows:

[TABLE]

where $Q\equiv Q^{\text{st}}/Q^{\text{t}}$ is the ratio of the configurational partition functions.

In order to clarify the physical meaning of $Q$ consider the potential energy of the system with a step: $U^{\text{st}}(\mathbf{x})$, where $\mathbf{x}$ is the $3N$-dimensional vector with the atomic coordinates. We can perform a canonical transformation on the atomic coordinates and separate the variables describing the step configuration from all other variables. With this transformation the potential energy can be written as $U^{\text{st}}=U^{\text{st}}(\mathbf{R},\mathbf{r})$, where $\mathbf{R}$ are the step degrees of freedom and $\mathbf{r}$ represents all other degrees of freedom (i.e., bulk and surface degrees of freedom). With this set of generalized coordinates the configurational partition function of the system with the step can be written as

[TABLE]

where

[TABLE]

is a coarse-grained potential energy which involves only the step degrees of freedom. Notice that, for convenience, we have performed the canonical transformation in such a way as to render $\mathbf{R}$ and $\mathbf{r}$ dimensionless quantities. Equation (4) implies that $U_{\text{cg}}(\mathbf{R})$ is the portion of the free energy associated with the bulk and surface configurational degrees of freedom $\mathbf{r}$. Alternatively, Eq. (3) suggests that $U_{\text{cg}}(\mathbf{R})$ can also be seen as the potential that generates the step dynamics on that system. Because of this last interpretation $U_{\text{cg}}(\mathbf{R})$ is also known as the potential of mean force Tuckerman (2008), i.e., it is the potential acting on the step that arises from the mean contribution of the bulk and surface degrees of freedom. In the limit of adiabatic decoupling between the step and the rest of the system $U_{\text{cg}}(\mathbf{R})$ becomes an effective potential on which the step degrees of freedom ($\mathbf{R}$) can be assumed to evolve in time independently from the other degrees of freedom ($\mathbf{r}$).

The step free energy can be written as a function of the coarse-grained potential energy $U_{\text{cg}}(\mathbf{R})$. First, we substitute Eq. (3) in the equation for $Q$ and use $Q^{\text{t}}=\exp\left[-\beta(F^{\text{t}}-3Nk_{\text{B}}T\ln\Lambda)\right]$:

[TABLE]

where we have defined the step effective Hamiltonian $\mathcal{H}(\mathbf{R})\equiv U_{\text{cg}}(\mathbf{R})-(F^{\text{t}}-3Nk_{\text{B}}T\ln\Lambda)$. Now the step free energy can be obtained from Eqs. (1) and (2):

[TABLE]

Notice that Eq. (5) does not involve any approximation, we have only separated and interpreted specific parts of the partition function $Q$. Hence, the calculation of the step free energy $\gamma^{\text{st}}$ using Eq. (5) still involves an integral over the phase space of all particles.

II.2 Capillary-wave model for steps

It is now possible to introduce a model for the step effective Hamiltonian, $\mathcal{H}(\mathbf{R})$, that simplifies the calculation of Eq. (5) but still includes all relevant physical properties that govern the step dynamics. A reasonable model that forms the basis for capillary-wave theory (e.g., Ref. Gelfand and Fisher (1990)), is to assume that a fluctuation of the step line that causes a change $\delta\ell$ in step length has an energetic cost of $\sigma\delta\ell$, where $\sigma$ is the step energy per unit length. With this physical picture the step effective Hamiltonian takes the form:

[TABLE]

where $\theta(x)=\tan^{-1}(y/x)$ is the step orientation with respect to the average step-line direction and the integral is over the curve $y(x)$ describing the step-line profile, as illustrated in Fig. 1. Thus, $\mathcal{H}$ is a functional of the step configuration $y(x)$ Kardar (2007).

Notice that $\sigma(\theta)$ defined in Eq. (6) is different from the step stress tensor \mbox{\boldmath\tau}^{\text{st}} as defined in, for example, Refs. Freitas et al. (2017); Li et al. (2011). The step stress tensor couples mechanically to the system strain and gives origin to a elastic deformation energy which can be directly measured in atomistic simulations Freitas et al. (2017). The physical interpretation of $\sigma(\theta)$ is more complicated, as discussed in detail in Ref. Gelfand and Fisher (1990). For example, $\sigma(\theta)$ reflects the energy per unit physical length of the step, while \mbox{\boldmath\tau}^{\text{st}} and $\gamma^{\text{st}}$ are defined per unit length of the average step direction, indicated as $L$ in Fig. 1 and Eq. (1). For our purpose in this paper we will refer to $\sigma(\theta)$ as the step tension in the line-fluctuation model given by Eq. (6).

To compute from Eq. (6) the equilibrium spectrum for the capillary fluctuations, and the resulting free energy, the traditional approach Buff et al. (1965); Gelfand and Fisher (1990); Kardar (2007) is to make use of the small slope approximation where $\theta(x)\approx y^{\prime}(x)$ and the terms inside the integral in Eq. (6) can be expanded in powers of $y^{\prime}(x)$. Collecting the terms with the same power and keeping only terms $O(y^{\prime 2})$ allow us to write the step effective Hamiltonian as

[TABLE]

where $\sigma\equiv\sigma(0)$ is the step tension of the state with a straight step in this model, and $\tilde{\sigma}\equiv\sigma(0)+\sigma^{\prime\prime}(0)$ is the step stiffness, where $\sigma^{\prime\prime}(0)$ denotes the second derivative of the step tension with respect to the orientation of the step normal evaluated in the state where the step is straight. The term $\sigma L$ in Eq. (7) is the energy of a straight step, while the second term is the energy penalty in having any curvature along the step line, i.e., the energy cost of step fluctuations.

The next step is to discretize the integral in Eq. (7) into a Riemann sum, resulting in a Hamiltonian that is quadratic in $y(x_{n})$, with $x_{n}=n\Delta x$ where $n=0,1,\ldots,M-1$ and $\Delta x=L/M$. In what follows we adopt a similar approach based on a Fourier representation of the step profile. This formulation, while equivalent, leads to expressions more aligned with the CFM analysis of computer simulation results.

The step line profile $y(x)$ shown in Fig. 1 can be decomposed in normal modes as:

[TABLE]

where the wavevectors $k_{n}$ are given by $k_{n}=n(2\pi/L)$ with $n=0,\pm 1,\pm 2,...,\pm(M-1)/2$ (assuming $M$ is odd). Using Eq. (8) we can compute the integral in Eq. (7) and obtain:

[TABLE]

where we have made use of the fact that $A^{*}_{n}=A_{-n}$ since $y(x)$ is real. In this system of coordinates the amplitudes $A_{n}$ of the normal modes are the step degrees of freedom since they define the step configuration $y(x)$ through Eq. (8). The step effective Hamiltonian given by Eq. (9) is quadratic in all its degrees of freedom and thus many properties of the system can be obtained exactly.

It is clear from Eq. (9) that the properties of this system depend on how the step is coarse-grained, i.e., how closely spaced ($\Delta x=L/M$) are the $M$ points describing the step line. This is a reflection of the number of degrees of freedom attributed to the step effective Hamiltonian Gelfand and Fisher (1990); Kayser (1986), Eq. (6), as discussed in Sec. II.1; $\Delta x$ determines the largest wavevector considered in the effective Hamiltonian of Eq. (9): $k_{\text{max}}\equiv\pi/\Delta x$. An extensive review of the consequences and interpretations of this dependency on $\Delta x$ and, consequently, on the capillary-wave wavelengths considered, is given in Ref. Gelfand and Fisher (1990). We return to this point in Sec. V when comparing the results of the CFM analysis to the values of the step free energy computed in Ref. Freitas et al. (2017).

Notice that the term $\sigma L$ in Eq. (9) is a simple shift in energy, thus the dynamics of the step is completely parametrized by the step stiffness $\tilde{\sigma}$. In the next section we review how $\tilde{\sigma}$ can be derived from atomistic simulations.

II.3 Step fluctuation spectrum

According to the equipartition theorem each quadratic degree of freedom in the Hamiltonian of a system at constant temperature $T$ contributes $k_{\text{B}}T/2$ to the system’s average energy. We can apply this theorem to Eq. (9) since each mode amplitude $A_{n}$ appears quadratically in the Hamiltonian. Notice that the real and imaginary parts of $A_{n}$ are independent and, thus, each part contributes with $k_{\text{B}}T/2$ to the total energy. Hence

[TABLE]

where $\langle\ldots\rangle$ indicates a canonical ensemble equilibrium average. This equation can be used to compute the step stiffness, $\tilde{\sigma}$, if we rewrite it as

[TABLE]

The normal-mode amplitudes can be obtained from atomistic simulations and used to adjust a curve of $\langle\left|A_{n}\right|^{2}\rangle\;\text{versus}\;k_{n}^{-2}$, from which $\tilde{\sigma}$ can be extracted.

When using Eq. (10) to compute $\tilde{\sigma}$ it is necessary to define the step profile, shown as $y(x)$ in Fig. 1. Hence, the value of $\tilde{\sigma}$ obtained can be sensitive to how $y(x)$ is determined, particularly if one relies on normal modes with wavelengths comparable to the atomic spacing. The physical origin of this arises because the distinction between interface and bulk degrees of freedom is not clear for atomistic systems Gelfand and Fisher (1990), i.e., the definition of the interface position from the atomic configuration is ambiguous. In Sec. III.3 we study the inherent ambiguity in defining the step configuration in atomic-scale simulations and discuss how $y(x)$ can be determined in such a way as to minimally affect the value of $\tilde{\sigma}$ obtained from Eq. (10).

III Methodology of atomistic simulations

III.1 Molecular dynamics simulations and system geometry

Molecular dynamics simulations of surface steps were performed using the LAMMPS Plimpton (1995) (Large-scale Atomic/Molecular Massively Parallel Simulator) software. The interatomic interactions were described by the embedded-atom method Daw and Baskes (1984) for a system of pure copper Mishin et al. (2001) and the Langevin thermostat Schneider and Stoll (1978) was used to sample the particles’ phase space according to the canonical ensemble distribution. The thermostat relaxation time was $\tau_{\text{L}}\equiv m/\gamma=2\,\text{ps}$, where $\gamma$ is the friction parameter and $m$ the atomic mass. The timestep ($\Delta t$) for the integration of the equations of motion was chosen based on the phonon spectrum of the system; we used $\Delta t=2\,\text{fs}$ which is approximately $1/60$th of the oscillation period of the highest-frequency normal mode of this system.

The geometry of the simulation box is illustrated in Fig. 1. Periodic boundary conditions were used for the directions parallel to the surface ($\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$) while free boundaries were used along $\mathbf{\hat{z}}$ to create the system surface. We kept the box length along the $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ directions fixed while the system fluctuates freely along $\mathbf{\hat{z}}$ (normal to the surface) in order to guarantee mechanical equilibrium with the vacuum. A non-orthogonal simulation box Shilkrot and Srolovitz (1996) was employed in such a way as to have only one step on the system surface.

We have chosen the $(111)$ surface of face-centered cubic copper as a representative metal surface to study steps. The surface properties of the interatomic potential employed have been studied extensively previously Frolov and Mishin (2009); Freitas et al. (2017). Of relevance for this study is the fact that the $\{111\}$ surfaces have been found to be faceted at all temperatures up to the melting point of this model ($T_{\text{m}}=1327\,\text{K}$). In order to study effects of anisotropy we will consider steps of two different orientations (Fig. 2), with step line directions along $[110]$ and $[211]$. The $[110]$ direction presents two distinct steps, $[110]$A and $[110]$B, which have different nearest-neighbor configurations on the layer immediately below the step, as illustrated in Fig. 2. The temperature for all simulations was $T=1300\,\text{K}$; it was chosen to be close to the melting point so that the step fluctuation timescales were compatible with the short physical times accessible to the MD simulations. In Sec. III.4 we present an analysis of the step fluctuation relaxation times to assess which modes are adequately sampled.

III.2 System dimensions

Surface steps deform the crystalline lattice around them creating an elastic field that can interact with other elastic fields present in the crystal. The total energy of an isolated step will be referred to as the self-energy ($U_{0}$). Because of the periodic boundary conditions applied in the simulations the step can interact with its periodic images, as well as with the surface at the bottom of the simulation box (Fig. 1). The effects of these interactions can be made negligibly small by choosing system dimensions such that the interaction energy is much smaller than $U_{0}$.

The step-step interaction energy decreases with the distance $d$ between the steps as Müller and Saúl (2004) $E_{\text{int}}\propto d^{-2}$ and can be attractive or repulsive depending on the orientation of the steps. To determine the magnitude of this interaction we followed the approach from Ref. Shilkrot and Srolovitz (1996) and computed the step-step interaction energy for different step separations. We have chosen the step-step separation distance for our simulations as the minimum distance such that $E_{\text{int}}/U_{0}\leq 10^{-4}$. In practice this resulted in distances $d\approx 60\,\mathring{\mathrm{A}}$ as shown in Table 1.

The bulk depth ($h$ in Fig. 1 and Table 1) was determined by considering the step elastic-field decay along $\mathbf{\hat{z}}$, the direction normal to the surface. The step elastic energy ($E^{\text{bulk}}_{\text{step}}$) decays with the bulk depth proportionally to $\exp(-h/\xi)$ where $\xi$ is a characteristic length which depends on the step orientation and length. After we verified this relationship we used it to impose the same energy tolerance used for the step-step interaction, i.e., $E_{\text{step}}^{\text{bulk}}/U_{0}\leq 10^{-4}$. The selected bulk depth values are shown in Table 1. In all simulations presented here the last six layers of $(111)$ planes at the bottom of the simulation box were frozen at their equilibrium position to guarantee that no bending of the structure would occur. The frozen layers were added beyond the values of bulk of depth $h$ shown in Table 1.

In order to determine the simulation system dimension corresponding to the step length $L$ it is necessary to consider the assumptions of the CFM, as presented in Sec. II. This model is not valid for the description of the step at scales smaller than its coarse-graining scale ($\Delta x$), thus we need $L\gg\Delta x$. However, it would be computationally unfeasible to have a step that is excessively large since the relaxation time of normal modes with long wavelength can be very long on MD time scales due the long-range atomic diffusion necessary to change the configuration of these modes. Therefore, it is necessary to study the normal-mode relaxation times before determining what is a satisfactory step length. The details of the analysis of these relaxation times is presented in Sec. III.4, and based on the results we have chosen $L\approx 100\,\mathring{\mathrm{A}}$ as shown in Table 1. This step length is equivalent to the largest relaxation time of the normal mode (with largest wavelength) being $\tau_{\text{max}}\approx 1\,\text{ns}$.

III.3 Step profile determination

Given any interface between two distinct phases there is no unambiguous approach to determine the interface position from the microscopic atomistic structure of the system Gelfand and Fisher (1990); Delgado-Buscalioni et al. (2008). Therefore, there is no algorithm that uniquely defines the step position, i.e., the one-dimensional interface separating two surface terraces. Here we compare two different algorithms Benet et al. (2014); Delgado-Buscalioni et al. (2008); Jorge et al. (2010); Davidchack et al. (2006); Davidchack and Laird (2000) that determine the step profile from the atomic configurations captured in MD simulations.

Based on the dimensions for $d$ and $h$ presented in Table 1 we have constructed a system with a $[110]$A step with step-step distance $d=60.5\,\mathring{\mathrm{A}}$, bulk depth $h=42.8\,\mathring{\mathrm{A}}$, and step length $L=403.6\,\mathring{\mathrm{A}}$. This system was equilibrated for $10\,\text{ns}$ and, afterwards, the step configuration was captured every $0.2\,\text{ps}$ for $100\,\text{ns}$. All results presented in this section were obtained from this simulation. From the MD snapshots the atoms belonging to the top surface layer could be readily identified by counting the number of $(111)$ planes and selecting all atoms with height above some threshold height based on the interplanar separation. One snapshot of the result of this selection is shown in Fig. 3(a); from snapshots like this one we want to define the step line profile, being careful to not select any adatom belonging to the surface and also to not accidentally exclude atoms belonging to the step.

The first algorithm used will be referred to as the “grid algorithm”. In this algorithm the direction along the step length, $\mathbf{\hat{x}}$, is divided in equally spaced bins or delimiting strips, as illustrated in Fig. 3(b). We further divide each of these strips along $\mathbf{\hat{y}}$, creating rectangular cells, and calculate the density of atoms in each cell. The step height is chosen as the average value of the height of the first bin with zero density and the bin immediately before it. The parameters chosen for the dimension of each bin was $4.5\,\mathring{\mathrm{A}}$ parallel to the step line and $6.7\,\mathring{\mathrm{A}}$ perpendicular to the step line.

The second algorithm used is referred to as the “cluster algorithm”, illustrated in Fig. 4. Once again we start with the atomic configuration of the first layer, then determine all atomic clusters on this layer. Two atoms are considered to belong to the same cluster if there is a path between them through a sequence of neighbor atoms, where we consider two atoms to be neighbors if the distance between them is equal or smaller than some maximum radius $r_{\text{max}}$. The step atoms are then defined as the largest cluster of atoms, as shown in Fig. 4(b). From the configuration of the atoms belonging to the step the surface can be readily divided into strips along the $\mathbf{\hat{x}}$ direction and the atom with the highest value of $y(x)$ within that strip is selected to be the step height at that point. The cluster algorithm has the advantage of having only one adjustable parameter, namely, the maximum nearest-neighbor distance $r_{\text{max}}$, which can be easily estimated by considerations of the crystal lattice geometry. We have taken $r_{\text{max}}=1.2r_{\text{n}}$, where $r_{n}$ is the distance between nearest neighbors in the lattice, and the discretization length along the step line was $2.7\,\mathring{\mathrm{A}}$.

We have optimized both algorithms with respect to the parameters involved to obtain step profiles that best adjust to the real atomic configurations. Then we performed the Fourier transform of the height profiles and calculated the power spectrum (i.e., $\left|A_{n}\right|^{2}\,\text{versus}\,k_{n}$). The comparison of the algorithms is shown in Fig. 5 along with a straight line of slope $-2$. The agreement of the power spectrum with the $k_{n}^{-2}$ behavior predicted by the CFM in Eq. (10) is observed, this is an indication that the theory is adequate to describe the step fluctuations at the wavelengths probed in the MD simulations. From Fig. 5 we also see that the long-wavelength modes (small $k_{n}$) are insensitive to the choice of coarse-graining algorithm, as observed before in this type of analysis Jorge et al. (2010); Benet et al. (2014). At intermediate values of $k_{n}$ the cluster algorithm is observed to follow the CFM prediction, Eq. (10), to higher wave numbers than the grid algorithm. Hence, the cluster algorithm is preferred for capturing the step profile details, and this algorithm has been used for analyses of step fluctuations in the remainder of this paper. We also see that for large $k_{n}$ both algorithms deviate from the $k_{n}^{-2}$ behavior. This happens because when the normal-mode wavelength becomes comparable to the interatomic distance the atomic vibrations start to interfere with the step oscillations. The CFM was proposed to describe the long-wavelength capillary waves but it does not account for the discrete nature of the atomic configuration and degrees of freedom, thus it is no surprise that when these effects start to become significant (large $k_{n}$) our results start to deviate from the CFM. We further discuss the validity of the CFM to describe surface steps in Sec. IV.

III.4 Step normal mode relaxation times

Each normal mode in Eq. (10) has a different relaxation time since mass transport is required for changes in the step configuration. Short wavelength modes can change their configuration quickly since they only require short-range diffusion to modify their amplitudes, while modes with small $k_{n}$ have long relaxation times that limit the statistics for their sampling in the MD simulations. Thus, these relaxation times ultimately place a limit on the wavelengths that can be probed in the simulations.

The relaxation time of the normal modes are analyzed from the simulation data employing a time autocorrelation function of the amplitudes, i.e., if $f_{n}(t)=\left|A_{n}(t)\right|^{2}$ then the autocorrelation function $C_{n}(t)$ is

[TABLE]

where $\tau_{n}$ is the relaxation time of the normal mode of wavevector $k_{n}$. Representative plots of $C_{n}(t)$ for selected normal modes are shown in Fig. 6, and from such data we can estimate the relaxation time of each mode, as shown in Fig. 7. The data used to obtain Figs. 6 and 7 was extracted from the simulation performed in Sec. III.3 for a $[110]A$ step.

Notice in Fig. 7 that the MD simulations resulted in $\tau_{n}\propto k_{n}^{-4}$. This result indicates that the step capillary fluctuations are predominantly governed by atomic diffusion Constantin et al. (2007) along the step line, as opposed, for example, to diffusion of adatoms on the terrace, which would lead Karma (1993) to $\tau_{n}\propto k_{n}^{-3}$, or adatom attachment or detachment to or from the step edge that would be consistent with $\tau_{n}\propto k_{n}^{-2}$.

From the relaxation times obtained in Fig. 7 the step length appropriate for the MD simulation cells is determined as follows. We have limited the maximum relaxation time to be $\tau_{\text{max}}=1\,\text{ns}$, with this limitation the shortest wavevector we can sample is $k_{\text{min}}\approx 0.11\,\mathring{\mathrm{A}}^{-1}$ and the shortest step to contain this wavevector has a length of $L\approx 57\,\mathring{\mathrm{A}}$. Because we are making conservative choices for the step length of all other orientations, and to increase the number of points used in the fitting of Eq. (10), we have chosen to use step lengths of approximately $100\,\mathring{\mathrm{A}}$. The dimensions for the simulation cells used to obtain the results presented below are listed in Table 1.

IV Results

Plotted in Fig. 8 are the results of the MD calculated fluctuation amplitudes $\langle\left|A_{n}\right|^{2}\rangle\,\text{versus}\,k_{n}$, obtained as described in Sec. III.3 with the cluster algorithm used to characterize the instantaneous step profile. Using the normal modes with wavevectors $k_{\text{min}}<k_{n}<k_{\text{max}}$, where $k_{\text{min}}=0.110\,\mathring{\mathrm{A}}^{-1}$ and $k_{\text{max}}=0.335\,\mathring{\mathrm{A}}^{-1}$ (shown as dotted lines in Fig. 8), we have verified the $\langle\left|A_{n}\right|^{2}\rangle\propto k_{n}^{-2}$ behavior predicted by Eq. (10) by adjusting a general power law to these points, as shown in the inset of Fig. 8.

As explained in Sec. III.4 the $k_{\text{min}}$ value was determined based on the largest relaxation time that can be adequately sampled in the MD simulations we performed. It is clear in Fig. 8 that as $k_{n}$ decreases below $k_{\text{min}}$ the error bars become larger due to the reduced sampling statistics, associated with the longer relaxation times and the total simulation time of $100\,\text{ns}$.

The value of $k_{\text{max}}$ was chosen by comparing the MD results in Fig. 8 to a curve with slope $k_{n}^{-2}$ and visually deciding at which point the data started to diverge from this behavior, $k_{\text{max}}$ was selected as the average of the $k_{n}$ value of that point and the one immediately before it. Although the choice of $k_{\text{max}}$ is not unique (e.g., it depends on the method used to characterize the step profile), it is not completely arbitrary either and reasonable estimates can be made by inspection of the MD results, as shown in Fig. 8. Specifically, it is clear from Fig. 8 that at large enough $k_{n}$ the MD data deviates from the $k_{n}^{-2}$ behavior predicted by the CFM, and the point where this discrepancy becomes statistically significant provides the basis for a reasonable estimate of a lower bound for $k_{\text{max}}$. Hence, an important criteron is to adjust the $k_{n}^{-2}$ curve in such a way that the data for small $k_{n}$ lies accurately on the $k_{n}^{-2}$ curve. With that reference, the only arbitrariness comes from deciding where the data starts to deviate from the curve location imposed from the small $k_{n}$ points data. Here the cutoff $k_{\text{max}}=0.335\,\mathring{\mathrm{A}}^{-1}$ corresponds to a wavelength of $\lambda\approx 7.2a$, where $a$ is the atomic distance along the step line for steps along $\left<110\right>$ directions. This result implies that it is necessary to average over approximately seven atoms to eliminate noise due to atomic vibrations and correlated atomic displacements, to correctly capture the expected capillary wave behavior, consistent with the coarse-graining over atomic degrees of freedom necessary to define the step effective Hamiltonian $\mathcal{H}$ described in Sec. II.

According to Eq. (9) the step effective Hamiltonian depends on two parameters: the step tension of a straight step ($\sigma$) and the step stiffness ($\tilde{\sigma}$). For each step orientation listed in Table 1 we have run a $100\,\text{ns}$ simulation preceded by a $5\,\text{ns}$ equilibration period and applied Eq. (10) to obtain the step stiffness. The result is shown in the first column of Table 2. We have employed three different box sizes to test for size convergence: one with twice as much bulk depth and another with twice as much step-step distance (increasing the bulk depth to account for the deeper penetration of the step elastic field). Also listed in Table 2 are the step energies $U_{0}$ at $T=0\,\text{K}$.

The zero-temperature values $U_{0}$ are observed to show a significant anisotropy between the $\langle 110\rangle$ and $\langle 211\rangle$ orientations: a difference of $\approx 14\%$ characterizes these values in Table 2. By contrast, at $1300\,\text{K}$ the MD results yield stiffness values that are isotropic within the statistical precision of the MD data. These results imply that $\tilde{\sigma}\approx\sigma$, i.e., that the contribution of $\sigma^{\prime\prime}$ to the stiffness is negligible. The observed decrease in anisotropy of the step tension is interpreted to reflect the fact that when the step fluctuation amplitudes become large on the scale of the atomic dimensions the effects of the lattice are averaged out.

V Discussion

In this section the present CFM results are compared with values of the step free energy for the same system obtained by thermodynamic-integration methods previously Freitas et al. (2017). In Fig. 9 we plot the results from the previous TI calculations, where the solid black line gives the temperature dependence of the step free energy up to the melting point for a $\langle 110\rangle$ step orientation and the red circles are independent results obtained from TI calculations using the Frenkel-Ladd method. Also plotted in Fig. 9 with the diamond symbol is the step tension $\sigma$, which is assumed isotropic based on the MD results presented in the previous section, and thus equal to the step stiffness. It can be seen that the value of $\sigma$ obtained from the CFM analysis of the present MD data is lower than the step free energy obtained in Ref. Freitas et al. (2017) by approximately $25\%$: at $T=1300\,\text{K}$ the TI results yield $\gamma^{\text{st}}=(50.4\pm 0.4)\,\text{meV/}\,\mathring{\mathrm{A}}$, while the CFM yields $\sigma=(37\pm 1)\,\text{meV/}\,\mathring{\mathrm{A}}$. We discuss this discrepancy in what follows in the context of statistical-mechanical theories of capillary fluctuations (e.g., Refs. Gelfand and Fisher (1990); Buff et al. (1965); Kayser (1986); Abraham (1981)).

We begin by noting the differences in the way the two quantities are defined. In the TI (thermodynamic) formalism, which is Gibbsian in spirit, step free energy is defined as an excess free energy per unit length of the system (simulation block in this case). This thermodynamic formalism does not rely on or characterize the physical length of the fluctuating step. The configurational free energy contributions associated with these fluctuations are naturally included in the TI method. By contrast, in the capillary-wave theory the step tension $\sigma$, introduced in Eq. (6), is defined as a free energy per unit physical length of the fluctuating step. These considerations alone suggest that the step free energy obtained by the TI method and the step tension defined in the CFM are inherently different. Moreover, in the TI and CFM simulations the average physical length of the step is larger than that of the system dimension along the step by approximately $30\%$ to $35\%$ for the system sizes considered in Ref. Freitas et al. (2017). In the TI formalism one could in principle introduce the total excess step free energy per unit of average physical length, which in this case would be approximately $23\%$ smaller that the value of $\gamma^{\text{st}}$ given above.

To formalize the difference between the step stiffness and the step free energy, we note that in the literature, starting with the work of Buff et al. (1965), there is a distinction drawn between $\sigma$, often referred to as the “bare” stiffness, and $\gamma^{\text{st}}$, the step free energy. In these theories, the latter differs from the former due to the configurational free energy contributions associated with the step fluctuations. Formally, this difference can be derived by the use of Eq. (9) in Eq. (5), resulting in an expression for the configurational free energy that depends on the number of modes (coarse-graining length), as well as a length scale that is used in Eq. (5) to make the partition function dimensionless. These issues are discussed at length by Kayser (1986), who derives an expression for the configurational contribution to $\gamma^{\text{st}}$ that depends on system size and geometry. Importantly, this contribution can be shown to be always negative, and thus lowers the magnitude of the step free energy ($\gamma^{\text{st}}$) relative to the bare stiffness ($\sigma$). We note that, as shown in Fig. 9, we find the opposite trend, such that this cannot be the explanation for our finding that $\sigma$ derived from CFM analysis of the MD data is lower than $\gamma^{\text{st}}$ derived from thermodynamic integration.

To further consider the origins for this difference, we note that one possibility is that the TI results in Ref. Freitas et al. (2017) could suffer from hysteresis effects associated with sharp changes in the excess quantities as the step configuration evolves from being straight at low temperatures to rough at higher temperatures. We have investigated these issues in detail in our previous work and concluded that for the system sizes and time scales considered in Ref. Freitas et al. (2017) the excess quantities behave smoothly and no evidence of artifacts that would bias the TI integration was found.

We consider then an alternative explanation for the discrepancy between the TI derived value of $\gamma^{\text{st}}$ and CFM derived value of $\sigma$ shown in Fig. 9. Specifically, as discussed by Gelfand and Fisher (1990) and also noted by Kayser (1986), the theoretical analysis of Abraham (1981) for one-dimensional line interfaces in the 2D Ising model shows that the stiffness that governs the growth the mean-square width of the step with system size is exactly equal to the interfacial free energy. This result suggests a limitation to the classical capillary-fluctuation theory where the two quantities $\sigma$ and $\gamma^{\text{st}}$ are distinct. This limitation is discussed by Kayser who argues that the stiffness that governs step fluctuations should depend on the wave number $k$ of the fluctuation, i.e., Kayser argues that large-wavelength fluctuations “see a ‘renormalized’ surface tension” that differs from the bare stiffness due to the configurational degrees of freedom associated with the smaller wavelength fluctuations. In this picture, it could be possible that the stiffness governing step fluctuations for the lowest wave numbers considered in our simulations, which dominate our fitting of the MD data to extract the stiffness values, is lower than the value of $\gamma^{\text{st}}$ obtained from TI in Ref. Freitas et al. (2017) which considered considerably smaller step lengths. These considerations suggest an interesting direction for future work that would involve detailed analysis of the size dependence of both TI and CFM results in the calculation of step free energies.

VI Summary and Conclusions

Molecular dynamics simulations have been employed in a study of capillary fluctuations of $\left<110\right>$A, $\left<110\right>$B, and $\left<211\right>$ steps on the $(111)$ surface of face-centered cubic copper at a homologous temperature of $0.98$. The simulation results were analyzed within the framework of the statistical-mechanical theory of capillary waves. Specifically, the mean-square fluctuation amplitudes $\langle\left|A(k)\right|^{2}\rangle$ derived from the simulation data were found to follow the inverse square dependence on wave number ($k^{-2}$) predicted from capillary fluctuation theory over a range $k_{\text{min}}$ to $k_{\text{max}}$. This range is bounded by values $k_{\text{min}}$, below which the fluctuation relaxation times were too long to be adequately sampled in an MD simulation of $100\,\text{ns}$, and $k_{\text{max}}$, above which the atomic-scale wavelengths of the fluctuations lead to deviations from from the theoretical $k^{-2}$ scaling. Over this range of wave numbers, the fluctuation relaxation times were observed to display a dependence on wave number consistent with a $k^{-4}$ scaling, corresponding to kinetics limited by step-edge diffusion. We notice that although there are no theoretical restrictions to the application of the CFM to temperatures much lower than $T_{\text{m}}$, it is possible to find a limit due to computational resources because lower temperatures will require longer simulation times to sample the step fluctuations adequately.

From the measured fluctuation amplitudes we derive step stiffness values ($\tilde{\sigma}$) for each of the step orientations considered, obtaining values for the largest system sizes of $(37\pm 1)\,\text{meV}/\,\mathring{\mathrm{A}}$ that are isotropic within the statistical precision of the simulation results. The values of $\tilde{\sigma}$ derived by this CFM approach are compared to recent results for the step free energy ($\gamma^{\text{st}}$) obtained for the same system using an alternative thermodynamic-integration approach Freitas et al. (2017). The TI values of $\gamma^{\text{st}}$ and CFM values of $\tilde{\sigma}$ show discrepancies at the level of $25\%$.

We discuss that the level of discrepancy can be considered within statistical-mechanical theories for step free energies (e.g., Ref. Gelfand and Fisher (1990)) that draw a distinction between values for the “bare” stiffness $\sigma$ and the step free energy $\gamma^{\text{st}}$ that arises from configurational free energy contributions to the latter. The theoretical considerations developed in this previous literature, discussed in the context of the present results, point to an opportunity to use the CFM analysis framework described in this paper and the TI formalism in Ref. Freitas et al. (2017) to derive more detailed insights into the connection between the (possibly $k$-dependent) values of $\sigma$ that govern fluctuations at intermediate length scales and the step free energy.

Acknowledgements.

The authors would like to thank Professors Jeffrey J. Hoyt and Alain Karma for valuable discussions. The authors also gratefully acknowledge the anonymous referee for the comments and details provided about the capillary-wave theory. The research of R.F. and M.A. at UC Berkeley was supported by the U.S. National Science Foundation (Grants Nos. DMR-1105409 and 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 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 DE-AC52-07NA27344.