Computer simulation and modeling the metal to insulating transition of liquid mercury via pair, empirical, and many-body potentials

In this study, Monte Carlo simulations were carried out to detect the metal–insulator transition in liquid mercury by changing the static properties at the microscopic scale. In the simulations pair, empirical, and many-body potentials govern metal–metal interactions in the canonical ensemble. The structural and thermodynamic changes over a wide temperature range from 773 to 1773 K are described in the first coordination shell and residual internal energy, respectively. The results reveal that during the simulated heating process, the properties undergo significant change at 1673 K, which is in connection with the metal-nonmetal transition in the liquid. These findings coincide with the experimental observations of this thermodynamic phenomenon. Notably, the free energy of association that renders the system to this thermodynamic state is estimated.


Introduction
As a unique feature of mercury (Hg) metal, it is liquid at room temperature. Despite its toxicity, Hg has found a wide range of applications in medicine and industry which is, for a large part, due to its ability of dissolving nearly all metals. Examples include dental amalgams [1], energy production [2], mercury-based self-assembled monolayers (SAMs) as semiconductor devices [2,3], or liquid Hg as a base for Langmuir films of imidazolium-based ionic liquids for electrochemical applications [4]. On the other hand, the removal of mercury as a contamination source has become a task of high environmental relevance in recent years, advanced by both experimental and computational studies [5][6][7][8][9]. To support these studies, further explorations of the properties and the structural phase transition of liquid Hg are needed, with special attention to the phase transition from liquid metal to liquid semiconductor or liquid insulator.
Due to related technical difficulties in dealing with experimental measurement, especially to the high temperatures, computer simulation currently is a powerful and convenient tool to study the properties at the microscopic level. It can remedy of experimental information.
This simulation study focused on the metal to nonmetal transition of the expanded liquid Hg in the framework of thermodynamic and structural properties. Hence, a temperature in the range of 773-1773 K consistent with their experimental densities is preferred.
The metal-nonmetal transition (M-NM) is a fascinating phenomenon that was observed firstly, in 1940 for fluid mercury [10,11]. Since then, researchers have constantly aimed for detecting of this phenomenon in both theoretical and experimental approaches. Briefly, the determination of the transition was probed by different methodologies based on a change in special properties. For instance, static properties; structure factor, S(Q), and pair correlation function, g (r), using classical MD [12] or ab initio MD [13,14], Reverse Monte Carlo (RMC) on X-ray diffraction of S(Q) [15], properties change through the equations of state [16][17][18][19], as well as dynamic properties; diffusion coefficient, shear viscosity by ab initio MD [14].
Notably, the conclusion of all these different methods, theoretical or practical confirmed that the M-NM transition occurred at a density nearness of 9 g/cm 3 . In the vicinity of this density, the transition appeared within the alteration of aforesaid properties.
In this study, Monte Carlo (MC) simulations were conducted to study the M-NM transition in liquid mercury. Changing in the structure and thermodynamic properties specifically coordination number and excess internal energy is used to acquire this transition.
Remarkably, probing this phenomenon through pair, empirical, and many-body potentials herein is for the first time.

Modeling and simulation
Monte Carlo (MC) simulation was performed in the NVT ensemble to probe the phase transition of metal to nonmetal of liquid mercury. The MC method is based on the Metropolis algorithm [25], in which the positions of the atoms are allowed to change between the old and new state [25,26]. This procedure is continued until equilibration is achieved by the convergence of the energy.
The cubic cell of side L corresponding to the experimental density ρ at the given temperature was chosen (see Table 1). The periodic boundary conditions were employed in all three directions. The pair potentials were truncated at the cut-off radius 3.5 times of atomic diameter. Of the 4 × 10 7 configurations generated in the simulation, the first half was propagated to reach equilibration and the rest were used for statistical analysis and ensemble average calculations.
The atomic interactions were modeled through the Lennard-Jones (LJ), empirical, and many-body potentials that are concisely described as follows:

A: Lennard-Jones potential model
where r is the distance between two interacting fluid atoms, σ denotes the LJ atom-atom collision diameter, and ε is the LJ fluid-fluid potential well-depth, in which their values are reported elsewhere [30].

B: Embedded atom model
The embedded atom model (EAM) is written as follows: The first term is many-body interaction, which is embedding energy required for the placement of i-th atom in the electron cloud density (ρ i ) that is made by j neighbor atoms at the location of atom i. The second term is pair-wise contribution. More details and a lot of potential parameters of the EAM are stated in the reference [31,32].

C: Empirical potential model
The form of the empirical potential is the following functional equation: The first term is Born-Mayer repulsion, while the latter is the attractive part, in which A 1 is a modeled as a polynomial function of the temperature (T) as: The potential parameters are given by Bretonnet et al. [33].

Results and discussion
In this study, the simulated static properties namely configurational internal energy and the first coordination number of the liquid mercury at different thermodynamic states are utilized for examining the M-NM transition that can be seen in the below section.

Static properties change
It should be noted that EAM is a complicated potential function in which possesses a lot of parameters in both pair and embedded parts. Therefore, initially to the accuracy, the program code that was used here, the properties such as excess internal energy (U exe ) at 773 K are reproduced from reference [31]. It is − 52.15335 kJ/mol that is the same as the value of − 52.16008 kJ/mol via the MD method in this reference [31]. The U exe from this reference is calculated by internal energy (U), i.e., U exe = U − 1.5 RT, in which R is the universal gas constant. As well as, this reproduced scheme is used for the empirical potential (C), for instance, the pair correlation function, g (r), at 773 K is duplicated from reference [33] (supplied by Bomont), as shown in Fig. 1. As seen, all g (r) features namely the peaks position and their heights are exact matching in the comparison. Then, the simulations are conducted and continued for the aim of this study for detecting how the thermodynamic phenomenon of metal-nonmetal transition affects on these properties at the atomic level.
The structural change of liquid Hg is depicted by the first coordination shell (Z) or liquid microstructure that it estimated using integration on pair correlation function, g (r) [34,35], as follows: where r p is the first peak position of g (r) (3.05 Ǻ) and ρ is the number of density.
The calculated Z values as a function of temperature (Z-T), using pair, empirical, and many-body potentials are illustrated in Fig. 2. Firstly, the overall feature of this curve coincides with previous studies [13,15,[36][37][38] and manifests there is an inhomogeneous expansion. Because, the first peak position in g (r) as the nearest neighbor distance based on three model potentials (A, B, C) mainly remains constant, while the trend of coordination number turns down.
The figure of Z-T indicates that at 1673 K (corresponding to the density of 9.25 g/cm 3 ), the slope of the curve or the rate of change in the decreasing of Z on two sides of this temperature is evidently different. The slope of the curve for T ≤ 1673 K (the range of 773-1673 K) slightly declines, while for T ≥ 1673 K (the variety of 1673-1773 K) significantly decreases. The slopes of the curve along with (4) Z = 2ρ ∫ r ρ 0 g(r)4πr 2 dr these two temperature ranges, cross at 1673 K, as clear in the figure. Consequently, at 1673 K, a distinctive change in the microstructure of the liquid occurs, which is connected to metallic to nonmetallic (M-NM) phenomenon. Hence, this instability in the liquid microstructure is a driving force for altering metal state into insulating ones. On the statistical mechanic language, as thermodynamic states change toward the nonmetallic phase, the local aggregated of Hg atoms (molecular-like) emerges in the system [13], accordingly, the phase transition (M-NM) arises. Therefore for observation of this issue, the average ensemble of neighboring distribution for each Hg atoms is calculated by having geometric locations of atoms in the simulation box, i.e., x 2 + y 2 + z 2 ≤ r p 2 . r p is the radius of the first coordination shell. The comparison of neighboring distribution in the metallic state and nonmetallic through EAM is plotted in Fig. 3. This plot characterizes a significant feature that in the Lines are trend lines nonmetallic state more fraction of atoms' tendency to adapt one neighboring (dimer-like) distribution. Similar curves resulted for the other potential models (A, C) (see Figs S1 and S2 in the Supplementary Information). Although the distributions generated by the three models exhibit sizeable quantitative differences, they all indicate that the nonmetallic state involves a higher dimer-like contribution than the metallic state. In more detail, both B and C (r p = 3.05 Ǻ) exhibit sizeable dimer-like contributions (see Fig. S1), while in the LJ model (A) (r p = 3.1 Ǻ), the dimer-like configuration becomes dominant (see Fig. S2).
Also, the existence of the dimer-like of Hg is experimentally detected [39] and coincides with the neighboring distribution feature in this study. In metallic bond language, the emergence of these dimer-like species leads to localized conduction electrons in the system and hence, insulating or semiconductor state consequences. This issue coincides with the value of free energy association (∆A ass ) between atoms in the system for metallic and nonmetallic states, which is estimated through pair, many-body, and empirical potentials, as listed in Table 2. The ∆A ass is estimated via the following formula [40,41]: where w(r) is the potential of the mean force, and V is a normalization factor for correlation length and it is the volume of the first coordination shell.
As can be seen from the values of ∆A ass in Table 2, the nonmetallic state (1673 K) exposes more positive free energy of association with respect to metallic (773 K) ones that points out the association among atoms in the entire system is declined. As the snapshots (see Fig. 4) show, in the nonmetallic state, the space among atoms is significantly increased with respect to the metallic state. Hence, the lack of sufficient packing of atoms or the absence of effective spatial distances between atoms in the system resulted in the nonmetallic phase.
Additionally, the excess internal energy as a function of temperature (U exe -T) is visualized in Fig. 5. As can be seen, there is a particular feature in U exe -T curve, in which it exposes upward curvature as a thermodynamic state approaching the nonmetallic phase (1673 K). The boundary of changing trends is visualized by the vertical dashed line in this figure. Furthermore, in this figure, a comparison of the results within EAM is made with the evaluable data of MD based on EAM [31]. Moreover, it is not detected yet as the metal to nonmetal phase evolutions, how the excess internal energy alters. It should be noted that in this figure, the significant change in excess energy that resulted from empirical potential (C) in the temperature range of 773-1473 K is not correlated to the phase transition. This characteristic relates to the temperature dependency of parameter A 1 (see Eq. 3) in this potential equation. The empirical potential in the nonmetallic domain (1673-1773 K) points out the same trend with the other potential models (A, B). Therefore, three model potentials in spite of their different numerical values due to their natures forecast the nonmetallic range with the same trend. In short, the MC simulation of this study suggests metal to the insulating transition occurs in 9.25 g/cm 3 . This result is in good agreement with both theoretical and experimental studies [14,21,22,37] that the transition is predicted around the density of 9 g/cm 3 .
This study differs from previous publications on the same topic in its analysis of the metal-nonmetal transition (M-NM) in terms of the coordination number, the excess internal energy, and the neighbor distribution in the thermodynamic states. Here, a computer experiment is performed where the liquid system expands (increasing temperature and decreasing the corresponding density), within this process the thermodynamic state (ρ, T) = (9.25 g/cm 3 , 1673 K) of the NM transition is modeled. The certain thermodynamic state is not pictured in earlier studies and a range of densities have been estimated for predicting NM transition. For instance, in the study by Delhommelle et al. [42], the nonmetallic state is modeled over a range of low densities, while the metallic state is modeled over a range of high densities only for one temperature point only (1473 K) by use of the empirical potential (C) via the Wang-Landau scheme.
In the reference [28], the depth of the minimum pair part of the EMA (B) is changed as analytically in the range 1273 K < T < 1803 K (10.98-6.8 g/cm 3 ) that could be related to the M-NM transition.

Conclusions
Metallic to the insulating transition of liquid mercury is investigated in a canonical ensemble by Monte Carlo simulation in the context of pair, empirical, and many-body interaction potentials.
The temperature dependencies of static properties are utilized to acquire the metal to nonmetallic transformation. The simulation results indicate that the first coordination number and configurational internal energy, differ markedly in their trend and features, at the temperature of 1673 K (the density of 9.25 g/cm 3 ). Therefore, in this thermodynamic state, a phase transition from metallic state to insulating occurs, which is consistent with experimental studies. Especially, the simulated neighboring distribution estimates more fraction of atoms' tendency to form a dimer-like species in a nonmetallic state.