Isomorphs in the phase diagram of a model liquid without inverse power law repulsion

It is demonstrated by molecular dynamics simulations that liquids interacting via the Buckingham potential are strongly correlating, i.e., have regions of their phase diagram where constant-volume equilibrium fluctuations in the virial and potential energy are strongly correlated. A binary Buckingham liquid is cooled to a viscous phase and shown to have isomorphs, which are curves in the phase diagram along which structure and dynamics in appropriate units are invariant to a good approximation. To test this, the radial distribution function, and both the incoherent and coherent intermediate scattering function are calculated. The results are shown to reflect a hidden scale invariance; despite its exponential repulsion the Buckingham potential is well approximated by an inverse power-law plus a linear term in the region of the first peak of the radial distribution function. As a consequence the dynamics of the viscous Buckingham liquid is mimicked by a corresponding model with purely repulsive inverse-power-law interactions. The results presented here closely resemble earlier results for Lennard-Jones type liquids, demonstrating that the existence of strong correlations and isomorphs does not depend critically on the mathematical form of the repulsion being an inverse power law.


Introduction
Recently a series of papers has been published concerning so-called strongly correlating liquids and their physical properties [1][2][3][4][5][6]. Liquids that exhibit these strong correlations have simpler thermodynamic, structural, and dynamical properties than liquids in general. A strongly correlating liquid is identified by looking at the correlation coefficient of the equilibrium fluctuations of the potential energy U (r 1 , . . . , r N ) and virial W (r 1 , . . . , r N ) ≡ −1/3 i r i · ∇ ri U (r 1 , . . . , r N ) [7] at constant volume: Here brackets denote averages in the NVT ensemble (fixed particle number, volume, and temperature), Δ denotes the difference from the average. The virial W gives the configurational part of the pressure [7], pV = N k B T (p 1 , . . . , p N ) + W (r 1 , . . . , r N ).
Strongly correlating liquids are defined [1] as liquids that have R ≥ 0.9. The origin of strong W U correlations was investigated in detail in references [3,4] for systems interacting via the a e-mail: tbs@ruc.dk Lennard-Jones (LJ) potential: The fluctuations of W and U are dominated by fluctuations of pair distances within the first neighbor shell, where the LJ potential is well approximated by an extended inverse power law (eIPL), defined as an inverse power law (IPL) plus a linear term [3]: The IPL term gives perfect U W correlations, whereas the linear term contributes little to the fluctuations at constant volume: when one pair distance increases, others decrease, keeping the contributions from the linear term almost constant (this cancellation is exact in one dimension). The consequence is that LJ systems inherit some of the scaling properties of the IPL potential -they have a "hidden scale invariance" [4,8]. Prominent among the properties of strongly correlating liquids is that they have "isomorphs", i.e., curves in the phase diagram along which structure, dynamics, and some thermodynamical properties are invariant in appropriate units [5,6]. The physics of strongly correlating liquids was briefly reviewed recently in reference [9]. Since the LJ system consists of two IPL terms, it is perhaps tempting to assume that a repulsive (inverse) power law is necessary for the hidden scale invariance described above. In the present paper we use the modified Buckingham (exp-six) pair potential to show that this is not the case. The Buckingham potential was first derived by Slater from first-principle calculations of the force between helium atoms [10]. Buckingham later used this form of the potential to calculate the equation of state for different noble gases [11]. The Buckingham potential has an exponential repulsive term, while the attractive part is given by a power law [12,13]: Here is the depth of the potential well and r m specifies the position of the potential minimum. The parameter α determines the shape of the potential well. The Buckingham potential is better able to reproduce experimental data of inert gasses than the LJ potential [14][15][16], but is also computationally more expensive (unless lookup tables are utilized [7]).
All simulation data in this paper were obtained from molecular dynamics in the NVT ensemble. The samples contained 1000 particles. The simulations were set up by instant cooling from a high temperature state point followed by an equilibration period, to ensure the simulations were independent from each other. The simulations were performed with the RUMD molecular dynamics package [17], which is optimized for doing computations on state-of-the-art GPU hardware.

Correlations in single-component Buckingham liquids
To compare the simulations with experiments [18], argon parameters from reference [14] were used; α = 14.0, r m = 0.3866 nm, /k B = 123.2 K. As can be seen in Figure 1a, the single-component Buckingham (SCB) liquid is strongly correlating (R ≥ 0.9) in parts of the phase diagram, particularly at high densities and/or temperatures. The correlation coefficients (Eq. (1)) of the Buckingham systems are very similar to those of argon and the LJ system (dotted line in Fig. 1). This is a first indication that the actual functional form of the repulsive part of the potential does not have to be an inverse power law in order for a system to exhibit strong W U correlations.
Another interesting property of the fluctuations is the quantity γ defined [4,5] as When a system is strongly correlating (R is close to one), ΔW ≈ γΔU . For IPL potentials γ is constant and equal to n/3 and R = 1. For non-IPL potentials, however, R < 1 and γ may change with temperature and density as seen in Figure 1b [1,19]. Especially for R < 0.9, we find γ changing rapidly. The curves are similar for the 20.0 mol/L SCB and SCLJ systems, except for a vertical offset.  [14], and R was calculated directly from equation (1). For argon R was calculated from experimental data [18] as described in references [1,3]. The correlations are strongest for state points with both high density and high temperature, and the difference between the Buckingham and the LJ potential is small. The correlation coefficient R > 1 for low-temperature 20.0 mol/L argon is of course unphysical and either caused by an uncertainty in the experimental data or the approximations applied in the calculation of R (see Refs. [3,20] for details). (b) The value of γ (Eq. (6)) plotted versus temperature for the same systems as in (a). For argon γ was calculated from experimental data [18] using equation (7). γ decreases slowly for increasing temperatures, except when the correlation coefficient is low (R 0.9).
Fluctuations in U and W are of course only directly accessible in simulations. For experimental systems one must revert to the use of thermodynamic quantities that reflect the fluctuations in U and W . For instance, the configurational part of the pressure coefficient β ex V = (∂(W/V )/∂T ) V and the configurational part of the isochoric specific heat per unit volume c ex V = (∂(U/V )/∂T ) V can be used to to calculate γ for argon as follows [1,3,5]: The values of γ for argon obtained in this way are plotted in Figure 1b, and the agreement with the Buckingham systems is good. This confirms that the Buckingham potential produces more accurate predictions of experimental argon data than the LJ potential. Interestingly, low density argon has a higher correlation coefficient than high density argon. This is the opposite of what is found for the Buckingham and the LJ potentials. Furthermore, the buckingham data are in better agreement with the argon data at low density than at high density. At the present we do not have any explanation for these observations.

Isomorphs in binary Buckingham mixtures
Strongly correlating liquids are predicted to have isomorphs, which are curves in the phase diagram along which structure, dynamics, and some thermodynamical properties are invariant in appropriate reduced units [5,6]. Introducing reduced coordinates asr i = ρ 1/3 r i , two state points (1) and (2) are defined to be isomorphic if pairs of microscopic configurations with same reduced coordinates (r Here the constant C 12 depends only on the two state points and equation (8) is required to hold to a good approximation for all physically relevant configurations [6]. An isomorph is a curve in the phase diagram for which all points are isomorphic (an isomorph is a mathematical equivalence class of isomorphic state points). The isomorphic invariance of structure, dynamics, and some thermodynamical properties -all in reduced units -can be derived directly from equation (8) [5]. Only IPL liquids have exact isomorphs, but it has been shown that all strongly correlating liquids have isomorphs to a good approximation (Appendix A of Ref. [5]).
Among the thermodynamical properties that are isomorphic invariant is the excess entropy, S ex ≡ S − S ideal , where S ideal is the entropy of an ideal gas at the same temperature and density. In the following, isomorphic state points were generated by utilizing that the quantity γ in equation (6) can be used to change density and temperature while keeping the excess entropy constant [5,6]: By choosing the density of a new isomorphic state point close to the density of the previous isomorphic state point, the temperature of the new state point can be calculated from the fluctuations by combining equations (6) and (9) [5]. In this way a set of isomorphic points can be obtained from one initial state point. The predicted isomorphic invariance of the dynamics is most striking in viscous liquids, where the dynamics in general depend strongly on temperature and density. To demonstrate that systems interacting via the Buckingham potential have isomorphs, we study what we term a Kob-Andersen binary Buckingham (KABB) mixture with potential parameters being the same as for the original Kob-Andersen binary LJ (KABLJ) mixture [21]: AA = 1.0, r m, AA = 6 √ 2, AB = 1.5, r m, AB = 0.8 6 √ 2, BB = 0.5, r m, BB = 0.88 6 √ 2. A 4:1 mixture (A:B) was used with α = 14.5. The potentials were truncated and shifted at r cut ij = 2.5r m, ij / 6 √ 2. One of the predicted invariants on an isomorph is the structure of the system. To test this prediction, the radial distribution function in reduced coordinates g(r) = g(ρ 1/3 r) was plotted for isomorphic state points (Fig. 2a). 1 1.  for simulations of the KABB mixture. Both graphs are in reduced units wherer = ρ 1/3 r. (a) g(r) for isomorphic state points and the three different particle combinations. The structure is invariant on the isomorph for the AA particle pairs, but for the AB and BB pairs the structure is less invariant. (b) g(r) for isothermal state points of smaller density variation. The structure is not invariant on the isotherm for any of the particle pairs.
The structure is invariant for the large (A) particle pair correlation function to a very good approximation. For the AB and BB pairs the structure is less invariant. However, when a comparison is made with Figure 2b, it is clear that g(r) for the AB and BB pairs is still more invariant on an isomorph than on an isotherm (note that the density variation on the isomorph is larger than on the isotherm). This situation is similar to what is found for the KABLJ system [5].
To investigate the dynamics of the systems, the incoherent intermediate scattering function F s (q, t) is plotted in reduced units in Figures 3a and 3b. The presence of a plateau in F s shows that the system is in a viscous state,  where the dynamics are highly state point dependent. The large difference in F s for the two isothermal state points confirms this (dashed lines). For the isomorph all F s data collapse more or less onto the same curve, showing that the dynamics are indeed invariant to a good approximation on an isomorph. In contrast to the radial distribution functions, the invariance holds well for both types of particles.
To investigate the invariance in dynamics further, the coherent intermediate scattering function was calculated (Fig. 4). The coherent intermediate scattering function was calculated from the spatial transform of the number density ρ(q) [7]. In order to obtain good results, it is necessary to average over time scales that are 10-15 times longer than what is usual for the intermediate scattering function. This is the reason that there are less state points shown for the coherent-, than for the incoherent intermediate scattering function. The data confirm that the dynamics are invariant on the isomorph, especially when compared to the isothermal density change (dashed lines). However, the invariance seems to hold slightly better for the AB and BB parts, which is the opposite of what is seen for the structural invariance.
For systems described by a generalized LJ potential consisting of two IPL terms, the invariance of the structure  The solid lines represent isomorphic state points, while dotted lines show isothermal density changes for comparison. Again, the q-vector is kept constant in reduced units: qAA = 7.34(ρ/1.2) 1/3 , qAB = 8.01(ρ/1.2) 1/3 and qBB = 6.01(ρ/1.2) 1/3 . Also F (q,t) is invariant on the isomorph. Contrary to what is seen for g(r), the invariance holds better for the AB and BB parts.
leads to a prediction for the shape of an isomorph when plotted in the U -W plane [6] (generalized LJ potentials are a sum of inverse power laws). Since the repulsive term in the Buckingham potential is described by an exponential function, it is not possible to derive an exact equation that describes the isomorph in terms of U and W . Figure 5a shows that isomorphs for the KABB system agree well with the prediction for the 12-6 LJ system if α = 14.5 (this value of α was chosen to demonstrate this feature). For α = 13.0, there is a significant difference with the predicted shape at higher density and temperature. Figure 5b shows the isomorphs for both values of α after scaling U and W by the same isomorph-dependent factor, demonstrating the existence of a master isomorph [6]. This shows that master isomorphs exist not only in generalized LJ systems where they can be justified from analytical arguments [6]. (Color online) (a) Plot of the potential energy per particle versus virial per particle for the KABB system. The solid lines are the predictions of the isomorph shape for the 12-6 LJ potential [6]. For these predictions the initial state points with ρ = 1.2 were used as reference point. Since a new value of γ was calculated for every state point, γ is not constant on the isomorphs, but changes approximately 10% along the isomorphs. The shape of the KABB isomorph agrees very well with the predicted shape for the 12-6 LJ potential for α = 14.5. For α = 13.0, the shape is different. (b) The same data now scaled with W * 0 defined as the virial at U = 0 [6]. The isomorphs scale onto each other, forming a so called master isomorph for each value of α.

The inverse-power-law (IPL) approximation
As mentioned in the introduction, a generic explanation [3][4][5] for the existence of strong correlations and isomorphs in non-IPL systems, is the fact that some pair potentials can be well approximated by an eIPL (Eq. (4)) as shown in Figure 6. Putting this explanation to the test, it was recently demonstrated that structure and dynamics of the KABLJ system can be reproduced by a purely  Figure 7. The difference of the IPL approximation and the Buckingham potential (dashed green line) is more or less linear in the first peak of g(r). By subtracting this linear term from the IPL term the eIPL approximation is found (dashed blue line).
repulsive IPL system even in the viscous phase [22]. In the following we demonstrate that this procedure works also for the KABB system, despite its non-IPL repulsion. Following Pedersen et al. [22], we assume that the Kob-Andersen IPL (KABIPL) system used to approximate the KABB system has the form where the parameters ij and σ ij are the Kob-Andersen parameters for the different types of particles and the constants A and n are independent of particle type. For IPL liquids it is known that W = (n/3)U , so in principle the value of n could be calculated from γ determined from the W U fluctuations (Eq. (9)). For non-IPL liquids however, there is a slight state point dependence of γ, so instead the slope of an isochore was used to determine n (Fig. 7a) making use of the identity [5] For the Buckingham potential with α = 14.5 we obtained γ = 4.904 and n = 14.71. This is lower than the γ = 5.16 which was found for the 12-6 LJ potential [22]. This is also consistent with the data in Figure 1b where the SCB system has a lower value of gamma than the SCLJ system. From equation (10) it follows that the total internal energy of the IPL system can be written as The method used to find the value of A of equation (12). U IP L = i>j ij (σij /rij ) n was calculated from configurations drawn from the KABB simulations and plotted against the energy obtained during the simulations; the value of A was then obtained from the slope of the mean energies (again marked by yellow crosses).
The scaling factor A was determined from the slope of the mean values of the energies in a U , U IP L plot (Fig. 7b), where U IP L is given by equation (12) evaluated on configurations from simulations of the KABB mixture [22]. Using these parameters, simulations of the KABIPL systems were performed and the results were compared with the results of the KABB system. In Figure 8 the incoherent intermediate scattering function of the two systems is plotted for comparison. The KABIPL reproduces the dynamics of the KABB system very well. It should however be noted that in spite of the good reproduction of the dynamics, the KABIPL had a stronger tendency to crystallize than the KABB system at the two lowest temperatures due to the absence of attractive forces. The good agreement shown in Figure 8 only holds if both systems are in the same (supercooled) state. From the fluctuations in the potential energy one can calculate the excess isochoric specific heat using [7]: In Figure 9a C ex V is plotted for different isochores calculated from KABB and KABIPL simulations. The heat capacities for the two systems follow each other closely, although there is a small and systematic difference increasing with density. This is similar to what was found for the KABLJ system [22], but the deviations are slightly larger for the KABB mixture. Figure 9b shows that the excess heat capacity to a good approximation obeys density scaling, C ex V /N = f (ρ γ /T ), and Rosenfeld-Tarazona scaling, The configurational part of the intensive isochoric specific heat C ex V /N as a function of temperature. Three isochores were simulated of the KABB and the KABIPL systems. At low density the agreement between the two systems is fairly good, but for higher densities the differences become larger. (b) The same data plotted versus T /ρ γ where γ = 4.904. The data collapse on a single curve, which shows that density scaling works. The function (1.42ρ γ /T ) 2/5 was fitted to the data (dashed line), showing that Rosenfeld-Tarazona scaling is also obeyed. C ex V /N = g(ρ)T −2/5 [23] -again in good agreement with results for the KABLJ system [22].

Conclusion
The Buckingham potential has been shown to be strongly correlating like the Lennard-Jones potential. In spite of its exponential repulsion, the Buckingham potential's dynamics and heat capacity can be closely approximated by a purely repulsive IPL system. In particular the system has good isomorphs in the phase diagram. These findings are very similar to those found for Lennard-Jones systems. We conclude that the existence of strong correlations and isomorphs is not dependent on the repulsion being an inverse power-law.