Possible phase transition in liquid Cesium at ambient pressure

TL;DR

This study uses molecular dynamics to investigate whether liquid cesium undergoes a phase transition at 590 K under ambient pressure, finding no evidence of such a transition despite some thermodynamic anomalies.

Contribution

The paper provides the first molecular dynamics analysis of liquid cesium at ambient pressure around 590 K to explore potential phase transitions.

Findings

Small thermodynamic kinks observed

No evidence of liquid-liquid phase transition

Thermodynamic anomalies do not indicate a phase change

Abstract

We report a molecular dynamics study of liquid cesium at ambient pressure intended to check the possibility of liquid-liquid phase transformation at $T=590$ K. We find the presence of small kinks on thermodynamic characteristics of the system, but no phase transition.

Equations

$$\gamma =\frac{1}{T}{\left(\frac{\partial T}{\partial P}\right)}_S.$$

Full text

I Introduction

Cesium is an element which demonstrates extremely complex behavior. The phase diagram of cesium contains several structural transitions in crystalline region tonkov . At low pressures it forms an bcc phase (Cs-I) which transforms into fcc (Cs-II) under pressure. Further pressurizing induces a transition into Cs-III phase which has very narrow region of stability. More structural transitions take place in Cs if the pressure is increased further.

Interestingly, both Cs-I and Cs-II phases demonstrate maxima at the melting line. It is widely believed that maxima on the melting line evidence of the complex behavior of the liquid phase. Moreover, from the common sense one expects that close to the melting line the local structure of liquids should be similar to the corresponding crystal. That is why it is possible that structural transformations in the solid part of the phase diagram could induce also transformations in the liquid.

Such changes in the properties of the liquid Cs were discovered several decades ago. For example, in Ref. exp-1990 the structure of liquid cesium at several pressures and temperatures slightly above the melting line maxima was measured by X-ray diffraction. The authors discovered even two transformations. The first one at the pressure $P_{1}\approx 2$ GPa corresponds to changes from bcc-like liquid to the fcc-like one. The second transition occurs at $P_{2}\approx 3$ GPa. This transition transforms the fcc-like liquid into even denser one. Further evidence of the liquid-liquid phase transition (LLPT) in cesium was reported in falconi . X-ray measurements of structure of liquid cesium at high pressure and high temperature were performed and it was found that the density demonstrated a jump at $P=3.9$ GPa. However, in a recent work hattori it is claimed that the procedure for measuring the density in Ref. falconi is not enough precise and more accurate determination of the density demonstrates no discontinuity.

Phase transitions in liquid cesium under pressure were studied in a large body of publications. However, all these publications were concerned to the transitions mentioned above. At the same time several years ago several publications reported the possibility of another phase transition in liquid cesium blag-97 ; blag-2000 ; blag-2008 ; blag-2012 . The authors of these publications performed measurements of adiabatic thermal pressure coefficient of liquid cesium at the ambient pressure and the temperatures up to 600 K. The adiabatic thermal pressure coefficient is defined as:

[TABLE]

From thermodynamic relations it can be shown that $\gamma=\frac{\alpha_{P}}{\rho c_{P}}$, where $\alpha_{P}$ is the thermal expansion coefficient and $c_{P}$ is the isobaric heat capacity. A small jump of $\gamma$ at the temperature $T_{0}=590$ K was observed in blag-97 ; blag-2000 ; blag-2008 ; blag-2012 . Although this jump is just about $2\%$, the authors claim that it is within the experimental uncertainties and make a conclusion that a phase transformation takes place.

An important point about this transition is that it takes place at ambient pressure and relatively low temperatures. Typically LLPT in metals are related to electronic transformations. This explanation was also used for explanation of possible LLPT in Ce at high pressure. It is not clear whether electronic transformation can take place at ambient pressure and temperature about 600 K. However, LLPT can appear even in simple systems of particles interacting via isotropic pair potentials (see, for instance, boat ; genvdw ; skib ; skib1 ; rs-pre ; rs-jetp ). In this case LLPT is induced by existence of two characteristic distances in the interaction potential. Therefore, even at the conditions of experiments from Refs. blag-97 ; blag-2000 ; blag-2008 ; blag-2012 LLPT appears as a transition from BCC-like liquid to FCC-like one. However, up to now only one group reported observation of this phase transition at $T=590$ K and independent measurements would be able to shed light on the problem. There is also a theoretical work ghatee which pretends to describe the transition in liquid Cs at ambient pressure by employing Lennard-Jones potential with parameters which depend on the temperature. However, this model was constructed specially to describe the results of experiments on Refs. blag-97 ; blag-2000 ; blag-2008 ; blag-2012 and because of this it cannot be considered as independent verification of the results of these experiments.

In this paper we for the first time perform a simulational study of the possible phase transition in liquid cesium at $T=590$ K and ambient pressure in order to carry out an independent exploration of the phenomenon by a different approach.

II System and Methods

In the present study we simulated a system of 3456 particles in a cubic box with periodic boundary conditions by means of molecular dynamics method. The interaction of the particles is described by embedded atom potential proposed by Belashchenko in Ref. eam-cs . This potential was specially designed for molecular simulation of liquid cesium. In particular, in Ref. eam-cs a possibility of LLPT at $3.9$ GPa was questioned and no indications of the transformation were observed. The initial configuration of the particles was taken to be bcc lattice. We studied the temperatures from $T_{min}=300$ K up to $T_{max}=1200$ K. At all these temperatures the system rapidly melts. The system was simulated at constant pressure (isothermal-isobaric ensemble). The pressure was set to 1 bar. Additional simulations at 100 and 1000 bars were also performed. The time step was set to 0.2 fs. The equilibration run took $5\cdot 10^{6}$ steps and the production run - $10\cdot 10^{6}$ steps. We calculated the equations of state of the system and the internal energies. The response functions (isobaric heat capacity $c_{P}$ and thermal expansion coefficient $\alpha_{P}$) were obtained by numerical differentiation of these functions. We also performed the calculations of structural properties (radial distribution functions, structure factors and angle distribution functions) and diffusion coefficients. The later were computed via the Einstein relation.

III Results and Discussion

First of all we calculate the equation of state along an isobar $P=1.0$ bar from $T_{min}=300$ K up to $T_{max}=1200$ K. We also performed calculations of enthalpy along this isobar and computed the isobaric heat capacity $c_{P}$ by numerical differentiation. The coefficient $\gamma$ was obtained as $\gamma=\frac{\alpha_{P}}{\rho c_{P}}$, where the thermal expansion coefficient $\alpha_{P}$ was obtained by numerical differentiation of the density $\alpha_{P}=-\frac{1}{\rho}\left(\frac{\partial\rho}{\partial T}\right)_{P}$. The results for $\alpha_{P}$, $c_{P}$ and $\gamma$ are given in Figs.  1 (a)-(c). One can see that in spite of rather good averaging the data are still noisy. The thermal expansion coefficient demonstrates two kinks at $T=620$ and $T=690$ K (Figs.  1 (b)). Also a small maximum and minimum are observed at $T=900$ and $1100$ K respectfully. The isobaric heat capacity looks smoother although it is also noisy in the interval $580-720$ K (Fig.  1 (c)).

The behavior of the coefficient $\gamma$ is mostly governed by the behavior of the thermal expansion coefficient (Fig.  1 (b)). It also demonstrates two small peaks at $T=620$ and $T=690$ K which are slightly above the assumed transition point.

A comparison with experimental values of $\gamma$ is also given in Fig.  1 (a). One can see that the data from simulations are systematically lower. However, the kinks at the $\gamma$ from simulations almost coincide with the experimental point ($T=620$ K from simulation vs. $T=590$ K from experiment) which allows us to assume that the simulation shows qualitative agrement with experiment.

From Fig.  1 (a) we cannot make a clear conclusion about the presence of transition. Although we observe some kinks in the adiabatic thermal coefficient they are just a bit above the range of the errors of calculations. However, the effect observed in experiments is also hardly visible. In order to see whether there is a phase transition in the system we study its properties on both sides of the possible transition.

We study the equations of state along isotherms (Fig. 2) and we do not observe anything which can be interpreted as a phase transition.

As the next step we study radial distribution functions and structure factors of the system (Fig. 3 (a) and (b)) and we again do not observe any significant changes in the structure. Some less pronounced structural changes can be observed from the triplet correlation function which is related to the angle distribution function, shown in Fig. 3 (c). However, this function also does not demonstrate any changes by passing the hypothetical transition point.

Phase transformations are usually accompanied by jumps of the nearest neighbor number (NN). The location of the first peak of radial distribution function also can experience a jump. The coordination number can be obtained as $NN=4\pi\rho_{N}\int_{0}^{r_{m}}r^{2}g(r)dr$, where $\rho_{N}=N/V$ is the number density of the system and $r_{m}$ is the location of the first minimum of $g(r)$. Fig. 4 shows the number of nearest neighbors at $P=1$ bar isobar. One can see that the NN is rather noisy. One can find some kinks next to the kinks of $\gamma$, however, they are within the accuracy of calculations. The position of the first peak of $g(r)$ changes smoothly with temperature.

We apply one more criterion which can establish the presence of the phase transforation in the case of phase transition one can expect that the mobility of particles can change which would lead to the changes in diffusion coefficient. The diffusion coefficient as a function of temperature along $P=1$ bar isobar is shown in Fig. 5. Again we do not observe any changes in its behavior. The temperature dependence of diffusion is smooth and does not indicate any kind of phase transition.

Finally, we study the behavior of the adiabatic thermal pressure coefficient along several isobars in order to see if it demonstrates some kind of jumps at other pressures. The results are shown in Fig. 6. One can see that at all studied pressures the curves are rather noisy which is, probably, related with numerical differentiation of the results of the simulations and does not allow to make any firm conclusion. Although we have checked that our simulation converged, there is always systematic noise in MD simulation which can be reduced by increasing the averaging. In order to reduce the noise we perform additional simulations. The simulation setup was the same as it was described in the methods section. Some of the additional simulations were performed starting from previously obtained configurations. Others were started from different initial positions. In total 10 different simulations at each point were performed and the averaging of the quantities of interest was done. The results are shown in Fig. 7. One can see that the curve becomes smoother, however, the kinks at $T=620$ and $690$ K preserved. Therefore, these kinks can be real feature of the system under investigation rather then numerical error, however, the physical reason of these kinks is unclear.

In conclusion, we perform molecular dynamics study of probable liquid-liquid phase transition in liquid cesium. The results of the simulations do not show any traces of phase transition in such quantities as equation of state, radial distribution function, structure factor, angle distribution function and diffusion coefficient. However, the adiabatic thermal pressure coefficient demonstrates two small peaks which cannot be eliminated by increasing the averaging. The origin of these peaks is not clear and requires further studies.

We thank D.K. Belaschenko and R.M. Khusnutdinoff for help with EAM potential construction and V.V. Brazhkin for useful discussions. This work was carried out using computing resources of the federal collective usage center ”Complex for simulation and data processing for mega-science facilities” at NRC ”Kurchatov Institute”, http://ckp.nrcki.ru, and supercomputers at Joint Supercomputer Center of the Russian Academy of Sciences (JSCC RAS). The work was supported by the Russian Foundation of Basic Research (Grant No 18-02-00981).