Effect of Surface Interatomic Potential on Thermal Accommodation

Subscriber access provided by Universitaetsbibliothek | Johann Christian Senckenberg

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Effect of Surface Interatomic Potential on Thermal Accommodation Coefficients Derived From Molecular Dynamics Timothy A. Sipkens, and Kyle J. Daun J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b06394 • Publication Date (Web): 13 Aug 2018 Downloaded from http://pubs.acs.org on August 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Effect of Surface Interatomic Potential on Thermal Accommodation Coefficients Derived from Molecular Dynamics T. A. Sipkens*, K. J. Daun Department of Mechanical and Mechatronics Engineering University of Waterloo, 200 University Ave. W., Waterloo ON Canada N2L 3G1 Corresponding author: Timothy A. Sipkens, [email protected]

Abstract This work investigates how the interatomic surface potential influences MD‐derived thermal accommodation coefficients (TACs). Iron, copper, and silicon surfaces are considered over a range of temperatures that include their melting points. Several classes of potentials are reviewed, including two‐body, three‐body, and bond‐order force fields. MD‐derived densities and visualization of the surfaces are used to explain the differences in the parameterizations of these potentials within the context of gas‐surface scattering. Finally, TACs are predicted for a range of gas‐surface combinations and recommended values of the TAC are selected that take into account the robustness and uncertainties of each of the considered parameterizations. Further, it is observed that there is a significant change in the TAC about phase changes that must be taken into account for applications with a large range of surface temperatures. 1 Introduction Free molecular conduction is important for certain gas sensors1, film coatings that enhance heat transfer2, hydrogen storage devices3, and gas‐phase nanoparticle synthesis4–6, as well as time‐resolved laser‐induced incandescence (TiRe‐ LII)7,8, which is used to characterize aerosolized nanoparticles. In the free molecular regime, the energy transfer is quantified using the thermal accommodation coefficient (TAC), , defined as



Eo  Ei

(1) where · denotes an average, and Ei and Eo are the incident and scattered energies of the gas molecule, respectively. The denominator represents the maximum energy transfer allowed by the Second Law of Thermodynamics, which would occur were the gas molecules to equilibrate with the surface. Both experimental and theoretical studies have been used to quantify this parameter for a broad range of applications over the past fifty years9. Classical treatments consist of simplified expressions for the TAC derived from first principles10. Most early studies considered only a single surface atom or approximated the surface atoms as hard spheres. The Baule model11, for example, considers a gas atom colliding with a solitary surface atom initially at rest. This model does not account for the thermal motion of surface atoms and requires the gas atom to be much less massive than the surface atom. Cube models, in which the gas molecule interacts with a moving surface cube that represents an ensemble of surface atoms, address some of these shortcomings12,13. Although these analytical models allow for a greater understanding of the physics underlying the trends in the TAC (c.f. Refs. 14,15), their simplicity significantly limits the accuracy of TAC estimates. Molecular dynamics (MD) relaxes the assumptions needed to obtain an analytically‐tractable solution. Early

Eo  Ei

max

simulations by Wachman et al.16 and Charita et al.17 used MD to evaluate the TAC but did not incorporate realistic surface interatomic potentials and were restricted in their gas‐ surface combinations. Since then, progressively more complex and realistic simulations have been performed for a broader range of gas‐surface pairs. Daun et al. considered monatomic18 and polyatomic19 gas molecules scattering from graphite. Kamat et al.20 used a more refined ReaxFF reactive force field to account for the scattering of multiple gas species and to model chemical reactions between those species and carbon nanoparticles. Several other studies21–23 used MD to model energy transferred when gases flow through nanochannels. Hu and McGaughey24 and Schiffres et al.25 performed similar simulations but within multiwalled carbon nanotubes. MD simulations have also been used to calculate TACs for gases scattering from Mo26, Fe26, Ni27,28, Si29, Au30,31, Pt32, and Al15,33 surfaces. MD‐simulations26,29 have also been used to interpret experimental TiRe‐LII data6,29,34. Uncertainties in MD‐derived TACs stem from various sources, including the uncertain boundary conditions and gas‐surface potential parameters. In response, several parametric studies have investigated the sensitivity of MD‐ derived TACs to gas temperature17,19,26, surface temperature19, and potential well depth24,26,27. One remaining uncertainty that has not been addressed is the sensitivity of the TAC to the surface potential parameterization, which often involve large sets of fitting parameters and diverse functional forms that make parametric studies challenging. This work compares the different potentials available in the literature to elucidate their differences and to investigate how uncertainties will propagate into MD‐derived TACs. Iron, copper, and silicon surfaces are studied, representing a diverse set of materials with different crystal structures and dynamics. The various classes of empirical interatomic potentials are first reviewed, including discussion of three‐body, embedded‐

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

atom model (EAM), and bond‐order potentials. Differences between the force fields are assessed through surface visualizations and calculating MD‐derived densities before proceeding with calculation of the TAC using a Monte Carlo procedure. This procedure is repeated over a range of surface temperatures, including the large temperature ranges important to TiRe‐LII. We finally give recommended values of the TAC for each of the gas‐surface pairs considered in this work. Most importantly, this comparative study provides insights into how the surface dynamics are related to the kinetics of gas‐surface scattering. 2 Empirical interatomic potentials One of the biggest challenges in any MD simulation lies in defining accurate interatomic potentials. The most ubiquitous are the Lennard‐Jones (LJ) 6‐12 potential 12 6 Uij  4D   rij     rij   ,  

(2)

and Morse potential

U ij  D exp   2 a  rij  R    2 D exp   a  rij  R   ,

(3)

where D, , R, and a are material constants. These pairwise potentials are often used to model van der Wals‐ and bonding‐type interactions, respectively. Unfortunately, such simple potentials cannot describe all possible material configurations, which establishes a need for more sophisticated treatments. 2.1

Embedded atom model force fields

The embedded atom model (EAM) class of potentials, originally attributed to Daw and Baskes35, are among the most widely used to represent metals. They operate on the principle of embedding atoms within an electron cloud, resulting in a potential energy of the form

  1  j  rij   , (4)  U 2  rij   F   2 j i  j i  where U2 is a pairwise potential between the ith and jth atoms, j is the electron charge density resulting from the jth surrounding atom, and F is the embedding function. The functional forms of U2, j, and F vary depending on the parameterization selected. The Finnis‐Sinclair (FS) potential, for example, uses an embedding function of F(x) = x1/2, based on tight‐binding theory. Modified embedded‐atom models (MEAM)36,37 modify the electron density presented above from a simple superposition of spherical densities to ones that depend on the relative position and direction of other atoms. This extends the EAM form to also include an angular dependence typical of covalent bonding arising from d‐electrons in transition metals. Ui 

2.2

Three-body potentials

Three‐body potentials are formed from the linear combination of two‐body, U2, and three‐body, U3, terms,

Ui 

Page 2 of 32

1 U 3  rij , rik , ijk  , U 2  rij    2 j i j i k  j

(5)

where ijk is the angle formed between the three atoms denoted as i, j, and k. Two‐body and three‐body terms are intended to represent bond stretching and bending, which is useful in modeling diamond cubic structures like silicon and carbon. The most common parameterization of this form is the Stillinger‐Weber (SW) potential38, which adopts a form similar to the LJ potential with an explicit cutoff function applied to the pairwise contribution. 2.3

Bond-order potentials

The bond‐order class of potentials, derived from tight‐ binding theory, are multibody potentials that model atomic interactions as bonds that are modified by their local environment, expressed mathematically as U ij  U R  rij   bijkU A  rij  ,

(6)

where UR and UA are the repulsive and attractive components of the potential representing the bond and bijk modifies the attractive component based on the environment. This treatment is well‐suited to covalently‐ bonded materials, like silicon and carbon, where interactions are dominated by interactions between closest neighbours. The form of bijk varies greatly depending on the parameterization. The aforementioned FS potential, for example, can be represented in this form39, where bijk is a function of the local electron density. The Tersoff potential40 is one of the simplest bond‐order potentials, having repulsive and attractive components similar to the Morse potential and an explicit cutoff function. The cutoff function is much sharper compared to those of other potentials, a characteristic noted later in this work. The bijk term in this case incorporates three‐body interactions, specifying an angular dependence that acts to weaken the bonds for unfavorable configurations. The environment‐dependent interatomic potential (EDIP)41 has a similar to the SW potential, but modifies the attractive term based on the coordination number, Zeff, that is the number of first neighbours to an atom. At low coordination numbers, the EDIP potential can reinforce bonding and predict sparse diamond crystal structures. However, as the coordination number increases and the atoms become more tightly packed, the EDIP potential allows for weaker angular dependences that favor metal‐like interactions (like liquid silicon) or sp2 hybridization (like that in graphite). This treatment differs from the Tersoff potential in that it considers interactions on an atom‐by‐ atom basis, without requiring symmetry. 3 MD simulations of gas-surface scattering Consider now MD simulations of gas scattering from iron, copper, and silicon surfaces. This work follows the procedure of Ref. 29, using LAMMPS42 with a timestep of 1 fs. Scattering is simulated over an infinite surface generated by applying periodic boundary conditions to the lateral edges of the simulation domain. Surfaces are first warmed to a specified temperature using the Nose‐Hoover thermostat for 30 ps. The warmed surfaces are then allowed to continue


Page 3 of 32

in the NVE ensemble for 5 ps, to ensure that steady state is reached before the atomic trajectories are written to a dump file for use in subsequent gas scattering simulations. Videos demonstrating the first 3.5 ps of this warming procedure for iron and silicon surfaces are included in the supplemental material. The scattering simulations consider gas atoms issuing from a Maxwellian gas. The TAC is defined as an average over all possible incident gas velocities and surface states,       v i ,    f v  v i ; Tg   f    ; Ts   d v i  d  , (7)  vi

where  represents the surface state, consisting of the instantaneous position and velocity of the surface atoms; v = [vn,vt]T is the gas velocity with normal, vn, and tangential, vt, components; vi is the initial velocity with corresponding normal and tangential components; (v,) is the value of the TAC for a single scattering event, which, for a monatomic gas is14





mv v o 2  v i 2 Eo  Ei   vi ,     ; Eo  Ei max 4k B Ts  Tg  2

2

tangential velocity components into Eq. (8) and by reducing the maximum energy accordingly14. The total TAC is the average of these modes,  = (n + t)/2. 3.1

Iron surfaces

Most of the available surface potentials for iron belong to the EAM class. Comparative studies are limited to evaluation of melting characteristics during laser heating53 and the presence of interstitials and defects54,55. Two significant conclusions of these studies are that Morse‐type potentials are generally incapable of describing the properties of iron and that the melting point and thermal expansion of iron are well‐described by the FS parameterized by Mendelev et al.46. In this work, we examine four parameterizations for iron, including those of Filippova et al.43 (LJ), Zhou et al.47 (EAM), Mendelev et al.46 (FS), and Müller et al.48 (Tersoff). Figure 2 shows the pairwise potential (or effective pairwise potential for the EAM and FS cases56,57) for the parametrizations considered in this work, along with several

(8)

0.2

fv(v;Tg) is the Maxwellian distribution for the velocity at a specified gas temperature, Tg.; and f(;Ts) is the distribution of surface states for a specified surface temperature, Ts. The integral in Eq. (7) is evaluated using a Monte Carlo procedure

0.1

1 N     k  vik ,  k  , N k 1

(9)

where k denotes the kth sample of any property and N is the number of simulated scattering events. To improve the randomness in the surface state, the TAC is calculated using six independently‐warmed surfaces for a single surface temperature. In the present case, 1500 samples, 250 per warmed surface, are used to generate estimates. A sample trajectory for Fe‐Ar at Ts = 2500 K following this procedure is shown in Figure 1. A video corresponding to this case is included in the supplemental material. The TAC is occasionally presented in terms of normal, n, and tangential, t, modes to facilitate discussion. These modes are obtained by substituting the normal and

M orse

Fe-Fe

Girifalco and Weizer

0 -0.1

U, U2eff [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-0.2

EAM Zhou et al.

-0.3

FS Mendelev et al.

-0.4 EAM

-0.5 -0.6 -0.7

LJ

-0.8 1.5

Pasianot et al.

Tersoff Müller et al.

FS Ackland et al.

FS

Filippova et al.

2.0

Finnis and Sinclair

2.5

3.0 3.5 r [Å]

4.0

4.5

5.0

Figure 2: Pair potentials for Fe‐Fe interactions from Filippova et al.43, Finnis and Sinclair44,45, Mendelev et al.46, Zhou et al.47, Müller et al.48, Pasianot et al.49, Ackland et al.50, and Girifalco and Weizer51. The effective pair scheme is used for the EAM and FS potentials. For the Tersoff potential, bijk is calculated considering the equilibrium lattice positions of a BCC structure. Faded lines correspond to potentials for which TACs were not calculated in this work.

Figure 1: A MD‐predicted trajectory of a helium atom scattering from an iron surface, for Ts = 2500 K, Tg = 300 K, and using a surface potential parameterization of the EAM potential from Zhou et al.47. This specific trajectory involved a single bounce on the surface followed by backscattering. A video demonstrating this scattering event is included in the supplemental material. Visualizations here and throughout the work are rendered in Ovito52.



others available in the literature44,45,50,51. For visualization purposes, the value of bijk for the Müller parameterization is set to unity. The shape of the Müller pair interaction is highly influenced by the cutoff function explicitly defined in the Tersoff potential and is the reason for the abrupt change in the slope of the potential around r = 3 Å. The Zhou and Mendelev parameterizations feature shallower potential wells and share similar equilibrium radii and well depths. The Müller and Filippova force fields have significantly deeper wells. In the latter case, the deeper potential well is likely required to maintain the solid phase in the absence of manybody effects. Figure 3 shows iron surfaces produced by each of the force fields at four surface temperatures (corresponding to , , , and liquid iron). The surface phases are generally well‐predicted with body‐centered cubic (BCC) structures for ‐ and ‐Fe, face‐centered cubic structures (FCC) for ‐ Fe, and no long‐range order for liquid iron. On the other hand, ‐ and ‐Fe surfaces obtained from the Filippova parameterization show noticeable strain, which is expected to disrupt phonon transport through the surface. This is not unexpected as the absence of multibody effects makes the parameterization better suited to closest‐packed structures, like FCC (such as the ‐Fe case, where this strain is not observed) and liquid structures. This parameterization also predicts a solid structure well above the melting temperature, likely a consequence of the larger potential well‐depth relative to the other parameterizations in the literature. This is consistent with the results of a similar Morse parameterization51, noted in Ref. 53. The Mendelev

Page 4 of 32

parameterization, in contrast, underpredicts the melting temperature, with no long‐range order observed in the surface realization at Ts = 1750 K. In spite of this, the melting temperature is still predicted within 200 K, consistent with Ref. 55 and presents an improvement over the original FS potential44,45 which greatly overpredicts the melting point58. In most solid cases, the expected long‐range order is disrupted at the surface, where the energy will be different than in the bulk. MD‐derived densities, m, for iron surfaces are plotted against experimental values in Figure 4. The MD‐derived densities are generally consistent with experimental values through the ‐Fe phase. Notable exceptions are MD‐derived densities for the Filippova and Müller force fields above 700 K, which do not exhibit the expected thermal expansion, likely due to their deeper potential wells. For the ‐Fe phase, the MD‐predicted densities across the parameterizations diverge significantly. The Filippova parameterization shows a contraction due to the change to a closest‐packed structure combined with the smaller equilibrium radius and deeper potential well (cf. Figure 2). The Mendelev and Zhou parameterizations are consistent with each other and predict volumetric expansion, likely due to the shallower potential well and multibody effects. The Müller parameterization is the most consistent with experimental measurements but does not exhibit any thermal expansion through the phase (again, likely due to the deeper potential well). The ‐Fe phase is generally poorly represented, with the Filippova and Müller parameterizations overpredicting the density and the Mendelev and Zhou potentials

Figure 3: Realizations of iron surfaces using a range of surface potentials, including the Filippova (LJ)43, Zhou (EAM)47, Mendelev (FS)46, and Müller (Tersoff)48 potentials. Surface temperatures are chosen so that all four phases (, , , and liquid) of iron should be realized. Images show a lateral surface, with the top surface used for scattering gas atoms.


Page 5 of 32

The Journal of Physical Chemistry 8500

M D / LJ

Fe

Filippova et al.

8000

M D / Tersoff Müller et al.

7500

Basinski et al. Mills

7000

M D / FS Mendelev et al.

M D / EAM

6500

Zhou et al.

6000 5500 0

BCC (α)

500

FCC (γ)

1000

Hixson et al.

BCC (δ)

ρ m [kg/m 3]

Liquid

1500 Ts [K]

2000

2500

3000

Figure 4: Trends in MD‐derived density of iron surfaces for the FIlippova (LJ)43, Zhou (EAM)47, Mendelev (FS)46, and Müller (Tersoff)48 potentials. Also shown are experimentally‐derived densities from Mills60, Basinski et al.61, and Hixson et al.62.

underpredicting the density. In many cases, this results in lattice distortion, which is visible in Figure 3. Liquid densities are well‐represented by the Mendelev and Zhou force fields, which correctly predict both the value of the density and the thermal expansion of the surface up to Ts = 3000 K. The other potentials predict significantly higher densities that are inconsistent with experiments. Scattering simulations are carried out using the pairwise Morse gas‐surface potential parameters from Refs. 26, 59, which are summarized in Table 1. They are derived from DFT calculations of the ground state energy for a system consisting of a single gas atom above a supercell of surface atoms. Videos of two sample trajectories of helium scattering from iron are included in the supplemental material. Figure 5 shows trends in the normal and tangential components of the MD‐derived TAC for Fe‐Ne scattering. In all cases, error bounds expand as the surface temperature approaches the gas temperature (Tg = 300 K), a consequence α

0.4

Norm al 0.3

Fe-Ne

γ

of amplification of stochastic and model errors in Eo – Ei as the denominator of Eq. (8) approaches zero. In general, energy is accommodated less efficiently into the tangential mode. For the Mendelev, Zhou, and Müller force fields, sharp increases are observed in the tangential mode about the ‐ Fe, ‐Fe, or liquid Fe transition (with the precise location of the change depending on the parameterization, e.g. the Mendelev parameterization47 does not predict a solid ‐Fe, resulting in an earlier rise in the TAC). Figure 3 indicates that this shift corresponds to a change from a relatively smooth solid surface to a liquid surface having considerable superatomic roughness (roughness on the order of several surface atoms). This allows the gas atom to interact with the lateral edges of surface atoms and increases the number of bounces, collectively resulting in increased tangential thermal accommodation. At surface temperatures below this transition, the tangential mode is effectively frozen. The Mendelev, Zhou, and Müller force fields also feature a subtler increase in the normal TAC about the ‐Fe to ‐Fe transition. This is a result of changes to the phonon dispersion in the surface following the transition from a BCC to FCC crystal structure. Between these phase transitions, both modes of the TAC exhibit only modest increases with surface temperature, most notably in the tangential component for the liquid phase. Remaining differences in the TACs predicted using the different force fields are related to the respective surface roughness scales, with the sharper surface features observed in Figure 3 for the Mendelev Table 1: Parameterizations of the Morse potential for iron‐gas pairs. The Fe‐He and Fe‐Ar parameterizations are taken directly from Ref. 26. D (meV)

a (Å‐1)

R (Å)

Fe‐He Fe‐Ne Fe‐Ar

2.483 1.860 2.238

1.290 1.340 1.204

4.281 4.160 4.779

Liquid

δ

α

0.4

FS

γ

Liquid

δ

FS

Tangent ial

Mendelev et al.

0.3

Mendelev et al.

Fe-Ne

EAM

0.2

0.2

0.1

0.1

αt

αn

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0

Zhou et al.

0 LJ Filippova et al.

-0.1

EAM

Tersoff

Zhou et al.

-0.2 500

1000

Müller et al.

1500 2000 Ts [K]

2500

3000

Tersoff

-0.1 -0.2 500

Müller et al.

LJ Filippova et al.

1000

(a)

1500 2000 Ts [K]

2500

3000

(b)

Figure 5: Trends in the (a) normal and (b) tangential TACs with Ts for Tg = 300 K and Fe‐Ne. Surface potential parameterizations include those of Filippova (LJ)43, Zhou (EAM)47, Mendelev (FS)46, and Müller (Tersoff)48. Both modes of the TAC are plotted on the same vertical scale for comparison. Error bounds are included for the Mendelev (red, dashed lines) and Filippova (yellow, shaded region) parameterizations and correspond to two standard deviations of the mean of six surface realizations. Error bounds for other cases are excluded for clarity but are similar in magnitude.


The Journal of Physical Chemistry α

Fe-He

γ

Liquid

δ

Mendelev et al.

0.15

0.30

FS

FS

Finnis and Sinclair

α

Fe-Ne

0.25 0.20

0.10

γ

Liquid

δ

0.35

FS

0.30

Mendelev et al.

EAM

α 0.10

0

EAM Zhou et al.

LJ

Tersoff

Filippov et al.

-0.05 500

Müller et al.

-0.05 1000 1500 2000 2500 3000 500 Ts [K]

Liquid

δ

FS Finnis and Sinclair

LJ Filippov et al.

FS EAM Zhou et al.

0.10

0.05 0

0.15

γ

Mendelev et al.

0.20

0.15

0.05

α

Fe-Ar

0.25

Zhou et al.

α

0.20

α

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

0.05 0 -0.05 500 1000 1500 2000 2500 3000 Ts [K] Tersoff

LJ

Müller et al.

Filippov et al.

(a)

Tersoff Müller et al.

1000 1500 2000 2500 3000 Ts [K]

(b)

(c)

Figure 6: Trends in the total TAC with surface temperature, Ts, for (a) Fe‐He, (b) Fe‐Ne, and(c) Fe‐Ar at Tg = 300 K. Surface potentials include the Filippova (LJ)43, Zhou (EAM)47, Mendelev (FS)46, and Müller (Tersoff)48 parameterizations. Also included are results using the original Finnis‐Sinclair (FS) parameterization for iron44 used by Daun et al.26. TACs are plotted on different vertical scales for each gas. Error bounds are included for the Mendelev (red, dashed lines) and Filippova (yellow, shaded region) parameterizations and correspond to two standard deviations of the mean of the six surface realizations. Error bounds for other cases are excluded for clarity but are similar in magnitude.

parameterization resulting in higher thermal accommodation and the smooth, rolling surface features of the Müller parameterization giving lower thermal accommodation. The TAC is also inversely correlated with the density of the surface and the surface potential well depth for both modes. The TACs derived using the Filippova parameterization stand out from those obtained using the other surface potentials. The TACs are significantly lower, particularly for Ts > 750 K, and the general trends with Ts are also different. This is consistent with the surface realizations and densities shown in Figure 3 and Figure 4, which also deviate significantly from the other parameterizations and experiments. One can thus conclude that this parameterization performs poorly over large temperature ranges, which is unsurprising since it was derived using the properties of iron at room temperature. These observations, along with the inadequacy of the Morse parameterization51 reported by Imamova et al.53, argue against using pairwise potentials to model iron surfaces. Figure 6 shows trends in the MD‐predicted total TAC with Ts for Fe‐He, Fe‐Ne, and Fe‐Ar. The TACs for Fe‐Ne and Fe‐Ar share many characteristics with the normal and tangential modes shown in Figure 5, including noticeable discontinuities about the phase changes and moderate increases with surface temperature through the liquid phase. The Filippova parameterization significantly underpredicts the TAC for Ts > 750 K. The TAC for Fe‐Ar is larger than the one for Fe‐Ne, following trends expected with mass ratio,  = mg/ms, where mg and ms are the molecular mass of the gas molecule and atomic mass of the surface atom, respectively7,10,34. The total TACs obtained with different parameterizations are inversely correlated with the scale of surface roughness. The Fe‐He case differs from the Fe‐Ne and Fe‐Ar results in several ways, most significantly in that the TACs are insensitive to surface temperature and do not change about the phase transitions. The Fe‐He TACs for Ts < 1000 K are also larger than those for the heavier gas molecules, contrary

to the trends expected with increasing . Much of this is due to the tangential mode, which is more accessible for Fe‐He interactions relative to the larger molecules. This is a consequence of the small diameter of helium atoms, which allows them to see the surface corrugation resulting from individual atoms in the surface, even at lower temperatures. This is particularly the case for the sparser BCC phase over the closest‐packed FCC and liquid phases. This atomic‐scale corrugation changes very little with surface temperature, resulting in the insensitivity of the TAC to surface temperature noted in Figure 6. Recommended values of the TAC at Ts = 750 K (solid) and 2500 K (liquid) are provided in Table 2. The TACs predicted by the Zhou potential47, which are bounded by the Mendelev and Müller potentials and seemingly predict the most reliable phase transitions, are taken as the most Table 2: Recommended values of the total TACs for various gases above iron, copper, and silicon for Tg = 300 K and for Ts = 750 K and liquid (Ts = 2500 K for iron and silicon and Ts = 2000 K for copper). Error bounds correspond to two standard deviations on the mean. Values for iron are a combination of predictions from Zhou (EAM)47, Mendelev (FS)46, and Müller (Tersoff)48 potentials, with the latter two potentials weighted half as heavily. The results for the original FS potential44,45 calculated by Daun et al.26 are also included in Fe‐He and Fe‐Ar values for the liquid case, also weighted half as heavily. Values for copper are the average of Zhou (EAM)63 and Foiles (EAM)64 potentials. Solid values for silicon are the average of the SW38, T2 (Tersoff)40, Justo (EDIP)65, and Jelinek (MEAM)66 potentials and the liquid values are the average of the SW, Justo, and Jelinek potentials.

 (Solid, Ts = 750 K)

 (Liquid)

Fe‐He Fe‐Ne Fe‐Ar Cu‐He Cu‐Ne Cu‐Ar Si‐He Si‐Ar

0.095 ±0.021 0.070 ±0.020 0.031 ±0.014 0.036 ±0.015 0.092 ±0.024 0.250 ±0.024 0.089 ±0.011 0.105 ±0.012

0.095 ±0.008 0.168 ±0.016 0.181 ±0.014 0.077 ±0.011 0.229 ±0.024 0.493 ±0.037 0.110 ±0.008 0.382 ±0.021


reliable in the context of predicting the TAC. The TACs from the remaining parameterizations are weighted half as significantly. The results for the original FS parameterization44,45, as calculated by Daun et al.28, are also included in the Ts = 2500 K case, again weighted half as heavily as the predictions from the Zhou parameterization. 3.2

Copper surfaces

A large number of candidate parameterizations have been proposed for copper, mostly of the EAM form. Previous comparative studies are based on predictions of high angle grain boundaries67 and some other material properties57. These studies highlight the need to include multibody effects when modeling the interactions between copper atoms. Figure 7 shows the LJ parameterization of Filippova et al.43 and the effective pair potentials for a range of EAM parameterizations50,63,64,68. Other EAM parameterizations not shown69,70 are bounded by the Universal 4 and Foiles potentials. The EAM parameterizations feature similar equilibrium distances and potential well depths. A comparison of Figure 2 and Figure 7 reveals several similarities. For instance, in both cases, the LJ parameterization has a much deeper potential well, while those of the Morse parameterizations lie between the LJ and effective EAM potentials and have equilibrium distances that are the greatest of the force fields considered. We here consider two EAM parameterizations by Zhou et al.63 and Foiles64, and the LJ parameterization of Filippova et al.43. The densities predicted by the EAM parameterizations, shown in Figure 8, are consistent with experimental values over a wide temperature range and with observations by Ryu and Cai74 for a range of similar parameterizations. The surface realizations shown in Figure 9 highlight the solid and liquid surface structures at 500 K and 1500 K, respectively. Without multiple phase changes, Figure 9 indicates that the force fields tend to be more

consistent with each other and do not predict significant strain. The Filippova parameterization presents an exception, predicting lower densities over much of the domain. Figure 9 also indicates that the Filippova (LJ) parameterization for copper predicts a crystalline structure well above the melting point, as was the case for iron. However, unlike iron, the Filippova parameterization for copper also results in a considerable number of free atoms that evaporate from the surface for Ts > 1000 K. These free atoms are likely to disrupt the initial placement of the gas atom and will interfere with the incident gas trajectories. Gas‐surface potentials for Cu‐He, Cu‐Ne, and Cu‐Ar are derived by fitting a pairwise Morse potential to the physisorption potential given by Chizmeshya and Zaremba75, using a static FCC lattice of copper atoms with a 9000

CRC Handbook

Cu Kurochkin et al.

8500

Brillo and Ergy

8000 M D / EAM Zhou et al.

7500

M D / LJ Filippova et al.

7000 M D / EAM Foiles

6500 6000 0

FCC Liquid

500

1000

1500 Ts [K]

2000

2500

3000

Figure 8: Density of copper as predicted by MD simulations using parameterization of the potential by Zhou et al.63 and Foiles64 and the LJ potential parameterized by Filippova et al.43. Also shown are experimental densities given by the CRC Handbook71, Kurochkin et al.72, and Brillo and Ergy73.

0.1

Cu-Cu 0 FS

-0.1

U, U2eff [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ρ m [kg/m 3]

Page 7 of 32

LAMMPS Mendelev

FS Ackland et al.

-0.2

EAM Universal 4 Adams et al.

-0.3

EAM

M orse

-0.4

Foiles

Wolf

EAM Zhou et al.

M orse

-0.5 -0.6 1.5

Gades and Urbassek

2.0

LJ Filippova et al.

2.5

3.0 3.5 r [Å]

4.0

4.5

5.0

Figure 7: Potentials describing Cu‐Cu interactions from Zhou et al.63, Foiles64, Ackland et al.50, Adams et al.68, LAMMPS42, Filippova et al.43, Wolf67, and Gades and Urbassek57. The Universal 369 and Mishin70 EAM potentials are generally bounded by the potential by Foiles and the Universal 4 potentials. Faded lines correspond to potentials for which TACs were not calculated in this work.

Figure 9: Realizations of copper surfaces using the Filippova (LJ) parameterization43, and the EAM potentials from Zhou et al.63 and Foiles64. Images show a lateral surface, with the top surface used for scattering gas atoms. For Ts = 1500 K and the LJ potential, considerable evaporation means that numerous atoms exist outside of the illustrated domain.



lattice constant of a0 = 3.615 Å. The Morse parameterization for each gas is provided in Table 3 and plotted as a function of interatomic distance for Cu‐He, Cu‐Ne, and Cu‐Ar in Figure 10. A video of a sample trajectory of argon scattering from copper is included in the supplemental material. Figure 11 shows the total TACs predicted for Cu‐He, Cu‐ Ne, and Cu‐Ar as a function of surface temperature for the chosen parameterizations. The Cu‐He case reveals similar features to the Fe‐He case. The modest increase in the TAC with temperature comes as a result of increases in atomic‐ scale corrugation as it is perturbed from its denser FCC Table 3: Parameterizations of the Morse potential for copper‐gas pairs. Values represent fits of the Morse form to the calculations of Chizmeshya and Zaremba75.

D (meV)

a (Å‐1)

R (Å)

Cu‐He Cu‐Ne Cu‐Ar

0.5255 1.5769 6.0557

1.1033 1.1082 1.0779

3.9372 3.4947 3.3464

structure (relative to the sparser BCC structure for iron). The TACs predicted using the Filippova parameterization change significantly with the onset of evaporation from the surface. As noted previously, the apparent increase in the TAC as the surface temperature approaches the gas temperature (300 K) is a numerical artifact resulting from amplification of stochastic and modeling errors in the numerator of Eq. (8) as the denominator approaches zero. For Cu‐Ne and Cu‐Ar, like the Fe‐Ne and Fe‐Ar cases above, the TAC changes abruptly about the melting point. This is again attributed to a change in the superatomic roughness of the surface upon melting (which can also be seen in Figure 9). TACs predicted by both EAM parameterizations are consistent across the entire temperature range. The Filippova parameterization suffers from many of the failings noted both here and for iron in the preceding section. The TACs for Cu‐Ne and Cu‐Ar follow the expected trends with mass ratio10,14,29. Recommended values are presented in Table 2. 3.3

8

Cu-Gas

Cu-He Zaremba and Kohn (ZK)

6

Cu-He

U [mev]

4

Chizm eshya and Zarem ba

Cu-Ne

2

Chizm eshya and Zarem ba

0 -2 -4 Cu-Ar Chizmeshya and Zarem ba

-6 -8 2.0

2.5

3.0

3.5

4.0 r [Å]

4.5

5.0

5.5

6.0

Figure 10: Inferred Morse potential for Cu‐gas interactions based on physisorption potentials given in Chizmeshya and Zaremba75. Also shown is a Morse potential derived using the same procedure, but from the physisorption potential given by Zaremba and Kohn76.

Zhou et al.

0.10

0.30

Cu-He

EAM EAM

α Filippova et al.

1000

Ts [K] (a)

1500

2000

Zhou et al.

0.4

EAM Foiles

0.15

0 500

Filippova et al.

Ts [K]

1500

(b)

EAM Zhou et al.

EAM Foiles

0.3 0.2

LJ

1000

FCC Liquid

Cu-Ar

0.5

0.05

LJ

0.6

EAM

0.10

0 -0.05 500

FCC Liquid

0.20

0.05

Many studies have compared the myriad of surface potentials available for silicon74,77–82, often based on how well the potentials predict the melting point (which can differ from the experimental value by 1000 K). In this work, we focus on: (i) the original Stillinger‐Weber (SW) (ii) the Tersoff 2 (T2) parameterization38; parameterization40; (iii) the EDIP parameterization by Justo et al.65; and (iv) the MEAM parameterization by Jelinek et al.66. Figure 12 shows the two‐body and angular components of these parameterizations, along with several others available in the literature; the angular functions are normalized by their maximum. For the Justo parameterization, Figure 12 shows both the ideal liquid (Zeff = 12) and diamond cubic (Zeff = 4) cases and indicates how transitioning to higher coordination numbers shifts the potential from enforcing bonding at specific angles to capturing only repulsive interactions. Also note that the two‐ body component of the T2 and T3 force fields are very similar, while the angular components are quite different.

Cu-Ne

0.25

Foiles

Silicon surfaces

α

FCC Liquid

0.15

α

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

2000

LJ Filippova et al.

0.1 0 500

1000

Ts [K]

1500

2000

(c)

Figure 11: Trends in MD‐derived total TACs for (a) Cu‐He, (b) Cu‐Ne, and (c) Cu‐Ar using two EAM parameterizations63,64 and the LJ parameterization by Filippova et al.43. Error bounds are included for the Zhou (purple, dashed lines) and Filippova (yellow, shaded region) parameterizations and correspond to two standard deviations of the mean of scattering from six realizations of the copper surface.


Page 9 of 32 2

1.0

Si-Si

ED IP Justo et al. Zeff = 12

Tersoff Erhart and Albe

0.8

ED IP Justo et al. Zeff = 4

SW

0 -1

ED IP

Tersoff

Justo et al. Zeff = 4

Tersoff 3

-2

geff (θ)

U2, Ueff 2 [eV]

1

Stillinger and Weber

0.6

Tersoff

0.4

Tersoff

0.2

ED IP

Tersoff 3

Tersoff 2

Tersoff SW

Tersoff 2

-3


Tersoff Erhart and Albe

-4 1.0

1.5

2.0

2.5 r [Å] (a)

3.0

3.5

4.0

0 0


30

60

90 θ [°] (b)

120

150

180

Figure 12: Pair potentials for Si‐Si interactions from Tersoff (T240 and T384), Stillinger and Weber (SW)38, Justo et al. (EDIP)65, and Erhart and Albe (Tersoff)85. For the Tersoff potentials, bijk = 1. For the Justo potential, two effective coordination numbers are shown: Zeff = 4 corresponding to an ideal diamond cubic phase and Zeff = 12 corresponding to an ideal liquid phase. The angular functions, geff(), are scaled against their maximum (except for the Zeff = 12 case of the Justo potential which is scaled by the maximum of the Zeff =4 case). Faded lines correspond to potentials for which TACs were not calculated in this work.

Figure 13 shows the MD‐ and experimentally‐derived densities of silicon as a function of surface temperature, and Figure 14 shows realizations of solid (Ts = 500 K) and liquid silicon surfaces (Ts = 1500 K). Below the melting point, the SW, T2, and Justo (EDIP) parameterizations predict very similar densities. The slight thermal contraction visible in the Justo case is the result of several atomic layers adjacent to the top and bottom surface collapsing away from the sparse diamond crystal structure in the bulk. This is likely a consequence of the drop in the coordination number at the surface and also corresponds to an increase in superatomic roughness, which affects TAC calculation. The Jelinek (MEAM) parameterization predicts thermal expansion, which Figure 14 indicates is a result of high surface energies and an early onset of melting at the surface. We note that the surface roughness appears on a smaller scale for the Jelinek 2600

Si

2500 2400

ρ m [kg/m 3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Rhim et al.

M D / SW

Oshaka et al.


Room tem perature density

2300 2200

M D / M EAM Jelinek et al.

M D/ Tersoff

2100

Tersoff 3

2000 1900 0

M D/ ED IP Justo et al.

c-Si Liquid

500

1000

1500 Ts [K]

2000

2500

3000

Figure 13: Trends in MD‐derived density of silicon surfaces for the SW38, T2 (Tersoff)40, Justo (EDIP)65, and Jelinek (MEAM)66 potentials. Error bars correspond to one standard deviation of the density of six realizations of a silicon surface. Also shown are experimental densities for liquid silicon from Rhim et al.86 and Oshaka et al.87, the room temperature density from Ref. 71, and the solid density at the melting point from Ref. 88,89.

parameterization (where it appears to be almost atomic) than the Justo parameterization (where the roughness is on the scale of the simulation box). The densities predicted by the various force fields diverge significantly at temperatures above the melting point. The SW parameterization predicts the melting point accurately, consistent with previous observations78,82,83 and is also the only parameterization to correctly predict that the density increases abruptly upon melting89. In contrast, the density predicted by the T2 potential is consistent with solid silicon at temperatures well above the melting point. This phenomenon has also been reported in the literature78,90 and is attributed to the strength of the bond bending term (partially manifested in the deeper potential well in Figure 12). This is further affirmed by Figure 14, which shows that the T2 potential predicts a solid structure at Ts = 2500 K. The Justo (EDIP) parametrization appears to predict a melting point within a couple hundred Kelvin of the experimental value. Unlike the SW case, however, the Justo parameterization predicts a relatively sparse structure, similar to what would be expected of amorphous silicon (e.g. Ref. 91,92). This expansion is likely a result of the bonds weakening as the coordination number increases for liquid silicon, allowing the material to expand until lower coordination numbers are achieved. Figure 14 indicates that these effects result in rougher surfaces than those predicted by the SW parameterization at Ts = 2500 K. The Jelinek (MEAM) parameterization also exhibits thermal expansion and underpredicts the melting temperature (around Ts ≈ 1500 K), which is roughly consistent with observations for similar force fields74,78. The gas‐surface potential for the scattering simulations are obtained following the same procedure as the iron case above, using a Morse parameterization fit to DFT calculations (in this case from Ref. 29). This improves on the treatment of Zhao et al.82, who used a purely repulsive potential for the gas‐surface interaction. A video demonstrating a sample trajectory of helium scattering from silicon is included in the supplemental material.



Figure 14: Realization of silicon surfaces using the SW38, T2 (Tersoff)40, Justo (EDIP)65, and Jelinek (MEAM)66 potentials at Ts = 1500 K (where a solid, diamond cubic structure is expected) and Ts = 2500 K (where a liquid structure is expected). Images shows a lateral surface, with the top surface used for scattering gas atoms.

Figure 15 shows trends with surface temperature in the normal and tangential modes of the TAC for Si‐Ar. The low temperature limit is consistent between all four potentials and for both modes, indicating similarity between the force fields close to room temperature. The same is true of the SW, Justo (EDIP), and Jelinek (MEAM) parameterizations in the high temperature limit for the normal mode. Small differences between the parameterizations in the high temperature limit of the tangential mode are well‐correlated with the magnitude of the surface roughness observed in Figure 14, which is most significant for the Jelinek parameterization for reasons noted above. At moderate temperatures, the TACs for both modes predicted by the SW and Jelinek parameterizations are relatively constant, with a sharp transition about the MD‐predicted melting point. The transition for the Justo parameterization is more gradual, consistent with observations that there are increases in the surface roughness before the surface fully melts. The T2 c-Si

0.7 0.6

Norm al Si-Ar

0.5 0.4

Liquid

0.7

M EAM

0.6

Jelinek et al.

c-Si

M EAM

Si-Ar

0.4

0.2

Liquid

Tangent ial

Jelinek et al.

0.5

ED IP Justo et al.

0.3

0.3

ED IP SW

Justo et al.


0.2

0.1

0.1

0 -0.1 -0.2 500

parameterization predicts a nearly constant TAC for both modes over the entire temperature range, consistent with the delayed onset of melting. The TAC begins to increase when the surface does start to melt (Ts > 2500 K). Figure 16 shows how the total TAC varies with surface temperature for Si‐He and Si‐Ar. The Si‐Ar case features the expected change at the melting point. The Si‐He case, in contrast, exhibits very little change over the entire temperature range considered. As with the iron and copper cases, the difference can be attributed to the smaller size of the helium atom, which is scattered by atomic‐scale surface corrugation, which is enhanced by the sparser diamond structure of silicon. Further differences between the force fields are correlated with the degree of superatomic roughness, which Figure 14 indicates is most significant for the Jelinek parameterization.

αt

αn

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32


Tersoff 2

1000

0

SW

Tersoff

1500 2000 Ts [K]

2500

3000

-0.1 -0.2 500

Tersoff Tersoff 2

1000

(a)

1500 2000 Ts [K]

2500

3000

(b)

Figure 15: Trends in MD‐derived (a) normal and (b) tangential TACs for Si‐Ar using the SW38, T2 (Tersoff)40, Justo (EDIP)65, and Jelinek (MEAM)66 potentials. Error bounds are included for the Justo (red, long‐dashed lines), T2 (blue, short‐dashed lines), and SW (yellow, shaded region) parameterizations and correspond to two standard deviations of the mean of scattering from six realizations of the silicon surface. The vertical axis is consistent between (a) and (b).


Page 11 of 32

c-Si

0.25

Liquid

Si-He 0.20 0.15

Si-Ar

M EAM

0.4

Jelinek et al.

ED IP

Liquid

M EAM Jelinek et al.

ED IP

Justo et al.

0.3

0.10

Justo et al.

0.2 0.1

0.05

Tersoff Tersoff 2

0 -0.05 500

c-Si

0.5

α

α

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

SW Stillinger and Weber

1000

1500 2000 Ts [K]

2500

3000

-0.1 500

1000

(a)

Tersoff

SW

Tersoff 2


1500 2000 Ts [K]

2500

3000

(b)

Figure 16: Trends in MD‐derived total TACs for (a) Si‐He and (b) Si‐Ar using the SW38, T2 (Tersoff)40, Justo (EDIP)65, and Jelinek (MEAM)66 potentials. Error bounds are included for the Justo (red, long‐dashed lines), T2 (blue, short‐dashed lines), and SW (yellow, shaded region) parameterizations and correspond to two standard deviations of the mean of scattering from six realizations of the silicon surface.

4 Conclusions While molecular dynamics is a useful tool for modelling gas‐ surface scattering, the role that the surface potential parameterization plays in these predictions is rarely considered. This work addresses these concerns by calculating the thermal accommodation coefficient (TAC) for monatomic gases scattering from iron, copper, and silicon surfaces using a range of surface parametrizations. This comparative assessment also provides important physical insights into how surface dynamics influence gas‐surface scattering. Scattering is highly influenced by the surface phase and scale of surface roughness. The latter quality has a significant impact on the tangential mode of the TAC, since increasing surface roughness allows gas atoms to interact with the lateral edges of surface atoms and increases the number of times the gas atom bounces on the surface. This effect is particularly pronounced for helium atoms due to their small size relative to atomic‐scale surface corrugation, which activates the tangential mode low temperatures and desensitizes the tangential TAC to surface temperature. Heavier gas atoms are only sensitive to superatomic surface roughness, which occurs upon melting. In the case of iron surfaces, the manybody potentials predict similar trends in the TAC, which increases when the surface melts and as the surface transitions from ‐Fe to ‐ Fe. Differences in the TACs obtained from the various parameterizations are related to the predicted scale of surface roughness observed in the surfaces. The Zhou parameterization (EAM)47, which best predicts the phase transitions and experimental densities, appears to be the most robust choice, while the pairwise Filippova parameterization (LJ)43 is not recommended. For copper surfaces, the Zhou63 and Foiles64 (EAM) parameterizations predict similar TACs and both show an increase in the TAC at the melting point. As with iron, the Filippova (LJ) parameterization43 fails to reliably reproduce key attributes over wide temperature ranges and is not recommended for calculating the TAC.

Finally, for silicon surfaces, the SW38 parameterization most accurately predicts the melting point, while the T2 (Tersoff) and Jelinek (MEAM)66 parameterizations overpredict and slightly underpredict the melting point, respectively. In all cases the TAC increases upon melting. The Justo (EDIP) parameterization65, in contrast, predicts high surface energies and significant superatomic roughness even for solid surfaces, resulting in a more gradual increase in the TAC. The TAC is reasonably consistent between all of the parametrizations at low temperatures, and between the SW, Justo, and Jelinek parameterizations at high temperatures. Future work will focus on extending this treatment to carbon, including a study of its graphitic and amorphous phases. This will further solidify observations regarding the effects that phase changes in the material have on the TAC. Supporting information A video of a helium atom scattering from an iron surface at Ts = 2500 K (FeHe_scatter1.avi) A second video of a helium atom scattering from an iron surface at Ts = 2500 K (FeHe_scatter2.avi) A video of an argon atom scattering from a copper surface at Ts = 1500 K (CuAr_scatter.avi) A video of a helium atom scattering from a silicon surface at Ts = 2500 K (SiHe_scatter.avi) A video of the warming procedure performed on an iron surface at Ts = 2500 K (Fe_warm.avi) A video of the first 3.5 ps of the warming procedure performed on a silicon surface at Ts = 2500 K (Si_warm.avi) Acknowledgements This research was supported by the Natural Science and Engineering Research Council (NSERC) RGPIN‐2018‐03756. The authors would also like to acknowledge the aid of Professor Mikko Karttunen and Dr. John Titantah.



References [1] Li, J.; Lu, Y.; Ye, Q.; Cinke, M.; Han, J.; Meyyappan, M. Carbon nanotube sensors for gas and organic vapor detection. Nano Lett., 2003, 3, 929‐933. [2] Kinefuchi, I.; Shiomi, J.; Takagi, S.; Maruyama, S.; Matsumoto, Y. Gas−surface energy exchange in collisions of helium atoms with aligned single‐walled carbon nanotube arrays. J. Phys. Chem., 2013, 117, 14254‐14260. [3] Dillon, A. C.; Jones, K. M.; Bekkedahl, T. A.; Kiang, C. H.; Bethune, D. S.; Heben, M. J. Storage of hydrogen in single‐walled carbon nanotubes. Nature, 1997, 386, 377‐379. [4] Lümmen, N.; Kraska, T. Investigation of the formation of iron nanoparticles from the gas phase by molecular dynamics simulation. Nanotechnology, 2004, 15, 525‐ 533. [5] Eremin, A.; Gurentsov, E.; Schulz, C. Influence of the bath gas on the condensation of supersaturated iron atom vapour at room temperature. J. Phys. D: Appl. Phys., 2008, 41, 1‐5. [6] Eremin, A. V.; Grentsov, E. V. Sizing of Mo nanoparticles synthesized by Kr‐F laser pulse photo‐dissociation of Mo(CO)6. Appl. Phys. A, 2015, 119, 615‐622. [7] Daun, K. J.; Smallwood, G. J.; Liu, F. Investigation of thermal accommodation coefficients in time‐resolved laser‐induced incandescence. J. Heat Trans., 2008, 130, 121201. [8] Michelsen, H. A. Derivation of a temperature‐ dependent accommodation coefficient for use in modeling laser‐induced incandescence of soot. Appl. Phys. B, 2009, 94, 103‐117. [9] Saxena, S. C.; Joshi, R. K. Thermal Accommodation and Adsorption Coefficients of Gases; Touloukian, Y. S., Ho, C. Y., Eds.; McGraw‐Hill Book Company: New York, 1981. [10] Goodman, F. O.; Wachman, H. Y. Dynamics of Gas‐ Surface Scattering; Academic Press: New York, 1976. [11] Baule, B. Phenomena in rarefied gases. Ann. Physik, 1914, 44, 145‐176. [12] Logan, R. M.; Stickney, R. E. Simple classical model for the scattering of gas atoms from a solid surface. J. Chem. Phys., 1966, 44, 195‐201. [13] Logan, R. M.; Keck, J. C. Classical theory for the interaction of gas atoms with solid surfaces. J. Chem. Phys., 1968, 49, 860‐876. [14] Sipkens, T. A.; Daun, K. J. Using cube models to understand trends in thermal accommodation coefficients at high surface temperatures. Int. J. Heat Mass Transfer, 2017, 111, 54‐64. [15] Mane, T. S.; Bhat, P.; Sundaram, D. S.; Yang, V. Energy accommodation under non‐equilibrium conditions for aluminum‐inert gas systems. Surf. Sci., 2018, 677, 135‐ 148. [16] Wachman, H. Y.; Greber, I.; Kass, G. Molecular dynamics computations of scattering from a surface using a Lennard‐Jones model of a solid. In Progress in Astronautics and Aeronautics; Vancouver, 1992.

Page 12 of 32

[17] Chirita, V.; Pailthorpe, B.; Collins, R. Molecular dynamics study of low‐energy Ar scattering by the Ni (001) surface. J. Phys. D: Appl. Phys., 1993, 26, 133‐142. [18] Daun, K. J.; Smallwood, G. J.; Liu, F. Molecular dynamics simulations of translational thermal accommodation coefficients for time‐resolved LII. Appl. Phys. B, 2009, 94, 39‐49. [19] Daun, K. J. Thermal accommodation coefficients between polyatomic gas molecules and soot in laser‐ induced incandescence experiments. Int. J. Heat Mass Transfer, 2009, 52, 5081‐5089. [20] Kamat, A. M.; van Duin, A. C.; Yakovlev, A. Molecular dynamics simulations of laser‐induced incandescence of soot using an extended ReaxFF reactive force field. J. Phys. Chem. A, 2010, 114, 12561‐12572. [21] Sun, J.; Li, Z. X. Molecular dynamics simulations of energy accommodation coefficients for gas flows in nano‐channels. Mol. Simul., 2009, 35, 228‐233. [22] Prabha, S. K.; Sathian, S. P. Molecular‐dynamics study of Poiseuille flow in a nanochannel and calculation of energy and momentum accommodation coefficients. Phys. Rev. E, 2012, 85, 041201. [23] Spijker, P.; Markvoort, A. J.; Nedea, S. V.; Hilbers, P. A. Computation of accommodation coefficients and the use of velocity correlation profiles in molecular dynamics simulations. Phys. Rev. E, 2010, 81, 011203. [24] Hu, L.; McGaughey, A. J. H. Energy accommodation between noble gases and carbon nanotubes. J. Phys. Chem. C, 2013, 117, 18804‐18808. [25] Schiffres, S. N.; Kim, K. H.; Hu, L.; McGaughey, A. J.; Islam, M. F.; Malen, J. A. Gas diffusion, energy transport, and thermal accommodation in single‐walled carbon nanotube aerogels. Adv. Funct. Mater., 2012, 22, 5251‐ 5258. [26] Daun, K. J.; Sipkens, T. A.; Titantah, J. T.; Karttunen, M. Thermal accommodation coefficients for laser‐induced incandescence sizing of metal nanoparticles in monatomic gases. Appl. Phys. B, 2013, 112, 409‐420. [27] Daun, K. J.; Titantah, J. T.; Karttunen, M. Molecular dynamics simulation of thermal accommodation coefficients for laser‐induced incandescence sizing of nickel particles. Appl. Phys. B, 2012, 107, 221‐228. [28] Daun, K. J.; Titantah, J. T.; Karttunen, M. Erratum to: Molecular dynamics simulation of thermal accommodation coefficients for laser‐induced incandescence sizing of nickel particles. App. Phys. B, 2013, 112, 599‐600. [29] Sipkens, T. A.; Mansmann, R.; Daun, K. J.; Petermann, N.; Titantah, J. T.; Karttunen, M.; Wiggers, H.; Dreier, T.; Schulz C. In situ nanoparticle size measurements of gas‐borne silicon nanoparticles by time‐resolved laser‐ induced incandescence. Appl. Phys. B, 2014, 116, 623‐ 636. [30] Liang, Z.; Keblinski, P. Parametric studies of the thermal and momentum accommodation. Int. J. Heat Mass Tran., 2014, 78, 161‐169. [31] Giri, A.; Hopkins, P. E. Analytical model for thermal boundary conductance and equilibrium thermal


Page 13 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


accommodation coefficient at solid/gas interfaces. J. Chem. Phys., 2016, 144, 084705. [32] Yamaguchi, H.; Matsuda, Y.; Niimi, T. Molecular‐ dynamics study on characteristics of energy and tangential momentum accommodation coefficients. Phys. Rev. E, 2017, 96, 013116. [33] Sha, H.; Tetiker, G.; Woytowitz, P.; Faller, R. Molecular simulation study of aluminum–noble gas interfacial thermal accommodation coefficients. AIChE J., 2018, 64, 338‐345. [34] Sipkens, T. A.; Singh, N. R.; Daun, K. J. Time‐resolved laser‐induced incandescence characterization of metal nanoparticles. Appl. Phys. B, 2017, 123, 1, 14‐30. [35] Daw, M. S.; Baskes, M. I. Embedded‐atom method: Derivation and application to impurities, surfaces, and other defects in metals. Phys. Rev. B, 1984, 29, 6443‐ 6453. [36] Baskes, M. I.; Nelson, J. S.; Wright, A. F. Semiempirical modified embedded‐atom potentials for silicon and germanium. Phys. Rev. B, 1989, 40, 6085‐6100. [37] Baskes, M. I. Modified embedded‐atom potentials for cubic materials and impurities. Phys. Rev. B, 1992, 46, 2727‐2742. [38] Stillinger, F. H.; Weber, T. A. Computer simulation of local order in condensed phases of silicon. Phys. Rev. B, 1985, 31, 5262‐5271. [39] Brenner, D. W. Relationship between the embedded‐ atom method and Tersoff potentials. Phys. Rev. Lett., 1989, 63, 1022‐1022. [40] Tersoff, J. New empirical approach for the structure and energy of covalent systems. Phys. Rev. B, 1988, 37, 6991‐7000. [41] Bazant, M. Z.; Kaxiras, E.; Justo, J. F. Environment‐ dependent interatomic potential for bulk silicon. Phys. Rev. B, 1997, 56, 8542‐8552. [42] Plimpton, S. Fast parallel algorithms for short‐range molecular dynamics. J. Comput. Phys., 1995, 117, 1‐19. [43] Filippova, V. P.; Kunavin, S. A.; Pugachev, M. S. Calculation of the parameters of the Lennard‐Jones potential for pairs of identical atoms based on the properties of solid substances. Inorg. Mater. Appl. Res., 2015, 6, 1‐4. [44] Finnis, M. W.; Sinclair, J. E. A simple empirical N‐body potential for transition metals. Phil. Mag. A, 1984, 50, 45‐55. [45] Finnis, M. W.; Sinclair, J. E. Erratum. Phil. Mag. A, 1986, 53, 161‐161. [46] Mendelev, M. I.; Han, S.; Srolovitz, D. J.; Ackland, G. J.; Sun, D. Y.; Asta, M. Development of new interatomic potentials appropriate for crystalline and liquid iron. Philos. Mag., 2003, 83, 3977‐3994. [47] Zhou, X. W.; Johnson, R. A.; Wadley, H. N. G. Misfit‐ energy‐increasing dislocations in vapor‐deposited CoFe/NiFe multilayers. Phys. Rev. B, 2004, 69, 144113. [48] Müller, M.; Erhart, P.; Albe, K. Analytic bond‐order potential for bcc and fcc iron—comparison with established embedded‐atom method potentials. J. Phys. Condens. Matter, 2007, 19, 326220.

[49] Pasianot, R.; Farkas, D.; Savino, E. J.; Empirical many‐ body interatomic potential for bcc transition metals. Phys. Rev. B, 1991, 43, 6952‐6961. [50] Ackland, G. J.; Bacon, D. J.; Calder, A. F.; Harry, T. Computer simulation of point defect properties in dilute Fe—Cu alloy using a many‐body interatomic potential. Phil. Mag. A, 1997, 75, 713‐732. [51] Girifalco, L. A.; Weizer, V. G. Application of the Morse potential function to cubic metals. Phys. Rev., 1959, 114, 687‐690. [52] Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool. Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012. [53] Imamova, S. E.; Atanasov, P. A.; Nedialkov, N. N.; Dausinger, F.; Berger, P. Molecular dynamics simulation using pair and many body interatomic potentials: ultrashort laser ablation of Fe. Nucl. Instr. Meth. Phys. Res., 2005, 227, 490‐498. [54] Willaime, F.; Fu, C. C.; Marinica M. C.; Dalla Torre, J. Stability and mobility of self‐interstitials and small interstitial clusters in α‐iron: ab initio and empirical potential calculations. Nucl. Instrum. Methods Phys. Res., Sect. B, 2005, 228, 92‐99. [55] Malerba, L.; Marinica, M.; Anento, N.; Björkas, C.; Nguyen, H.; Domain, C.; Djurabekova, F.; Olsson, P.; Nordlund, K.; Serra A.; Terentyev, D. Comparison of empirical interatomic potentials for iron applied to radiation damage studies. J. Nucl. Mater., 2010, 406, 19‐ 38. [56] Mishin, Y. Interatomic potentials for metals. In Handbook of Materials Modeling. Springer: Dordrecht, 2005, pp. 459‐478. [57] Gades, H.; Urbassek, H. Pair versus many‐body potentials in atomic emission processes from a Cu surface. Nucl. Instrum. Methods Phys. Res., Sect. B, 1992, 69, 232‐241. [58] Shibuta, Y.; Suzuki, T. Melting and nucleation of iron nanoparticles: A molecular dynamics study. Chem. Phys. Lett., 2007, 445, 265‐270. [59] Sipkens, T. A.; Daun, K. J; Titantah, J. T.; Karttunen, M. Quantifying the thermal accommodation coefficient for iron surfaces using molecular dynamics simulations. In ASME 2015 International Mechanical Engineering Congress and Exposition, Houston, Texas, 2015. [60] Mills, K. C. Recommended values of thermophysical properties for selected commercial alloys. Woodhead Publishing: Cambridge, 2002. [61] Basinski, Z. S.; Sutton, A. L. The lattice expansion of iron. Proc. R. Soc. Lond. A, 1955, 229, 459‐467. [62] Hixson, R. S.; Winkler, M. A.; Hodgdon, M. L. Sound speed and thermophysical properties of liquid iron and nickel. Phys. Rev. B, 1990, 42, 6485‐6491. [63] Zhou, X. W.; Wadley, H. N. G.; Johnson, R. A.; Larson, D. J.; Tabat, N.; Cerezo, A.; Petford‐Long, A. K.; Smith, G. D. W.; Clifton, P. H.; Martens, R. L.; Kelly, T. F. Atomic scale structure of sputtered metal multilayers. Acta Mater., 2001, 49, 4005‐4015.



[64] Foiles, S. M. Calculation of the surface segregation of Ni‐ Cu alloys with the use of the embedded‐atom method. Phys. Rev. B, 1985, 32, 7685‐7693. [65] Justo, J. F.; Bazant, M. Z.; Kaxiras, E.; Bulatov, V. V.; Yip, S. Interatomic potential for silicon defects and disordered phases. Phys. Rev. B, 1998, 58, 2539‐2550. [66] Jelinek, B.; Groh, S.; Horstemeyer, M. F.; Houze, J.; Kim, S. G.; Wagner, G. J.; Moitra A.; Baskes, M. I. Modified embedded atom method potential for Al, Si, Mg, Cu, and Fe alloys. Phys. Rev. B, 2012, 85, 245102. [67] Wolf, D. Effect of interatomic potential on the calculated energy and structure of high‐angle coincident site grain boundaries—II.(100) Twist boundaries in Cu, Ag and Au. Acta Metall., 1984, 32, 735‐748. [68] Adams, J. B.; Foiles, S. M.; Wolfer, W. G. Self‐diffusion and impurity diffusion of fee metals using the five‐ frequency model and the embedded atom method. J. Mater. Res., 1989, 4, 102‐112. [69] Foiles, S. M.; Baskes, M. I.; Daw, M. S. Embedded‐atom‐ method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys. Phys. Rev. B, 1986, 33, 7983‐7991. [70] Mishin, Y.; Mehl, M. J.; Papaconstantopoulos, D. A.; Voter, A. F.; Kress, J. D. Structural stability and lattice defects in copper: Ab initio, tight‐binding, and embedded‐atom calculations. Phys. Rev. B, 2001, 63, 224106. [71] Rumble, J. R., Ed., CRC Handbook of Chemistry and Physics, CRC Press: Boca Raton, FL, 2018. [72] Kurochkin, A. R.; Yagodin, D. A.; Borisenko, A. V.; Okhapkin, A. V. Density of copper‐aluminum alloys at temperatures up to 1400°C determined by the gamma‐ ray technique. High Temp., 2013, 51, 197‐205. [73] Brillo, J.; Egry, I. Density determination of liquid copper, nickel, and their alloys. Int. J. Thermophys., 2003, 24, 1155‐1170. [74] Ryu, S.; Cai, W. Comparison of thermal properties predicted by interatomic potential models. Modell. Simul. Mater. Sci. Eng., 2008, 16, 085005. [75] Chizmeshya, A.; Zaremba, E. Interaction of rare gas atoms with metal surfaces: A pseudopotential approach. Surf. Sci., 1989, 220, 443‐470. [76] Zaremba, E.; Kohn, W. Theory of helium adsorption on simple and noble‐metal surfaces. Phys. Rev. B, 1977, 15, 1769‐1781. [77] Bolding, B. C.; Andersen, H. C. Interatomic potential for silicon clusters, crystals, and surfaces. Phys. Rev. B, 1990, 41, 10568‐10585. [78] Cook, S. J.; Clancy, P. Comparison of semi‐empirical potential functions for silicon and germanium. Phys. Rev. B, 1993, 47, 7686‐7699. [79] Keblinski, P.; Bazant, M. Z.; Dash, R. K.; Treacy, M. M. Thermodynamic behavior of a model covalent material

Page 14 of 32

described by the environment‐dependent interatomic potential. Phys. Rev. B, 2002, 66, 064104. [80] Samela, J.; Nordlund, K.; Keinonen, J.; Popok, V. N. Comparison of silicon potentials for cluster bombardment simulations. Nucl. Instrum. Methods Phys. Res., Sect. B, 2007, 255, 253‐258. [81] Godet, J.; Pizzagalli, L.; Brochard, S.; Beauchamp, P. Comparison between classical potentials and ab initio methods for silicon under large shear. J. Phys. Condens. Matter, 2003, 15, 6943–6953. [82] Zhao, J.; Singh, V.; Grammatikopoulos, P.; Cassidy, C.; Aranishi, K.; Sowwan, M.; Nordlund, K.; Djurabekova, F. Crystallization of silicon nanoclusters with inert gas temperature control. Phys. Rev. B, 2015, 91, 035419. [83] Nurminen, L.; Tavazza, F.; Landau, D. P.; Kuronen A.; Kaski, K. Comparative study of Si (001) surface structure and interatomic potentials in finite‐ temperature simulations. Phys. Rev. B, 2003, 67, 035405. [84] Tersoff, J. Empirical interatomic potential for silicon with improved elastic properties. Phys. Rev. B, 1988, 38, 9902‐9905. [85] Erhart, P.; Albe, K. Analytical potential for atomistic simulations of silicon, carbon, and silicon carbide. Phys. Rev. B, 2005, 71, 035211. [86] Rhim, W. K.; Ohsaka, K. Thermophysical properties measurement of molten silicon by high‐temperature electrostatic levitator: density, volume expansion, specific heat capacity, emissivity, surface tension and viscosity. J. Cryst. Growth, 2000, 208, 313‐321. [87] Ohsaka, K.; Chung, S. K.; Rhim, W. K.; Holzer, J. C. Densities of Si determined by an image digitizing technique in combination with an electrostatic levitator. Appl. Phys. Lett., 1997, 70, 423‐425. [88] Gabathuler, J. P.; Steeb, S. Über die Struktur von Si‐, Ge, Sn‐und Pb‐Schmelzen. Zeitschrift für Naturforschung A, 1979, 34, 1314‐1319. [89] Mills, D.; Kolasinski, K. W. Solidification driven extrusion of spikes during laser melting of silicon pillars. Nanotechnology, 2006, 17, 2741‐2744. [90] Dozhdikov, V. S.; Basharin, A. Y.; Levashov, P. R. Two‐ phase simulation of the crystalline silicon melting line at pressures from–1 to 3 GPa. J. Chem. Phys., 2012, 137, 054502. [91] Hara, S.; Izumi, S.; Kumagai, T.; Sakai, S. Surface energy, stress and structure of well‐relaxed amorphous silicon: A combination approach of ab initio and classical molecular dynamics. Surf. Sci., 2005, 585, 17‐24. [92] Dalton, A. S.; Seebauer, E. G. Structure and mobility on amorphous silicon surfaces. Surf. Sci., 2004, 550, 140‐ 148.


Page 15 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic


The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

0 fs

450 fs

900 fs

1350 fs

ACS Paragon Plus Environment

Page 16 of 32

1800 fs

2250 fs

0.2

Page 17 of 32

Fe-Fe

Morse

0.1

Girifalco and Weizer

0 -0.1 U, U2eff [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


-0.2

EAM

Zhou et al.

-0.3

FS

Mendelev et al.

-0.4 -0.5

EAM Tersoff

-0.6

Müller et al.

-0.7

LJ

-0.8 1.5

Pasianot et al.

FS

Ackland et al.

FS

Filippova et al.

2.0

Finnis and Sinclair

2.5

3.0 3.5 ACS Paragon Plus Environment r [Å]

4.0

4.5

5.0

LJ 6-12

Filippova et al.

Zhou et al.

FS

Mendelev et al.

Page 18 of 32

Tersoff

Müller et al.

Ts = 500 K

Solid, α-Fe, BCC

Increasing temperature

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


EAM

Ts = 1500 K

Solid, γ-Fe, FCC

Ts = 1750 K

Solid, δ-Fe, BCC

Ts = 2500 K Liquid


Page 19 of 32

MD/LJ

Fe

Filippova et al.

8000

MD/Tersoff Müller et al.

7500

Basinski et al. Mills

7000

MD/FS

Mendelev et al.

MD/EAM

6500

Zhou et al.

6000 5500 0

BCC (α)

500

FCC (γ)

1000

BCC (δ)

ρm [kg/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

8500


1500 ACS Paragon Plus Environment Ts [K]

Hixson et al. Liquid

2000

2500

3000


0.4 0.3

α

Normal Fe-Ne

γ

Liquid

δ

0.4

FS

Mendelev et al.

0.3

0.2

0.2

0.1

0.1

αt

αn

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

0

Filippova et al.

EAM

Tersoff

Zhou et al.

-0.2 500

1000

Müller et al.

1500 2000 Ts [K]

2500

α

γ

Liquid

δ

FS

Tangential

Mendelev et al.

Fe-Ne

EAM

Zhou et al.

0

LJ

-0.1

Page 20 of 32

3000

Tersoff

-0.1

-0.2 500

(a)

Müller et al.

Filippova et al.

1000

1500 2000 Ts [K] (b)


LJ

2500

3000

Page 21 of 32


0.30 0.25

α

γ

Fe-Ne

0.20

Liquid

δ

0.35

FS

0.30

Mendelev et al.

EAM

α

0.10

0.15

γ

Liquid

δ

FS

Finnis and Sinclair

LJ

Filippov et al.

FS

Mendelev et al.

EAM

Zhou et al.

0.10

0.05 -0.05 500

Fe-Ar

0.20

0.15

0

α

0.25

Zhou et al.

α

α

1 2 3 4 5 6 7 8 9 10 11 Liquid α γ δ 120.20 FS Fe-He 13 FS Finnis and 14 Mendelev et al. Sinclair 15 0.15 16 17 18 190.10 20 21 220.05 23 EAM 24 Zhou et al. 25 0 LJ 26 Tersoff Filippov et al. 27 Müller et al. 28 -0.05 29 500 1000 1500 2000 2500 3000 30 31 Ts [K] 32 33 (a) 34 35 36 37 38 39 40 41

0.05 0 -0.05 500 1000 1500 2000 2500 3000 Ts [K] Tersoff

Müller et al.

LJ

Filippov et al.

(b)


Tersoff

Müller et al.

1000 1500 2000 2500 3000 Ts [K] (c)


0.1

Cu-Cu

0 FS

-0.1 U, U2eff [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

Page 22 of 32

LAMMPS Mendelev

FS

-0.2

EAM

Universal 4 Adams et al.

-0.3

EAM

Morse

-0.4 -0.5 -0.6 1.5

Ackland et al.

Foiles

Wolf

EAM

Morse

Zhou et al.

Gades and Urbassek

2.0

LJ

Filippova et al.

2.5


4.0

4.5

5.0

Page 23 of 32

9000

CRC Handbook

Cu Kurochkin et al.

8500

Brillo and Ergy

ρm [kg/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


8000

MD/EAM Zhou et al.

7500

MD/LJ

Filippova et al.

7000 MD/EAM

Foiles

6500

6000 0

FCC

500

1000

Liquid


2000

2500

3000


LJ

Zhou et al.

Solid

Ts = 500 K

Filippova et al.

EAM

EAM Foiles

Liquid

Ts = 1500 K


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

Page 24 of 32


Page 25 of 32

8

Cu-Gas

Cu-He

Zaremba and Kohn (ZK)

6

Cu-He

4 U [mev]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


Chizmeshya and Zaremba

Cu-Ne

2


0

-2 -4

Cu-Ar


-6

-8 2.0

2.5

3.0

3.5


5.0

5.5

6.0


Liquid

0.30

Cu-He

EAM

Zhou et al.

0.25

EAM

Foiles

FCC

0.5

EAM

Zhou et al. Foiles

0.15

1000

Ts [K] (a)

1500

2000

0 500

Filippova et al.

1000

Cu-Ar

Ts [K]

1500

(b)


Liquid

EAM

Zhou et al.

EAM

Foiles

0.3

0.2

LJ

0.05

Filippova et al.

FCC

0.4

EAM

0.10 LJ

0.6

Cu-Ne

0.20 α

Liquid

α

FCC

α

1 2 3 4 5 6 7 8 9 10 11 120.15 13 14 15 160.10 17 18 19 200.05 21 22 23 24 0 25 26 27 28 -0.05 29 500 30 31 32 33 34 35 36 37 38 39 40 41

Page 26 of 32

2000

LJ

Filippova et al.

0.1

0 500

1000

Ts [K] (c)

1500

2000

Page 27 of 32

2 Justo et al. Zeff = 12

EDIP

Tersoff


Tersoff 3

-2

Tersoff

SW

Tersoff 2

-3


Tersoff

Erhart and Albe

-4 1.0

1.5

2.0

2.5 r [Å] (a)

3.0

Erhart and Albe

0.8

0 -1

Tersoff

g(θ)/max[g(θ)]

1

1.0

Si-Si

EDIP

U2, Ueff 2 [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


3.5

4.0

EDIP


SW


0.6 0.4

Tersoff

0.2

EDIP

0 0


Tersoff Tersoff 3

Tersoff 2


30

60

90 θ [°] (b)

120

150

180

2600 2500

2400 ρm [kg/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


Si

Page 28 of 32

Rhim et al.

MD/SW


Oshaka et al.

Room temperature density

Gabathuler and Steeb

2300 2200

MD/MEAM Jelinek et al.

MD/Tersoff

2100

Tersoff 3

2000 1900 0

MD/EDIP

Justo et al.

c-Si Liquid

500

1000


2000

2500

3000

Page 29 of 32

SW


Tersoff Tersoff 3

EDIP

Justo et al.

MEAM

Jelinek et al.

Ts = 1500 K Solid, c-Si


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


Ts = 2500 K Liquid



0.7 0.6 0.5

c-Si

Normal Si-Ar

0.4

Liquid

0.7

MEAM

0.6

Jelinek et al.

Justo et al.

0.2 0.1

Liquid

Tangential

MEAM

Si-Ar

0.4

0.3

Jelinek et al.

0.3

EDIP

SW

Justo et al.


0.2 0.1

0 -0.1 -0.2 500

c-Si

0.5

EDIP

αt

αn

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

Page 30 of 32

Tersoff Tersoff 2

1000

1500 2000 Ts [K]

0

SW


2500

3000

-0.1 -0.2 500

(a)

Tersoff Tersoff 2

1000

1500 2000 Ts [K] (b)


2500

3000

Page 31 of 32

0.25

c-Si

Liquid

Si-He

0.20 0.15

0.5

MEAM

0.4

Jelinek et al.

EDIP

Justo et al.

0.3

0.10

α

α

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


0.05

-0.05 500

Tersoff 2


1500 2000 Ts [K]

2500

MEAM

Jelinek et al.

EDIP

Justo et al.

0.2

0

SW

1000

Si-Ar

Liquid

0.1

Tersoff

0

c-Si

3000

-0.1 500

(a)

Tersoff Tersoff 2

1000


1500 2000 Ts [K] (b)


SW

2500

3000

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


Page 32 of 32

Effect of Surface Interatomic Potential on Thermal Accommodation

Recommend Documents