Subatomic-Level Solid / Fluid Boundary of Lennard-Jones Atoms: A Molecular Dynamics Study of Metal-Inert Fluid Interface

: At the molecular scale, the deﬁnition of solid / ﬂuid boundary is ambiguous since its deﬁning precision is comparable to the size of the electron orbitals. It is important to ﬁgure out the sub-atomic-level solid / ﬂuid boundary as the deﬁnition of the solid / ﬂuid interface is related to estimating various properties such as slip length, Kapitza resistance, conﬁned volume, thermodynamic properties, and material properties. In this work, molecular dynamics (MD) simulations were conducted to show the e ﬀ ects of the solid / ﬂuid boundary on estimating thermodynamic properties. Our results reveal that the di ﬀ erent deﬁnitions of solid / ﬂuid boundary can cause a considerable impact on quantitative analysis and even qualitative analysis of a nanoscale system. The solid / ﬂuid boundary for Lennard-Jones atoms is determined within sub-atomic precision via heat transfer MD simulations and microscopic heat ﬂux relation. The result shows that solid / ﬂuid boundary is slightly shifted to the ﬂuid regime as the temperature increase. We suggested a mathematical expression of solid / ﬂuid boundary of LJ atom that is theoretically estimated by ignoring the thermal vibration. The results presented in this work are expected to improve the accuracy of analyzing nanoscale phenomena as well as the continuum-based models for nanoscale heat and mass transport.


Introduction
The recent advances of nanotechnology have motivated the need for understanding molecular-level physics of nano-devices such as energy storage [1], water purification [2,3], electric power generators [4,5], biochips [6], and integrated fuel cells. For those nano-devices, the scale of modeling is comparable to the size of electron orbitals, and therefore, the solid/fluid boundary is ambiguous (see Figure 1). Currently, the solid/fluid boundaries have been defined within one atomic diameter discrepancy. Depending on where the boundary defined, the size of the channel or nanoparticle [7], slip length [8], Kapitza length [9], dynamic properties [10], and even material properties can be differently estimated. Furthermore, defining solid/fluid boundary is related to the location of the boundary condition (i.e., hydrodynamic slip and interfacial thermal resistance) of continuum modeling [10][11][12][13]. These errors may have a minor influence on large systems (~larger than 100 nm) but have a significant effect on the sub-10 nm system. As semiconductor production processes recently reduced to the single digit nanometer sizes [14] and novel nanochannel fabrication techniques have been developed [15,16], a proper definition of the solid/fluid boundary becomes increasingly important. However, there has been little attention on the atomic-level definition of solid/fluid boundary. fabrication techniques have been developed [15,16], a proper definition of the solid/fluid boundary becomes increasingly important. However, there has been little attention on the atomic-level definition of solid/fluid boundary. For small systems, interfacial phenomena play a crucial role due to the large surface-to-volume ratio. However, it is difficult to access the molecular-level details of the interfacial phenomena via experiments. Because of the difficulty, molecular dynamics (MD) simulations have become an important tool to investigate the interfacial phenomena. The MD simulation provides all atomic trajectories, which can be converted into measurable macroscopic properties. Also, the thermodynamic and transport properties of both solids and liquids are well reproduced by MD simulations [17][18][19][20][21][22]. As MD simulation allows us to investigate molecular-level details of the interfacial phenomena, it is necessary to define the solid/fluid boundary with atomic-level accuracy. The atomic-level solid/fluid boundary involves the thermal motion, the discrete nature, the interatomic force penetration. Hence, defining a boundary with atomic precision is complicated a problem.
The followings are the definitions of solid/fluid boundary have been used during the past decades: i) innermost solid layer (center to center) [23][24][25][26][27], ii) midpoint between the innermost solid layer and the first adsorbed layer [9], iii) first zero fluid density (accessible area) [28,29], and iv) the first peak of the adsorption layer [12,30,31]. The problems arising from the absence of atomic-level boundary definition have been reported in several documents. Nagayama et al. [8] point out that the determination of velocity slip is sensitive to the defined boundary. Similarly, Han et al. [9] show the interfacial thermal resistance considerably depends on the solid/fluid boundary definition. Kim et al. [10] report that the evaluations of the local properties adjacent to the solid/fluid interface are affected by the boundary definition. Ramos-Alvarado et al. [19] documented that a proper definition of the liquid/gas interface is necessary for accurate measurement of the contact angle. Despite these reported problems, there is no classical MD study on elucidating the atomic-level definition of solid/fluid boundary.
The rest of the paper is organized as follows: we first present the details of the MD simulation. we then show the effects of the defined boundary on thermodynamic properties and material. The atomic-level boundaries for Ag-Ar, Ag-He, and Cu-Ar are estimated by heat transfer MD simulation and microscopic heat transport equation. We further study the temperature effects on the boundary definition. At the end of the article, a summary of the findings is presented.

Simulation Detail
To achieve the objectives of the present work, MD simulations were performed. For thermodynamic properties, equilibrium MD simulation is conducted for various diameters under the canonical ensemble (NVT). Non-equilibrium MD simulation is conducted to investigate a proper location of solid/fluid boundary. The metallic materials (Ag, Cu) are considered as a channel. For the fluid molecules, inert species (Ar, He) were considered. By testing the two solid models and two For small systems, interfacial phenomena play a crucial role due to the large surface-to-volume ratio. However, it is difficult to access the molecular-level details of the interfacial phenomena via experiments. Because of the difficulty, molecular dynamics (MD) simulations have become an important tool to investigate the interfacial phenomena. The MD simulation provides all atomic trajectories, which can be converted into measurable macroscopic properties. Also, the thermodynamic and transport properties of both solids and liquids are well reproduced by MD simulations [17][18][19][20][21][22]. As MD simulation allows us to investigate molecular-level details of the interfacial phenomena, it is necessary to define the solid/fluid boundary with atomic-level accuracy. The atomic-level solid/fluid boundary involves the thermal motion, the discrete nature, the interatomic force penetration. Hence, defining a boundary with atomic precision is complicated a problem.
The followings are the definitions of solid/fluid boundary have been used during the past decades: (i) innermost solid layer (center to center) [23][24][25][26][27], (ii) midpoint between the innermost solid layer and the first adsorbed layer [9], (iii) first zero fluid density (accessible area) [28,29], and (iv) the first peak of the adsorption layer [12,30,31]. The problems arising from the absence of atomic-level boundary definition have been reported in several documents. Nagayama et al. [8] point out that the determination of velocity slip is sensitive to the defined boundary. Similarly, Han et al. [9] show the interfacial thermal resistance considerably depends on the solid/fluid boundary definition. Kim et al. [10] report that the evaluations of the local properties adjacent to the solid/fluid interface are affected by the boundary definition. Ramos-Alvarado et al. [19] documented that a proper definition of the liquid/gas interface is necessary for accurate measurement of the contact angle. Despite these reported problems, there is no classical MD study on elucidating the atomic-level definition of solid/fluid boundary.
The rest of the paper is organized as follows: we first present the details of the MD simulation. we then show the effects of the defined boundary on thermodynamic properties and material. The atomic-level boundaries for Ag-Ar, Ag-He, and Cu-Ar are estimated by heat transfer MD simulation and microscopic heat transport equation. We further study the temperature effects on the boundary definition. At the end of the article, a summary of the findings is presented.

Simulation Detail
To achieve the objectives of the present work, MD simulations were performed. For thermodynamic properties, equilibrium MD simulation is conducted for various diameters under the canonical ensemble (NVT). Non-equilibrium MD simulation is conducted to investigate a proper location of solid/fluid boundary. The metallic materials (Ag, Cu) are considered as a channel. For the fluid molecules, inert species (Ar, He) were considered. By testing the two solid models and two fluid models, the effect of the potential parameters was investigated. We assumed that there is no chemical reaction as we utilized classical MD simulation. The interatomic potentials are modeled by 12-6 Lennard-Jones (LJ) potential. The potential parameter used in this work is summarized in Table 1. Lorentz-Berthelot combining rules are utilized for the LJ parameters between different atoms. All MD simulations were conducted using LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [32] with GPU (Graphics Processing Unit) package [33], while OVITO (Open Visualization Tool) [34] was used for the visualization.  Figure 2 shows the detailed simulation setup for the MD simulation with heat transfer used in this study. This heat transfer MD setup has extensively studied for interfacial thermal resistance [37][38][39]. In this work, three different channel widths were considered. The periodic boundary condition is applied to the x and y directions, and the fixed boundary condition is utilized to the z-direction. The number of fluid molecules is determined so that it can reproduce a proper bulk density at the center of the fluid block (bulk region). An ideal flat solid surface is considered. The outmost second to fourth layers were locally thermostated. The thermostat generates and removes the same amount of heat by manipulating the velocity magnitude conserving the net momentum. Hence, no heat accumulation is allowed during the simulation and a constant temperature gradient can be achieved. Therefore, 0.01 eV/ps z-directional heat flux is maintained during the steady state. To reaches the desired temperature, the system is first integrated under the canonical ensemble (NVT) for 0.5 ns1 .0 ns without the local thermostats. Then, micro-canonical ensemble (NVE) is used and the two local thermostats are activated. The heat flux data is collected during long time spans (200 ns~400 ns) to ensure reliable statistics. For the equilibrium MD simulation, the simulated system is the same as the heat transfer MD, but it is only integrated under NVT ensemble at 100 K without the local thermostats during 10 ns.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 11 fluid models, the effect of the potential parameters was investigated. We assumed that there is no chemical reaction as we utilized classical MD simulation. The interatomic potentials are modeled by 12-6 Lennard-Jones (LJ) potential. The potential parameter used in this work is summarized in Table  1. Lorentz-Berthelot combining rules are utilized for the LJ parameters between different atoms. All MD simulations were conducted using LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [32] with GPU (Graphics Processing Unit) package [33], while OVITO (Open Visualization Tool) [34] was used for the visualization.  is center-to-center distance. The solid blocks have a width of seven lattice constants except for the thermostated layers and nine lattice constants along x and y directions. Figure 2 shows the detailed simulation setup for the MD simulation with heat transfer used in this study. This heat transfer MD setup has extensively studied for interfacial thermal resistance [37][38][39]. In this work, three different channel widths were considered. The periodic boundary condition is applied to the x and y directions, and the fixed boundary condition is utilized to the z-direction. The number of fluid molecules is determined so that it can reproduce a proper bulk density at the center of the fluid block (bulk region). An ideal flat solid surface is considered. The outmost second to fourth layers were locally thermostated. The thermostat generates and removes the same amount of heat by manipulating the velocity magnitude conserving the net momentum. Hence, no heat  Figure 3 shows a typical density and pressure distribution near the metallic surface [23]. The fluid molecules form a dynamic structure near the surface called density layering. This non-homogeneous distribution of fluid atoms results in a high local pressure variation near the surface [40]. Several terminologies are defined in order to closely investigate the solid/fluid interface. The center-to-center height (H cc ) is defined as the distance between the innermost solid layers. The innermost solid layer is defined as the most probable location of innermost solid atoms. The parameter δ b is the distance between the innermost solid layer and the defined solid/fluid boundary as shown in Figure 3. As the definition of solid/fluid boundary is assumed as a plane in this work, a new definition of solid/fluid boundary can be easily applied to the continuum analysis [12]. Then, the effective channel height is a function of

Defining Terminologies in a Confined System
Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 11 Figure 3 shows a typical density and pressure distribution near the metallic surface [23]. The fluid molecules form a dynamic structure near the surface called density layering. This non-homogeneous distribution of fluid atoms results in a high local pressure variation near the surface [40]. Several terminologies are defined in order to closely investigate the solid/fluid interface. The center-to-center height ( ) is defined as the distance between the innermost solid layers. The innermost solid layer is defined as the most probable location of innermost solid atoms. The parameter is the distance between the innermost solid layer and the defined solid/fluid boundary as shown in Figure 3. As the definition of solid/fluid boundary is assumed as a plane in this work, a new definition of solid/fluid boundary can be easily applied to the continuum analysis [12]. Then, the effective channel height is a function of ; ≡ − 2 .

Effects of the Defined Boundary on the Average Density and Pressure
The thermodynamic properties of the confined system are calculated with various boundary definitions to demonstrate the impact of the boundary definitions on the property estimation. The average fluid density, , and pressure, , are used to interpret the nanochannel transport phenomena [30,[41][42][43]. These properties are coupled with the volume of fluid and the volume can be differently estimated depending on solid/fluid boundary. Assuming the microscopic solid/fluid boundary as a plane that is away from the innermost solid layers, the volume of confined fluid can be expressed as follows: where is the cross-section area of confined fluid (xy plane), is the center to center height and is the boundary location.
Then, the average fluid density between the two boundaries is described as follows: where is the number of fluid between the defined boundary. Similarly, the microscopic definition of the pressure of the confined fluid with boundary definition can be expressed as follows:

Effects of the Defined Boundary on the Average Density and Pressure
The thermodynamic properties of the confined system are calculated with various boundary definitions to demonstrate the impact of the boundary definitions on the property estimation. The average fluid density, ρ avg , and pressure, P avg , are used to interpret the nanochannel transport phenomena [30,[41][42][43]. These properties are coupled with the volume of fluid and the volume can be differently estimated depending on solid/fluid boundary. Assuming the microscopic solid/fluid boundary as a plane that is δ b away from the innermost solid layers, the volume of confined fluid can be expressed as follows: where A c is the cross-section area of confined fluid (xy plane), H cc is the center to center height and δ b is the boundary location. Then, the average fluid density between the two boundaries is described as follows: where n fluid is the number of fluid between the defined boundary. Similarly, the microscopic definition of the pressure of the confined fluid with boundary definition can be expressed as follows: where F ij is an interatomic force vector between i and j atoms and r ij is a relative positional vector from i to j atom.
As can be seen in Figure 4, the average fluid densities and pressure show non-linear variations with respect to the defined boundary. The non-linearity of the average properties are resulted by the dynamical structuring of fluid molecules near the solid surface, which induce the local variation of density and pressure profile. The different crystal orientations show a similar variation of ρ avg and P avg . It is noted that the largest discrepancy of the average fluid density is 16% for the smallest channel studied (i.e., 3.47 nm). Considering the fact that the dense fluid generally has tiny compressibility, this discrepancy can substantially affect the interpretation and understanding of the physical behavior of the confined fluid. Also, the average properties for smaller channels are more sensitive to the definition of solid/fluid boundary indicating that defining solid/fluid boundary becomes more important for a small system. This is because the smaller channel has a larger surface-to-volume ratio and it enhances the effect of the structural inhomogeneity of fluid near the surface. Thus, this size effect can be roughly estimated by non-dimensional number: 2σ sf/H cc , which indicates the ratio of boundary uncertainty over the channel size.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 5 of 11 where is an interatomic force vector between i and j atoms and is a relative positional vector from i to j atom.
As can be seen in Figure 4, the average fluid densities and pressure show non-linear variations with respect to the defined boundary. The non-linearity of the average properties are resulted by the dynamical structuring of fluid molecules near the solid surface, which induce the local variation of density and pressure profile. The different crystal orientations show a similar variation of and . It is noted that the largest discrepancy of the average fluid density is 16% for the smallest channel studied (i.e. 3.47 nm). Considering the fact that the dense fluid generally has tiny compressibility, this discrepancy can substantially affect the interpretation and understanding of the physical behavior of the confined fluid. Also, the average properties for smaller channels are more sensitive to the definition of solid/fluid boundary indicating that defining solid/fluid boundary becomes more important for a small system. This is because the smaller channel has a larger surface-to-volume ratio and it enhances the effect of the structural inhomogeneity of fluid near the surface. Thus, this size effect can be roughly estimated by non-dimensional number: 2σ , which indicates the ratio of boundary uncertainty over the channel size.

Local Thermal Equilibrium Near the Solid/Fluid Boundary
The local thermodynamic equilibrium near the solid/fluid boundary is tested using MD simulation. The slab bins are divided into 0.05 Å intervals, which is small enough to capture the full resolution of the structural inhomogeneity. The local thermodynamic equilibrium must be satisfied to use the mathematical relations of thermodynamic properties driven from statistical mechanics [44]. This can be tested by comparing local velocity distributions to the Maxwell-Boltzmann distribution and quantified by the correlation coefficient. The atomic velocity of the slab bean is collected for 200 ns, which is a relatively long time considering the typical data acquisition time of 2 ns. As shown in Figure 5a, the nucleus of solid and fluid atoms can occupy 0 < z < 0.1σ sf and 0.73σ sf < z as a result of the thermal motions and interactions between solid and fluid molecules. The correlation coefficient, R, represents how well the velocity distribution of MD simulation matches with Maxwell-Boltzmann distribution indicating a perfect match for R = 1 and no correlation for R = 0. Near the edge of the thermal oscillations, the statistical significance is rapidly broken (see Figure 5b). Statistical Depletion Region (SDR) is defined on the basis of 99% correlations as depicted in Figure 5a. The range of SDR can vary with the thermodynamics states and molecular characteristics as it is closely related to the density distribution. It is speculated that the solid/fluid boundary is located at the SDR where no nucleus occupies. SDR is occupied by electrons and an explicit investigation of electrons in SDR must be carried out by quantum mechanical methods. It is noted that the effect of electrons is implicitly considered by LJ potential in this study.

Local Thermal Equilibrium Near the Solid/Fluid Boundary
The local thermodynamic equilibrium near the solid/fluid boundary is tested using MD simulation. The slab bins are divided into 0.05 Å intervals, which is small enough to capture the full resolution of the structural inhomogeneity. The local thermodynamic equilibrium must be satisfied to use the mathematical relations of thermodynamic properties driven from statistical mechanics [44]. This can be tested by comparing local velocity distributions to the Maxwell-Boltzmann distribution and quantified by the correlation coefficient. The atomic velocity of the slab bean is collected for 200 ns, which is a relatively long time considering the typical data acquisition time of 2 ns. As shown in Figure 5a, the nucleus of solid and fluid atoms can occupy 0 < < 0.1 and 0.73 < as a result of the thermal motions and interactions between solid and fluid molecules.
The correlation coefficient, R, represents how well the velocity distribution of MD simulation matches with Maxwell-Boltzmann distribution indicating a perfect match for R = 1 and no correlation for R = 0. Near the edge of the thermal oscillations, the statistical significance is rapidly broken (see Figure 5b). Statistical Depletion Region (SDR) is defined on the basis of 99% correlations as depicted in Figure 5a. The range of SDR can vary with the thermodynamics states and molecular characteristics as it is closely related to the density distribution. It is speculated that the solid/fluid boundary is located at the SDR where no nucleus occupies. SDR is occupied by electrons and an explicit investigation of electrons in SDR must be carried out by quantum mechanical methods. It is noted that the effect of electrons is implicitly considered by LJ potential in this study.

Determination of Solid/Fluid Boundary with Microscopic Heat Flux Relations
To examine the solid/fluid boundary with atomic accuracy, heat transferred MD simulation is conducted. The mathematical description of microscopic heat flux expression is expressed as follows [45]: Where is the velocity of atom i. , is the total energy. is the position vector of atom i and j. is the force vector between atom i and j. Equation (4) includes a volume term, which is determined by the solid/fluid boundary. Inserting Equation (1) to Equation (4), we can derive Equation (5).

Determination of Solid/Fluid Boundary with Microscopic Heat Flux Relations
To examine the solid/fluid boundary with atomic accuracy, heat transferred MD simulation is conducted. The mathematical description of microscopic heat flux expression is expressed as follows [45]: where v i is the velocity of atom i. E i , is the total energy. r ij is the position vector of atom i and j. F ij is the force vector between atom i and j. Equation (4) includes a volume term, which is determined by the solid/fluid boundary. Inserting Equation (1) to Equation (4), we can derive Equation (5).
From the setup of MD simulation, the average heat flux in a steady state is given as 0.01 eV/ps, which is the rate of thermal energy generating and removing by the local thermostats. Thus, the only unknown variable in Equation (5) is the solid/fluid boundary location, δ b , and thus, the solid/fluid boundary can be estimated.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 11 From the setup of MD simulation, the average heat flux in a steady state is given as 0.01 / , which is the rate of thermal energy generating and removing by the local thermostats. Thus, the only unknown variable in Equation (5) is the solid/fluid boundary location, , and thus, the solid/fluid boundary can be estimated.  Figure 6 shows the calculated heat flux with respect to for Ag-Ar, Ag-He, and Cu-Ar system at 300 K. The calculated heat flux increases as the boundary definition shift to the liquid region. The relative difference between the calculated and the given heat flux is between −8% to +12% within shifting. This observation confirms the importance of the definition of the solid/fluid boundary for the quantitative analysis of heat transport. 0.01 eV/ps heat flux is archived at approximately / = 0.445 for Ag-Ar interface regardless of the channel heights. This result agrees with the location that recovers the local viscosity near the wall (i.e. / = 0.44) report by Kim et al. [10]. In the case of Ag-Ar interface at T = 100 K, this location of solid/fluid boundary can considerably reduce the discrepancies of average density from 7.94% to 1.89% and average pressure from 8.81% to 1.75%. It is  Figure 6 shows the calculated heat flux with respect to δ b for Ag-Ar, Ag-He, and Cu-Ar system at 300 K. The calculated heat flux increases as the boundary definition shift to the liquid region. The relative difference between the calculated and the given heat flux is between −8% to +12% within σ sf shifting. This observation confirms the importance of the definition of the solid/fluid boundary for the quantitative analysis of heat transport. 0.01 eV/ps heat flux is archived at approximately δ b /σ sf = 0.445 for Ag-Ar interface regardless of the channel heights. This result agrees with the location that recovers the local viscosity near the wall (i.e., δ b /σ sf = 0.44) report by Kim et al. [10]. In the case of Ag-Ar interface at T = 100 K, this location of solid/fluid boundary can considerably reduce the discrepancies of average density from 7.94% to 1.89% and average pressure from 8.81% to 1.75%.
It is noted that the difference of LJ diameter of Ar and He is approximately 29%, but the estimated boundaries show only 9% difference (See Figure 6a,b) when it is normalized by σ sf . Similarly, the LJ diameter of Cu and Ag show 13% difference but the estimated boundary shows about 2% difference (See Figure 6a,c). This result implies that the LJ diameter of solid-fluid pair is one of the key factors for solid/fluid boundary. This variation of heat flux affects the estimation of thermal conductivity using Fourier's law (so-called direct method [18,39]). Then, not only quantitative analysis but also qualitative analysis can be disturbed by boundary definition. The two boundary definitions δ b = 0 and δ b = σ sf lead to two contradicting conclusions on size effects of the thermal conductivity: as the channel size increase, the thermal conductivity increases (δ b = 0) or decreases (δ b = σ sf ) (See Figure S1).

Effects of Thermal Motion on Solid/Fluid Boundary
To gain a better understanding of solid/fluid boundary estimated by the previous section, the temperature effects of solid/fluid boundary is studied. The solid/fluid interface involves thermal motions, the structure of solid, potential energy between solid and fluid atoms. In this study, the quantum mechanical temperature effect on the potential energy surface does be not considered, which may have an impact on solid/fluid boundary. An identical LJ potential is used in all temperature range tested.
The solid/fluid boundary location is calculated with LJ potential, the lattice constant of FCC (Face Centered Cubic) crystal by ignoring thermal motion. It is known that monoatomic fluid molecules adjacent to the FCC (001) surface forms another layer of the FCC-like structure to minimize the energy [46]. If the thermal energy is very small, the fluid atoms may form a perfect FCC structure near the surface as shown in Figure 7a. The solid/fluid boundary is represented as the red line that has a three-dimensional corrugation pattern (See Figure S2). For simplicity, the solid/fluid boundary is assumed as a plane in this work. Then, the averaged location of solid/fluid boundary can be described as follows: where a is lattice constant and r m is the location of a minimum 12-6 LJ potential between solid atom and fluid atom (i.e., r m = 2 1/6 σ sf ).
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 11 noted that the difference of LJ diameter of Ar and He is approximately 29%, but the estimated boundaries show only 9% difference (See Figure 6a,b) when it is normalized by . Similarly, the LJ diameter of Cu and Ag show 13% difference but the estimated boundary shows about 2% difference (See Figure 6a,c). This result implies that the LJ diameter of solid-fluid pair is one of the key factors for solid/fluid boundary. This variation of heat flux affects the estimation of thermal conductivity using Fourier's law (so-called direct method [18,39]). Then, not only quantitative analysis but also qualitative analysis can be disturbed by boundary definition. The two boundary definitions = 0 and = lead to two contradicting conclusions on size effects of the thermal conductivity: as the channel size increase, the thermal conductivity increases ( = 0) or decreases ( = ) (See Figure  S1).

Effects of Thermal Motion on Solid/Fluid Boundary
To gain a better understanding of solid/fluid boundary estimated by the previous section, the temperature effects of solid/fluid boundary is studied. The solid/fluid interface involves thermal motions, the structure of solid, potential energy between solid and fluid atoms. In this study, the quantum mechanical temperature effect on the potential energy surface does be not considered, which may have an impact on solid/fluid boundary. An identical LJ potential is used in all temperature range tested. The solid/fluid boundary location is calculated with LJ potential, the lattice constant of FCC (Face Centered Cubic) crystal by ignoring thermal motion. It is known that monoatomic fluid molecules adjacent to the FCC (001) surface forms another layer of the FCC-like structure to minimize the energy [46]. If the thermal energy is very small, the fluid atoms may form a perfect As the temperature increases the solid/fluid boundary shifts toward the fluid regime. A similar shifting is observed during the study of the temperature dependence of ionic structure [47]. It is noted that the minimum temperature tested via MD simulation is 35 K because below this temperature, He atoms are solidified at the Ag interface and thus the method used in this paper is no longer works well. At below 100 K, the solid/fluid boundary shows a linear behavior with respect to the temperature and the linear extrapolation of those data show a good agreement with Equation (6) at T~0 K. At a temperature larger than 100 K, the behavior is non-linear, and it seems saturated after a certain temperature. Additional research is needed to elucidate the cause of this saturation.

Conclusions
In this work, the microscopic definition of solid/fluid boundary is investigated using classical MD simulation with metallic solid and monoatomic fluid. It is demonstrated that the average density and pressure of the nanochannel vary depending on where the solid/fluid boundary is defined. The variation in property estimation is more severe in smaller channels. This size effect can be approximately evaluated by the dimensionless number 2σ sf /H cc . The local thermodynamic equilibrium is tested near the solid/fluid boundary. It is found that the local equilibrium collapses near the edge of thermal oscillations of the nucleus, creating a statistically depleted region. This region is not occupied by both solid and fluid nucleus. This implies that the solid/fluid boundary is located within the depletion region. Solid/fluid boundary is determined with a sub-atomic precision through the microscopic heat flux equation and non-equilibrium MD simulation. The results show that the microscopic heat flux equation is met in a nearly identical solid/fluid boundary regardless of the channel size. The determined location of solid/fluid boundary for the different materials shows small deviations when it is normalized by σ sf . Therefore, it is implied that the LJ diameter of solid/fluid, σ sf , plays a key role in the microscopic solid/fluid boundary. The temperature effect on the solid/fluid boundary is further studied for Ar/He interface. It is shown that the solid/fluid boundary of LJ atom can be theoretically estimated by ignoring the thermal vibration. This theoretical estimation shows a good agreement with the extrapolation of the solid/fluid boundary estimated by MD simulation. The results in this study are expected to contribute to understanding nanoscale fluid transport by improving the accuracy and consistency of property calculation. Moreover, it is expected that the definition of solid/fluid boundary improves the continuum-based models for nanoscale heat and mass transport. Future works are required for molecules with complex geometry such as a water molecule. Also, the solid/fluid boundary at the charged surface must be further studied; the solid/fluid boundary could be very different to the non-charged surface as the charged surface creates the Electrical Double Layer (EDL). It is necessary to study how much this new definition of boundary improves the continuum-based model. Also, a quantum mechanical approach may need to elucidate the role of electrons in the solid/fluid boundary.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-3417/9/12/2439/s1, Figure S1: The thermal conductivity calculated by the Fourier's law of conductance with a various solid/fluid boundary for (a) Ag(001)/Ar interface and (b) Ag(111)/Ar interface. The errors are estimated as the summation of the standard error of the heat flux equation and the uncertainty of temperature gradient of the linear least square fitting. The black dash line represents the thermal conductivity calculated by Green-Kubo Method with a homogeneous LJ argon system, which contains 12,195 argon atoms with the approximately same thermodynamic state with the bulk region of nanoconfinement, Figure S2: Density distributions of the Ag(001)/Ar nanoconfinement: (a) two-dimensional contour, (b) one-dimensional profile for convex (α − α ) and concave (β − β ) corrugation. (c) semi-three-dimensional density profile of the adsorbed fluid near the Ag surface.