Structure of free surface crystallization of Hg : a first attempt

Oscillatory density profiles (DP) are a well-known feature occurring at the free surfaces of liquid metals. We analyse first the layered structure of Hg at or below the melting temperature with a simple interaction model [J.-M. Bomont and J.-L. Bretonnet, J. Chem. Phys. 124, 054504 (2006)], by molecular dynamics simulation in slab geometry. Then, by investigating the in-plane structure, some state transitions are observed. In particular, crystalline planes, whose axes are tilted by several degrees from the surface normal, appear and seem to be signalized by flat DP.


Introduction
The structure of the free liquid surface has attracted a substantial amount of experimental and theoretical works [1,2], with emphasis on the possible existence of liquid surface layers.This is a long standing problem whose goal is to understand how surface induced-layering depends on the interactions between atoms [3].It had long been believed that layering was absent at the liquid surface of simple fluids with isotropic interactions, and consequently thought to be unique to metallic liquids, many of them being characterized by a high critical temperature T C and a low melting temperature T M .However, the previous statement has been called into question recently.Chacon et al [4] have pointed out that metallic bonding did not play a crucial role in the formation of a layered surface and shown that layering would be a universal feature of any substance (metallic or not) provided that it has a sufficiently low T M /T C ratio.According to the authors, such a low ratio is the result of a broad and shallow potential interaction, typical of socalled "cold-liquids".In that framework, surface layering would be observed below a threshold that is T/T C ~0.2.Recent ab-initio simulation [5] results for sodium provide an indirect support to the finding of Chacon et al by concluding that neither Friedel oscillations, that are characteristic of liquid metals, nor rearrangements induced by under-coordinated atoms at the surface have any influence on the surface layering.In addition, the emphasis on the existence of the T/T C threshold has been straightforward supported by very recent X-ray reflectivity experiments, on the one hand, showing surface layering in molecular non-liquid-crystalline, nonmetallic liquids with low T M /T C ratios [6], and by computer simulation studies on the other hand, predicting a significant structure, extending several atomic layers into the bulk liquid-Hg (l-Hg) at room and melting temperatures [7].In the latter study, an accurate broad and shallow pair interaction model, presenting no Friedel oscillations has been used [8].
In this paper, we report a study on the liquid-vapour (L/V) interface of pure l-Hg at thermodynamic conditions at and below the melting point.This has been carried out by using molecular dynamics method, where the forces acting on the atoms are computed from a simple pair interaction model [8], resulting from a low freezing point (T M ~ 235K) and a high critical point (T C ~ 1750K).By investigating the behaviour of both DP and in-planestructure in conditions that are not yet accessible to experimental measurements, some state transitions are met.In addition, tilted crystalline planes are shown to be signalized by flat DP along the surface normal.

Pair potential model
In a metal the nature of the interactions changes markedly across the interface: the bonding changes from metallic (with delocalized valence electrons in the bulk liquid) to van der Waals type (with localized electrons on the atoms in the vapor), leading to a metal-nonmetal (M-NM) transition.A distinctive feature of l-Hg is that the M-NM transition does not occur at the liquid-vapor critical density (ρ C = 5.3 g/cm³) but rather gradually, in the range from 11 to 8 g/cm³.This means that the M-NM transition occurs across the interface, making its study a difficult and challenging problem.Indeed, we are aware of only two previous theoretical studies [9,10], based on Monte Carlo simulations involving pseudopotential EPJ Web of Conferences formalism, which were performed long before the significant advances mentioned in the previous section.To date, the x-ray reflectivity data for l-Hg at both melting and room temperatures [11], have only been reproduced recently by theoretical calculation including self-consistent classical simulations [9].For the present investigation, we have used an interaction model, developed in a previous work [8] that reads ( This model has already shown its capability to properly predict the liquid branch of the coexistence L/V phase diagram in terms of temperature and density (implying that it effectively incorporates subtle many-body effects).It presents a broad and shallow shape, and, contrary to pseudopotential theory, does not exhibit long-range Friedel oscillations.Excepting A 1 which satisfies an empirical law for T  T M , the four remaining parameters are constant along the whole liquid-vapour coexistence curve.Additional details are given in Ref. [8].In the present study, forces acting on the atoms have been computed with A 1 (T  T M ) = A 1 (T = T M ), that is to say we did not account for the interactions changes across the interface in this first attempt.

Simulation details
The calculations were performed in slab geometry.In order to start the simulation from a well equilibrated bulk configuration we proceeded in two ways.It is usual to equilibrate the system, whose initial dimensions are L x , L y and L z , where L x = L y , containing a chosen number of atoms so that the initial density of the slab is equal to the known bulk density at the considered temperature.Following this prescription, we equilibrated a liquid of 10 3 atoms at a temperature of 273K and a density of 13.6 g/cm³ (the expected density of l-Hg at coexistence).The simulation cell had the following dimensions: L x = L y = 24.56Å and L z = 40.61Å.However, it is important to emphasize that this usual method somehow looses a predictive character as seen in Ref. [12], in the sense that it can hide some deficiencies of the pair interaction model in predicting (once the interface has been created) the expected bulk density at the considered temperature.According to us, in a liquid-vapour interface study, the bulk density (the density in the inner part of the slab) has to be envisaged much more as a necessary final result rather than an input data of the simulation.That is the reason why we preferred to equilibrate the system with a lower number of ions than the one required by the usual method.Namely, we kept the dimensions of the simulation box unchanged, but containing 864 ions only.This second choice allows us to judge the accuracy of the employed pair interaction model in predicting the expected bulk density in the inner part of the slab after the interface has been created.In both cases, the equations of motion were integrated using a fourth order Gear predictor-corrector algorithm.The timestep was chosen as 1 fs.Temperature control was achieved through the use of a Gaussian thermostat which fixes the total kinetic energy and periodic boundary conditions were applied in the x and y directions of space.In the present calculation, the interfaces are allowed to relax for 2 ns and properties of interest are averaged over another run of 6 ns.

Results and discussion
We investigated the behaviour of the DP for temperatures at and below 235K.The initial equilibrated liquid slab has been cooled down at several temperatures, namely 235, 230, 220, 210 and 205K.As can be seen in figure 1, the specific shape of the DP is qualitatively similar to those calculated from simulations for other potential models [12].The first peak is always lower than the second one.At T = 235K, the density in the central region is equal to 13.7 g/cm 3 , in excellent agreement with the measured bulk density of l-Hg at this temperature.
The oscillations penetrate four atomic diameters into the bulk region.As the temperature is decreased down to 210K, changes are observed in the shape of DP, like a small increase in number and amplitude of the density oscillations.Moreover, the density increases slightly in the bulk region, indicating some structural modifications.In contrast, at T = 205K, a spectacular change in the behaviour of density profile occurs.The slab shrinks abruptly, the number of oscillations decreases substantially and the obtained DP is flat.In addition, it should be noticed that a brutal increase in the density (up to 15.8 g/cm 3 ) accompanies this drastic change.As far as we know, no flat DP has been yet reported in the literature.This feature prompted us to perform a structural analysis across the layered interface.In order to examine the in-plane structure of our model system, we have computed the transverse pair correlation function, g T (r), to study thin sections sliced parallel to the interface.As can be seen from Fig. 2, a gradual change in the surface structure is observed with decreasing temperature.While the transverse pair correlation function is still typical of the liquid state in the two outermost layers (1 and 2) at T = 230K and T = 210K, this is no more the case at T = 205K.In fact, g T (r) seems to indicate a face-centered cubic (FCC) structure in the second layer, while, in layer 1, some vestiges of the crystalline intralayer order of layer 2 remain.It should be 01026-p.2LAM14 noted that such behaviour has previously been reported for Ni [13].In the bulk region, while g T (r) reflects a more pronounced crystalline structure at T = 205K, it exhibits at higher temperatures a shoulder, located at its second peak (indicated by the arrow), that could be interpreted in terms of an amorphous state.This feature could be the result of vestiges of the crystalline structure observed at the lowest temperature [13].Let us turn now to the correlation between the observed crystalline structure and the flat DP appearing at T = 205K.Usual surface crystallisation studies predict the formation of solid planes, parallel to the surface, accompanied by sharp and narrow peaks in the DP across the interface.As can be observed from figure 3, molecular simulation reveals the formation of crystalline planes, parallel to each other, whose axes are preferentially titled by several degrees from the surface normal.Such a result has been reported in early observations for liquid n-alkanes.In addition, as can be seen in the inset of figure 3, atoms in the titled crystalline planes are quasi hexagonally packed, in good agreement with previous observations in surface crystallization [14].It turns out that a flat DP should be the signature of tilted crystalline planes in the inner region of the slab.

Conclusion
In this paper, we have investigated the free surface crystallisation of l-Hg using molecular dynamics simulations and employing an accurate pair interaction model.Our results may be summarized as follows.As temperature decreases, some usual transitions (liquid/amorphous and amorphous/solid) are observed.The structural properties predicted by our simulations are corroborated and supported by previous works.At temperatures below the melting point, a flat DP has been shown to reveal a posteriori the existence of crystalline planes, parallel to each other, in which atoms are quasi hexagonally packed, that are tilted from the surface normal.Further work along the lines of this work will include the interactions changes across the interface.

Fig. 3 .
Fig. 3. Typical snapshot of the system at T = 205K.In the inset is shown the in-plane atoms arrangement.