Layered structure of Lennard-Jones particle systems confined in a step-shaped gap

We investigate changes in the layered structure of particles confined between flat and step-shaped substrates. Using the Monte Carlo method, the density profiles of argon atoms interacting through a Lennard-Jones potential near a silicon step are calculated for different separation distances. Two different layered structures parallel to the surface of the substrate are observed far from the edge; the transition between the structure takes place within an interval of approximately 1 nm from the edge of the step. The particle distribution in the transition region reflects the formation of additional layers parallel to the contour of the Lennard-Jones potential generated near the edge. Although spatial changes in the layered structure of the nearest layer to the flat substrate across the step edge are small, they induce a non-uniform force on the substrate. If the substrate is flexible, the generated force acts to bend the substrate near the edge. The dependence of the layered structure on the temperature and the density is also evaluated. otherwise a


I. INTRODUCTION
In unconfined liquid or gas phases, particles such as atoms and molecules are randomly distributed. In contrast, when the particles are confined between and interact with parallel substrates separated by a small gap, the isotropy of the particle distribution breaks down and new structures appear. 1 In particular, layered structures are ubiquitous in confined systems, and significantly affect the physical and chemical properties of the corresponding systems. 2,3 For example, an anomalously low dielectric constant of water was recently reported for a confined water system. 4 If the Lennard-Jones interaction between a particle and a substrate is much stronger than that between the particles, the particles are accumulated in two layers, which are located at potential minima between the parallel substrates, except for very small separation gaps. This two-layered structure may become unstable upon increasing the strength of the interaction between the particles, because the Lennard-Jones potential energy rapidly increases as the particles approach. Thus, the total energy can be reduced by forming new layers between the existing ones and decreasing the particle density in each layers. Accordingly, the number of layers usually increases with the separation distance between the substrates. The interlayer distance and the density in each layer have been extensively investigated using Monte Carlo 5 and density functional theory 6-8 simulation methods for parallel configurations.
Because there are many steps on real surfaces, in this study we consider layered structures confined between an infinite flat and a step-shaped substrates. The separation distance between the two substrates shows an abrupt change at the step, and the layered structure far from the step edge is parallel to the surface of the substrate. However, the number of layers and the interlayer distance between layers at the two ends may be different. Therefore, the crossover between the two layered structures must occur near the step edge, we focus on the distribution of particles in this region.
Layered structures are often observed in a thin liquid film between two materials. For example, thin water layers with a total thickness of several nanometers exist between a lipid bilayer and substrate. 9 Recently, supported lipid bilayers (SLBs) have attracted increasing attention in biotechnology. 10,11 A SLB is fabricated by forming a lipid bilayer over a silicon well filled with an aqueous solution. In this system, the distance between the lipid bilayer and the silicon substrate changes abruptly at the edge of the well. Accordingly, the SLB is a good practical example of the system considered system in this study.
Although a large number of studies have focused the layered structure of a liquid between parallel plates, limited information is available on the structure of liquids sandwiched between irregular surfaces. Thus, we begin by considering the distribution of particles interacting with a simple and well-known potential near the edge of a step. This paper is structured as follows: in section II, we analyze the density profiles of argon atoms 12 confined between parallel silicon plates separated by different gaps. Each argon atom interacts with other argon atom, as well as with silicon surfaces through the Lennard-Jones potential. We examine the relationship between the number of layers and the potential energy between parallel substrates showing that the number of layers increases with the separation distance. In section III, the Lennard-Jones potential near the step edge is formulated under the assumption of continuum approximation to investigate the density profiles of argon atoms near step edge. In section IV, The force on the substrate generated by argon atoms interacting through the Lennard-Jones potential is considered. In section V, we examine the dependence of the layered structure near the edge on the temperature and the mean density of argon atoms. Finally, in Sec. VI, we summarize the results of the Monte Carlo simulations and discuss several issues that should be addressed in future studies.

II. LAYERED STRUCTURE OF PARTICLES BETWEEN PARALLEL PLATES
We consider the distribution of particles confined in the space shown in Fig. 1. The lower substrate can be created by forming a groove on its surface. The shaded area indicates the space where the particles are sandwiched, between two plates with a flat and a rectangle wave-shaped surface, respectively. Both the upper and lower substrates are fixed. The minimum and maximum distances between the plates are d and d + H, respectively, where, d is the width of a smaller gap and H is the depth of the groove.
The origin of the chosen coordinate system is located at a d/2 distance from the bottom surface, halfway between the two edges indicated by C and C ′ : l denotes the distance between C and C ′ , i.e., the width of the groove. We assume that the substrates contain Ns atoms per unit volume.
We first calculate the distribution of particles between infinite parallel plates, which are obtained by setting l = 0 as a special case in the considered system. We assume that two particles separated by a distance r interact through a Lennard-Jones potential, expressed as where σ and are material-dependent parameters. For example, the interaction c atoms are σ ArAr = 3.428 Å and ArAr = 0.043 eV. 12 The interaction potential between a confined particle and a substrate is normally calculated by adding the potentials corresponding to all possible particle-substrate atom pairs. However, in this study, the continuum approximation is employed to reduce the computational time in the simulations. The interaction potential energy Vp(h) between an infinite substrate containing Ns atoms per unit volume and a particle at (0, 0, h), where h( > 0) is the height from the surface of the substrate, is given by where r = √ x 2 + y 2 + (h − z) 2 is the distance between the particle and a volume element at (x, y, z). By performing the integration, the potential energy can be expressed as The interaction potential energy between a particle and a substrate is determined only by the height h due to the translation invariant along x − y plane. Rappé et al. 12 presented a table of Lennard-Jones parameters for various elements. The Lennard-Jones parameters between different elements A and B can be approximately determined by the following relationship where σX−X and X −X are the Lennard-Jones parameters between a pair of elements X. Figure 2 shows the potential energy of an argon

ARTICLE
scitation.org/journal/adv atom above a silicon plates. The parameters of the Lennard-Jones potential between an argon and a silicon atom are σ ArSi = 3.636 Å and ArSi = 0.0118 eV. The lattice constant a of the silicon crystal is 5.43 Å, 13 while the number density Ns is 0.05 Å −3 . The interaction potential between an argon atom and a silicon plate reaches a minimum at r = (2/5) 1/6 σ (= 3.12 Å).
Using Monte Carlo simulations, we calculate the distribution of argon atoms sandwiched between parallel silicon plates, separation distance dp for a fixed number of particles (N), the volume (Vs), and temperature (T). First 679 argon atoms are randomly distributed between square silicon plates of 5 nm side and the Metropolis method is applied. The initial number density of argon particles is 0.020 Å −3 , which is approximately equal to the density of argon at the normal boiling point at (87. 3 K), 0.021 Å −3 in bulk. An argon atom interacts not only with the silicon plates but also with other argon particles. However, we initially neglect the interaction between argon particles, and consider the dependence of the distribution on the number of iterations in the Metropolis procedure. This is because the equilibrium distribution can be accurately calculated by ignoring the interaction between the argon particles without MD simulation. By comparing the accurate distribution obtained by numerical integration and that by MD simulation, the accuracy of the MD results, which depends mainly on the iterations in the Metropolis procedure, is evaluated. The interaction between argon atoms is correctly taken into account in the subsequent simulations. In the absence of interactions between particles, the potential energy of a particle at position z, Vpp(z), is given by Figure 3 shows the z position dependence of the potential energy of an argon particle between parallel silicon plates separated by a distance dp of 13.6 Å ( ≈ 2.5 a). Two minima exist at zp = 3.66 Å and zm = − zp. The potential depth Vpp(zm) − Vpp(0) is 0.047 eV, which corresponds to the thermal energy at 545 K. Therefore, a particle trapped near the local minimum can hardly escape to another local minimum at the normal boiling point. The potential depth increases as the separation distance is increased from the 13.6 Å value, and approaches 0.062 eV in dp → ∞ limit. Figure 4 shows the probability density function along the z-axis of argon atoms confined between silicon plates separated by distance dp = 13.6 Å at 90 K, obtained for different number of iterations (Nr = 10, 10 2 , and 10 3 ). In our simulations, when Nr = α, the total number of randomly selected particles is N × α, and the Metropolis algorithm is applied for each selected particle. We employ hard wall boundary conditions along the sides of the space sandwiched by silicon plates; the potential energy is infinite for |x| > 5 nm or |y| > 5 nm. The solid line in Fig. 4 shows the canonical ensemble distribution: where k B is the Boltzmann constant. By increasing the number of iterations, the distribution profiles obtained from the Monte Carlo simulations show a better agreement with the canonical ensemble distribution P(z, 90K). Two peaks are observed at zp and zm in the distribution profile, which indicate the presence of two argon layers parallel to the silicon plates. The density within the layers shows a significant increase compared with its initial value. We then calculate the mean value of the minimum distance between two argon particles, as defined by where rij(t) is the distance between the i-th and j-th particle at the t-th trial in the Monte Carlo simulation. The distance lmin obtained for T 1 = 5 × 10 3 and T 2 =10 4 is 1.48 Å, which is much lower than the equilibrium distance between two argon particles (R ArAr ≡ σ ArAr 2 1/6 = 3.848 Å). Thus, if the interaction between argon atoms is taken into ARTICLE scitation.org/journal/adv account, a strong repulsive force acts between many nearest pairs of atoms.
We now consider the effect of the interaction between argon particles on the distribution profile. Figure 5 shows the density profile along the z-axis for three separation distances (d = 13.6, 16.3 and 21.7 Å, approximately corresponding to 2.5a, 3a, and 4a, respectively). The corresponding numbers of argon particles are N = 679, 815, and 1086, respectively. The comparison of Fig. 4 and Fig. 5(a) reveals the appearance of an additional peak at the origin ascribed to the total potential decrease when some particles near r = rm and rp move to the origin. As mentioned before, the mean minimum distance l min is shorter than the equilibrium distance R: hence, when the interaction between argon particles is taken into account, the repulsive interactions lead to a considerable increases in the total energy. Moreover, the distance between the peak in Fig. 4 and the origin is 3.66 Å, which is almost equal to R ArAr . Thus, the interaction between the particles in the layer at z = zm and those in the newly created layer at the origin is weak. As the separation distance increases, the number of layers also increases, as shown in Figs. 5(b) and (c). The distance between adjacent peaks is 3.3 Å for dp = 16.3 Å, and 3.1 Å for dp = 21.7 Å. Figure 6 shows the number density along the x-axis. Due to the hard wall boundary conditions employed, a vertical layered structure can be observed near the sides of the plot. This layered structure becomes negligible when periodic conditions are employed. Despite slight fluctuations near the origin, the number density is approximately 0.02 (Å −3 ). This implies that if the system size is larger than 5 nm, the density profile near the origin for the distance gaps considered here is approximately equal to that between infinite silicon plates.

III. LAYERED STRUCTURE OF PARTICLES NEAR A STEP EDGE
Here we examine the layered structure of the particles near a step edge. The Lennard-Jones potential between a particle at the (x, 0, z) position and a semi-infinite substrate occupying the space Ω = {(xs, ys, zs)|0 ≤ xs ≤ ∞, −∞ ≤ ys ≤ ∞, −∞ ≤ zs ≤ 0} is given by Vs(x, z) = Ns ∫ Ω d⃗ vV(r(x, z, xs, ys, zs)), where d⃗ v is a volume element and r(x, y, xs, ys, zs) is the distance between (x, 0, z) and (xs, ys, zs). The Lennard-Jones potential consists of the attractive and repulsive − 4 (σ/r) 6 and 4 (σ/r) 12 terms, respectively. The integral of the attractive term, Ia(x, z) can be expressed analytically as where r = √ x 2 + z 2 is the minimum distance between the particle and the edge. Using circular coordinates (x, y) = (r cos θ, r sin θ), Ia can be expressed as The magnitude of Ia decreases proportionally to r −3 , reaching a minimum at θ = 3/4π for a constant r. Similarly, the repulsive term can be expressed in analytical form; however, the calculation of this term near θ = 0 is difficult, due to loss of significant digits. Thus, we employ the Padé approximation near θ = 0. The magnitude of the repulsive term decreases proportionally to r −9 and reaches a minimum at θ = 3/4π. Figure 7 shows a contour plot of the Lennard-Jones potential energy near the edge of the semi-infinite substrate. Darker regions denote smaller potential energies. We then calculate the distribution of argon atoms confined in a step-shaped gap at 90 K for separation distances d = 13.6, 16.3 and 21.7 Å. The step height H and the distance between steps are 3.84 and 10 nm, respectively (see Fig. 1). The system sizes along the

ARTICLE scitation.org/journal/adv
xand y-axes are 30 and 1 nm, respectively, for any separation distance. The number of argon particles is also fixed to 1328 for any d value. By counting the number of particles between z = − d/2 and d/2 at different iterations in the Metropolis procedure, we find that the number density in the narrow gap is approximately 0.02 Å −3 for any value of d. Figure 8 shows the two-dimensional profile of the argon density near the edge of the step in the x−z plane. The density is calculated by averaging the number of particles in a volume element with dx = 0.1 Å and dz = d/500. The ensemble averaging is carried out for the configurations obtained between Np = 5 × 10 4 and 10 5 , starting from 100 different initial configurations. Layered structures parallel to substrates are observed between x = 50 Å and 100 Å, and the numbers of layers are the same as those in Fig. 5. Parallel layered structures are also observed near x = 30 Å, but they are different from those in Fig. 5. Thus, it may be useful to divide the system into three domains, according to the x value: (L) x < 40 Å, (C) 40 Å ≤ x ≤ 55 Å, and (R) 55 Å < x. In general, the layered structure within the (L) domain is almost identical to that near a single infinite substrate, while the layered structure in the (R) domain is almost the same as that observed between infinite parallel substrates separated by a gap d. 14 The most interesting aspect is the transition of the layered structure in domain (C). The layered structure in this region can be considered as a combination of horizontal layers and layers parallel to the contour of the Lennard-Jones potential shown in Fig. 7. To consider the impact of the interactions between argon atoms on the layered structure, we simulate without the interactions between argon atoms. Figure 9 shows two-dimensional profile of the argon density at d = 13.6 Å. The layered structure clearly disappears. Figure 10 shows the mean number densities, averaged from z = − d/2 to z = d/2 as a function of x, for different separation distances. Except near the edge, the density values are around 0.02 Å −3 , which is almost equal to the value considered in the previous section.  The sharp increase near 50 Å for d = 13.6 Å is due to the vertical layered structure formed as a result of the interaction with the vertical side of the lower substrate.
Vertical density profiles obtained for different horizontal positions are shown in Fig. 11. Transitions between the layered structures occur for x between 45 and 55 Å, as shown in 11 Fig. (c), (h), and (m). The top layer interacts strongly with the upper substrate; hence, the value of the corresponding peak is almost independent on the horizontal position. This implies that the influence of the step on the upper layer is small. On the other hand, the peak values corresponding to the lower layers in domain (R) are significantly smaller than those in domain (L). This change occurs within an interval of approximately 1 nm along the x-axis. In summary, the layered structures away from the step can be described in terms of those between parallel plates. Specific changes in layered structure are observed in the region facing the vertical surface of the substrate, where the presence of the edge affects the structure over a range of 1 nm, depending on the parameter σ in the Lennard-Jones potential.

IV. SPATIAL CHANGES IN THE FORCE ACTING ON THE SUBSTRATE
The mean force exerted on a substrate atom by the particles located between parallel substrates does not depend on the horizontal position of the atom. In contrast, if a step exists on the lower substrate, the force acting on the upper substrate can be altered by changes in the layered structure. Here, we define the normalized force acting on the atom positioned at (x, 0, z) in the upper substrate, where xe = 50 Å is the x coordinate of the edge, and Here, ri is the distance between the i-th particle and (x, 0, z), while the bracket denotes an ensemble average. Figure 12 shows plots of . This change is primarily due to the slight difference between the upper substrate-top layer mean distances in domains (L) and (R). As can be seen in Fig. 8(a), the central layer in domain (R), which is located near z = 0, appears brighter and thick than the second layer from the top in domain (L). Thus, the top layer in domain (R) is more strongly attracted to the central layer. As a result, the mean distance between the top layer in domain (R) and the upper substrate increases, resulting in a weaker force.

V. DEPENDENCE OF LAYERED STRUCTURE ON TEMPERATURE AND NUMBER OF PARTICLES
At high temperatures, the layered structure of the Ar particles becomes less pronounced. This is because, to increase entropy, the space occupied by argon particles increase with the temperature. To examine this effect, we calculate the density profile at 100 K, which, however, shows only small changes. The difference between the densities at 90 and 100 K for d = 13.6 Å is shown in Fig. 13(a) (dashed line). The negative difference in density near z = 3.6 Å, where the top layer is located, indicates that the peak value of the density at 100 K is smaller than that at 90 K. On the other hand, the density increases in the 1 -3 Å interval, corresponding to the space between the top and the central layer.

AIP Advances
Conversely, the peak value of the density at 80 K increases. These changes are also observed for other separation distances. In short, increasing the temperature leads to a broadening of each layer without significantly changing the number of particles contained in the layer. Fig. 14 shows the density variation caused by a 10% increase or decrease in the number of particles at 90 K. The peak value of the density increases upon increasing the number of particles; however, the overall shape of the peaks is hardly affected.

VI. CONCLUSIONS
We investigated change in the particle distribution between two substrates separated by a narrow gap. In the first system considered, both substrates are flat, while in the second system a groove was created on one substrate. We conclude that the observed differences in the particle distributions of the two systems are mainly due to two geometrical differences: 1) the abrupt change in the separation between the substrates in the second system; 2) the interactions between the particles and the vertical surface, which are present only in the second system. The first effect leads to a concentration of the particles between substrates with a smaller gap and a vertical shift of the layered structure. The second effect results in the bending of the layered structure near the edge. If the different geometry of the gap is introduced through a gentle ramp function instead a step function, the second effect may be smaller.
The second effect produces structures similar to an interference pattern, as shown in Fig. 8. The flat substrate forces particles to align parallel to its surface. On the other hand, the edge of the substrate ARTICLE scitation.org/journal/adv forces particles to align to the contour of the Lennard-Jones potential generated near the edge. The interference pattern is generated by the combination of these effects. A higher density of particles is produced in the region where the layers generated by the two effects overlap.
In the bulk, the phase of a Lennard-Jones system is determined by a few parameters such as temperature and density. 15 In contrast, in confined systems, geometrical aspects are more important than the above parameters. 14 Increasing the temperature leads to an increased width of the layers. However, the macroscopic structure hardly changes within the range explored in the simulations. Although the microscopic structure depends on the substrate material, the macroscopic structure is almost identical to the layered structure obtained in this study if the interaction is expressed with the Lennard-Jones potential and the ratio d/σ is equal.
We have shown that the force acting on the upper substrate changes across the edge. In the simulations, the substrates are fixed. However, a flexible substrate such as a lipid bilayer will deform in the presence of a uniform force. When the substrate is deformed, the layered structure also changes: thus, the substrate deformation and the layered structure must be determined at the same time. In addition, the thermal fluctuations of the substrate, which results in a less defined layered structure, must also be considered.
The analysis of confined water near the step edges is an important subject of future studies. 16,17 To simulate layered water, we need to incorporate the effects of electric forces and hydrogen bonding. It is likely that the presence of an electric double layer near the edge may significantly affect the layered structures. However, the layered structures obtained in this study may be used as initial configurations of Monte Carlo simulations of confined water near an edge, to reduce the computational time with respect to that corresponding to random initial configurations.