Structural behavior of supercritical fluids under confinement
Kanka Ghosh C. V. Krishnamurthy
cond-mat.stat-mech








 kankaghosh@physics.iitm.ac.in.

 cvkm@iitm.ac.in

Department of Physics, Indian Institute of Technology Madras, Chennai-600036, India  
  




 

 
 


 



 
 





             

The existence of the Frenkel line in the supercritical regime of a Lennard-Jones fluid shown through molecular dynamics (MD) simulations initially and later corroborated by experiments on Argon opens up possibilities of understanding the structure and dynamics of supercritical fluids in general and of the Frenkel line in particular. The location of the Frenkel line, which demarcates two distinct physical states, liquidlike and gaslike within the supercritical regime, has been established through MD simulations of the velocity auto-correlation (VACF) and Radial distribution Function (RDF). We, in this article, explore the changes in the structural features of supercritical LJ-fluid under partial confinement using atomistic walls for the first time. The study is carried out across the Frenkel line through a series of MD simulations considering a set of thermodynamics states in the supercritical regime (P = 5000 bar, 240K ⩽ T ⩽ 1500K) of Argon well above the critical point. Confinement is partial, with atomistic walls located normal to z and extending to "infinity" along the x and y directions. In the "liquidlike" regime of the supercritical phase, particles are found to be distributed in distinct layers along the z-axis with layer spacing less than one atomic diameter and the lateral RDF showing amorphous-like structure for specific spacings (packing frustration),  and non amorphous-like structure for other spacings. Increasing the rigidity of the atomistic walls is found to lead to stronger layering and increased structural order. For confinement with reflective walls, layers are found to form with one atomic diameter spacing and the lateral RDF showing close-packed structure for the smaller confinements. Translational order parameter and excess entropy assessment confirms the ordering taking place for atomistic wall and reflective wall confinements. In the "gaslike" regime of the supercritical phase, particle distribution along the spacing and the lateral RDF exhibit features not significantly different from that due to normal gas regime. The heterogeneity across Frenkel line, found to be present both in bulk and confined
systems, might cause the breakdown of the universal scaling between structure and dynamics of fluids necessitating the determination of a unique relationship between them.




  PACS numbers
47.11.Mn, 05.20.Jj, 65.20.De








                             

Structural behaviour of supercritical fluids under confinement
    C.V.Krishnamurthy
    
==============================================================







§ INTRODUCTION

With its high density features like liquids, large diffusivity like gases and excellent dissolving power, supercritical fluids are playing a significant role in purification and extraction processes of various industries <cit.>. Around 40 years back in their seminal paper, M.E.Fisher and B.Widom discussed liquid and gaslike supercritical states by observing the decay behaviour of pair correlation function at large distances using linear continuum and lattice models and challenged the existing description of supercritical fluid as a single homogeneous phase like other states of matter <cit.>. Since then, many experimental studies had been executed to validate the heterogeneous nature of the supercritical fluids <cit.>. G.G.Simeoni et al. <cit.> carried out inelastic X-ray scattering and molecular dynamics simulation to find out a demarcation line between two dynamically different regime ("liquidlike" and "gaslike") in supercritical fluid around critical point called Widom line. In a recent study, however, it has been found that this Widom line doesn't obey the corresponding states principle and the transition lines differ with different fluids <cit.>. Few years back, the discovery of Frenkel line in the phase diagram further adds to the current state of knowledge regarding heterogeneity of the supercritical state of a fluid. It indicates that there exists two distinct phases on either side of the Frenkel line: the "liquidlike" and the "gaslike" states at any arbitrary high temperature and pressures <cit.>. Thus the universal and dynamic Frenkel line is qualitatively different from the Widom line, which unlike Frenkel line, exists near to the critical point only <cit.>. The transition between these two regimes along the Frenkel line occurs when the liquid relaxation time (τ) becomes nearly equal to the Debye vibration period (τ_D), when the system becomes unable to support shear modes<cit.>. Nevertheless, the more convenient approach to detect Frenkel line on the phase diagram from atomistic simulation is to monitor the disappearance of oscillations in Velocity auto-correlation function (VACF) as proposed by V.V.Brazhkin et al. <cit.>.

Structural and thermodynamic properties, associated with Frenkel line crossover, have been deduced from MD simulations in terms of RDF and specific heat capacity <cit.>. D.Bolmatov et al. experimentally proved the presence of thermodynamic boundary associated with Frenkel line from a diffraction experiment on supercritical Argon in a diamond anvil cell (DAC) <cit.>. Extensive investigations have been done on the dynamic crossover of supercritical phases of water <cit.>, Iron <cit.>, CO_2 <cit.>, Argon <cit.>, CH_4 <cit.> etc. 

In a recent review article, J.M.Stubbs covers a wide range of molecular simulation studies of supercritical fluids (SCF) <cit.>. 

Transport and structural behaviour of normal fluids under confinement has been of interest within the physics community due to their unusual properties with respect to the bulk fluid systems. J.Mittal et al. reported a series of studies on self-diffusion, modification of the dynamics and layering of confined hard-sphere fluids through Molecular Dynamics and Monte Carlo simulations <cit.>. Recently, in a study of self-diffusion and radial distribution function of a strongly confined LJ fluid, N de Sousa et al.<cit.> found that in the solid-liquid phase transition region radial distribution function (RDF) corresponding to both the liquid and the solid phases are essentially indistinguishable.

In the present study, we choose atoms interacting with Lennard-Jones(LJ) potential mimicking Argon and consider a set of thermodynamics states in the supercritical regime (P = 5000 bar, 240 K ⩽ T ⩽ 1500 K) of the bulk phase and determine the Frenkel line cross-over through VACF and RDF calculations to be at around T∼ 600-700K. We explore the changes in the structural features across the Frenkel line through a series of MD simulations of supercritical fluid under partial confinement using atomistic walls. Though the structural crossover in bulk supercritical fluid characterized by RDF peaks have been well established <cit.>, the structural behaviour of confined supercritical fluids across the Frenkel line has not been studied yet. 

We use atomistic boundaries on a pair of parallel sides along z direction to simulate the partially confined systems. The simulation domain is taken to be a cuboid. Layering phenomenon under confinement has been observed and studied in detail before crossing , after crossing and in the close vicinity of the Frenkel line. The parallel and perpendicular components of the radial distribution function have been systematically studied for a wide range of confined spacings and structural ordering due to confinement has been understood through pair-excess entropy and translational order parameter calculations normal to the walls. Further, the differences in structural properties of supercritical fluid under both smooth, purely reflective and atomistic walls have been studied. The details of the MD simulation method have been presented in Section.II . Results are discussed in Section.III. Section.IV provides the summary and conclusions. The results of the confinement under purely reflective walls have been discussed in the appendix.



§ MODEL AND SIMULATION DETAILS


We carry out molecular dynamics calculations on LJ fluid using LAMMPS software package <cit.>. To determine the Frenkel line and consistency checks, we model 10^5 particles of LJ fluid fitted to Argon properties (LJ potential, ϵ/k_B=120K, σ=3.4 Å) in a number of  isothermal-isobaric (NPT) ensembles in bulk supercritical phase with periodic boundary conditions imposed along each of the three dimensions. We move on the P-T phase diagram of Argon <cit.> along an isobaric path with a constant pressure of 5000 bar and temperature ranging from 240 K to 1500 K using a Nose Hoover thermostat and barostat (Fig.<ref>). 

We use a cut-off of 20 Å (used previously for Ar in supercritical regime <cit.>) and shift the potential to make the potential and force continuous at the cut-off. In order to shift the potential to zero at the cut-off we have added Δϵ = -9.94×10^-7ev to the potential, which is numerically, too small to affect the critical point significantly. After an energy minimization, standard velocity-verlet algorithm with a time-step (Δt) of 0.0001 picosecond (ps) has been used to equilibrate the system up to 50 ps followed by a 200 ps production run to calculate and analyse the properties of interest and to perform consistency checks.

In our MD simulations for different P,T state points (240K to 700K at 5000bar) the difference between MD and experimental density from NIST database <cit.> is found to be less than 1 %. Our simulation range extends upto temperature 1500 K, but the NIST data <cit.> is available only upto 700 K at 5000 bar. However, the systems are quite well behaved for these state points (T⩾ 700K) with stabilized energy, density and the velocity distributions follow the Maxwellian distribution, with the value of the standard deviation being consistent with the analytical value (√(k_BT/m)) over the entire temperature range.

Partially confined systems of LJ supercritical fluids are simulated in a cuboid with parallel walls facing each other normal to the z axis at z = ± H/2, H being the separation between the walls. We employ atomistic walls for introducing atomistic roughness to the  boundaries (Fig.<ref>).



For comparison we also simulate smooth, purely reflective walls at z = ± H/2. States with two different temperatures, 300K and 1500K at P = 5000 bar are chosen. These two (P,T) state points in confinement lie on either side of the Frenkel line: while the point A(5000 bar, 300K) lies in the "liquidlike" regime, the point B(5000 bar, 1500K) lies in the "gaslike" regime. For each of these two (P,T) state points we have simulated different partially confined systems with different spacing between the walls which are varied from 1 atomic diameter to about 21 atomic diameters, keeping the density same as that of the bulk phase for the corresponding P,T point. Periodic boundary conditions are applied along both the x and y axes. To investigate the changes in structural behaviour of confined supercritical fluid (Ar) close to the Frenkel line we also consider a range of P,T points spanning both sides of the Frenkel line.

The solid, atomistic walls are made of 3 layers of the face-centered cubic (fcc) lattice. The number of wall-atoms are varied from 1710076 (H = 3.4 Å) to 83668 (H = 70 Å) for the state point at 300K and from 3164188 (H = 3.4 Å) to 154588 (H = 70 Å) for the state point at 1500K. For modelling supercritical fluid confined between these two atomistic walls, 10^5 Argon particles are simulated using NVT ensemble. The wall atoms are attached to the lattice sites by harmonic springs. We set the spring constant (k) for these springs as 1000 ev/ Å^2 for 300K temperature to make the wall-atoms behave sufficiently rigid. We calculate the average mean squared displacement (MSD) of the wall-atoms and the average root-mean squared displacement (RMSD) of the wall atoms is found to be around 100 times smaller than a typical atomic displacement of fluid particles (∼ 1 Å) along confined direction (z). This confirms the sufficient rigidity of the walls. At 1500K temperature, keeping same spring constant (k = 1000 ev/Å^2) the RMSD of the wall atoms is found to be around 40 times smaller than a typical atomic displacement of fluid particles which assures moderate rigidity of the wall-atoms. To find out the implications of the rigidity of the walls on the structure, we also study the structural behaviour of supercritical Argon by varying the k-values for a specific confined spacing.

LJ potential has been used to model the interactions between both wall-fluid and wall-wall atoms. To model simple yet realistic walls, values of mass and the size of the wall-particles are taken from the calcium crystal data <cit.>, where mass of the each wall-atom is taken as 40.078 a.m.u. The LJ interaction parameters between the wall-atoms, ϵ_wall-wall = 0.2152ev  and σ_wall = 3.6 Å are used. We use two different fluid(Ar)-wall LJ interaction strengths (ϵ_wall-Ar) for our model: one same as the Ar-Ar interaction (0.0103ev) and the other a relatively stronger wall-Ar interaction (ϵ_wall-Ar = 0.0471ev), obtained using Lorentz-Berthelot (LB) mixing rule (ϵ_ij = √(ϵ_iϵ_j) ) <cit.>. The σ_wall-Ar is taken as 3.5 Å (LB-mixing rule: σ_ij = (σ_i + σ_j)/ 2). The cut-off distance for fluid-wall interaction has been taken as purely repulsive type (r_c:wall-fluid = σ_wall-Ar). The motion of wall-atoms is coupled to a thermostat of Nose-Hoover type to maintain same temperature as that of the supercritical fluid. This avoids any unnecessary heat flow though the fluid.  



§ RESULTS AND DISCUSSIONS



 §.§  Identifying Frenkel line from VACF and RDF of bulk supercritical fluid


The Velocity autocorrelation function (VACF) is generally defined as 

    VACF(t)=⟨v⃗(0)v⃗(t) ⟩

, where v⃗(0) and v⃗(t) are velocity vectors of particles at initial and some later time respectively and ⟨ ... ⟩ denotes the ensemble average. It is quite well understood that VACF is a monotonically decaying function for gases but it shows both oscillatory and decaying behaviour for liquids and solids <cit.>. M.E.Fisher and B.Widom introduced long back the idea of using oscillatory and monotonous decay of pair-correlation function 

as a marker to distinguish liquidlike and gaslike supercritical states using their analysis on linear continuum and lattice models <cit.>. Recently, V.V.Brazhkin et al. <cit.> showed that the minimum of VACF disappears when the supercritical fluid crosses the Frenkel line and goes from a "liquidlike" to a "gaslike" phase. We have chosen an isobaric line at 5000 bar in the P-T phase diagram <cit.> of supercritical Argon (Fig.<ref>) and examined systems over a range of temperatures from 240K to 1500K. We observe the expected gradual loss of oscillation of the bulk VACF as we increase the temperature. Increasing the temperature helps the particles to overcome the transient cage like environment created by nearest neighbours and diffuse in a "gaslike" manner. We observe that this change to occur around T ≈ 600-700K (Fig.<ref>).

Along with VACF, the structural properties of bulk SCF also undergo a well-defined change on crossing the Frenkel line as seen in the temperature evolution of the radial distribution function (RDF) at 5000 bar pressure and depicted in Fig.<ref>. We observe the gradual decrement of the height of the first peak and near disappearance of 2nd and 3rd peak of bulk RDF after crossing the Frenkel line. The pronounced local ordering (giving rise to the oscillatory feature in VACF) gradually dies out after crossing the Frenkel line and the fluid undergoes a transition from a "rigid liquidlike" to "gaslike" state (Fig.<ref>). It is also found that the rate of decay of the first peak height is faster than the 2nd peak height decay rate. It may be recalled that g(r) is related to the coordination number, the total number of particles found within a given sphere, in the following way 

    N(r) = ∫_0^rg(r) 4π r^2 dr


where N(r) is the number of particles found within the radius r. As the total number of particles is kept constant, the change in the local structure of the first coordination shell affects the structures in other neighbouring shells to satisfy the above equation. The integral implies that not all peaks can decay at the same rate with increasing temperature. While the first peak decays rapidly, the second peak decays at a relatively slower rate. We calculate the peak-heights and the first derivatives of the first and second peaks of RDF and vary them as a function of temperature. 


Two regimes can be noted from Fig.<ref>, which shows the derivative of the first peak-height as a function of temperature: one with a fast change (liquidlike regime) and the other with a slow change (gaslike regime)

of RDF peaks as a function of temperature. A structural crossover may be identified, by extrapolating the rates as shown in Fig.<ref>, to occur at around 600-700K. The change across the Frenkel line can be seen to be gradual. Identical feature has been observed for the second peak-height of the RDF (see the appendix, Fig.<ref>). Similar observations have also been made by D.Bolmatov et al. <cit.>. 

Two-body excess entropy is an alternate way of looking at how the degree of ordering changes with temperature. The presence of two distinct decay rates in RDF is manifested in the behaviour of the excess entropy as a function of temperature. This order parameter is defined <cit.> as

    s^(2) = -2 πρ k_B ∫[g(r)lng(r)-g(r)+1]r^2dr

where ρ is the number density and g(r) is the radial distribution function. It is well established that this two-body excess entropy (s^(2)) contributes approximately between 85 % and 95 % of the total excess entropy (s_ex)for a wide range of thermodynamic states for LJ fluids <cit.>. Fig.<ref> shows the variation of the negative two body excess entropy (-s^(2)/k_B) in a logarithmic scale as a function of temperatures for bulk supercritical Argon. The gradual decay of ln(-s^(2)/k_B) from a higher to lower values with increasing temperatures at a pressure of 5000 bar indicates the gradual loss of ordering as we go from a low temperature(high density) to a high temperature(low density) regime.

Fig.<ref> shows two distinct regimes characterized by distinct slopes (-0.003 for liquidlike regime and -0.0008 for the gaslike regime). The crossover of ordering is found to occur at around 600-700K temperature.



 §.§  Number Density Fluctuations and Compressibility in bulk supercritical fluid

We investigate another structural aspect, namely density fluctuations of the bulk supercritical LJ-fluid along the Frenkel line. It is well known that the number density fluctuations are closely related to the isothermal compressibility of the system through the relation

    κ_T = (V/k_BT)⟨ (Δρ_N)^2 ⟩/ρ_N^2
 
where, κ_T is the isothermal compressibility, ρ_N is number density (ρ_N = N/V), k_B is Boltzmann constant and ⟨...⟩ is the ensemble average <cit.>. The term ⟨ (Δρ_N)^2 ⟩/ρ_N^2 is nothing but the square of the standard deviation of the normalized number density (σ_ρ_N^2). We calculate the density fluctuations of supercritical Argon for each temperature ranging from 240K to 1500K in the form of distribution of densities. We find the σ_ρ_N^2 = ⟨ (Δρ_N)^2 ⟩/ρ_N^2 by calculating standard deviation of each of these distribution plots corresponding to each temperature. In Fig.<ref>.(a), we evaluate this temperature evolution of the density fluctuations for bulk supercritical Argon. We observe that at 5000 bar the widths of the normalized number density fluctuations are increasing with increasing temperature. The values of corresponding standard deviations of the density distributions as well as isothermal compressibility (κ_T) calculated from the above formula of bulk supercritical Argon are given in tabular form in the appendix (Table.<ref>).


Thermodynamically, isothermal compressibility is defined as <cit.>

    κ_T ≡ -1/V(∂ V/∂ P)_N,T


It is quite well known that as a small change in pressure leads to a larger volume change in gases due to it's low density and packing with respect to liquids, the gases have higher compressibility than liquids. We found that the isothermal compressibility(κ_T), evaluated from the standard deviation for the density fluctuations for bulk supercritical Argon, increases along the 5000 bar isobaric line with increasing temperature from 240K to 1500K. In Fig.<ref>.(b) we have plotted κ_T in logscale as a function of temperature, ranging from 240K to 1500K across the Frenkel line.








We find that the slope of the fitted straight line changes as we go from lower to higher temperatures (at fixed pressure 5000 bar). The crossover or changing of the slope of the fitted line happens around 600-700K temperature which can be interpreted as a Frenkel line transition of isothermal compressibility from a liquidlike to gaslike phase. Thus, the isothermal compressibility also undergoes a gradual crossover from a low compressible ("liquidlike") state to a highly compressible ("gaslike") state on crossing the Frenkel line.




 §.§ Confinement studies of supercritical Argon along the Frenkel line

To simulate partial confinement of supercritical LJ fluid, we employ atomistic LJ walls in a cuboid at z=± H/2, H being the spacing between the walls and periodic boundary conditions along x and y directions. We investigate structural features of the system by considering wall separations in the range 3.4 Å ⩽ H ⩽ 70 Å. 



  §.§.§ Distribution of LJ particles in supercritical phase normal to the walls :

As a first signature of the structural behaviour of the LJ fluid in supercritical phase under partial confinement, we investigate the distribution of the particles normal to the walls. Fig.<ref> shows some general features of the particle distribution normal to the walls for different wall-spacings at each of the two P,T state points A (P =5000 bar, T =300K) and B (P = 5000 bar, T =1500K). Thus A represents "liquidlike" phase before crossing the Frenkel line and B represents "gaslike" phase after crossing the Frenkel line. Fig.<ref>.(a) shows particle distributions in phase "A" exhibiting distinct layering normal to the walls. For all the spacings considered, there is a depleted region close to the walls. For larger wall spacings the layering is seen to be prominent only near the walls, with the featureless central region describing the average bulk density. In sharp contrast, layering is absent in phase "B" as can be seen from fig.<ref>.(b). However, depletion region can be seen for all the spacings. In both phases ("A" and "B") the depletion region arises from repulsive interactions between the fluid and wall particles. The extent of the depletion region in all cases is of the order of 1/2 atomic diameter at both ends. As the total number of particles is conserved, the particle number distribution exhibits (i) strong oscillatory features in the liquidlike high-density phase "A", (ii) weak oscillatory features in the gaslike phase "B".










  §.§.§ The Structural features normal to the walls : Before crossing the Frenkel line

It is intriguing to note that in phase "A", the distance between two successive peaks in the number distribution profile is found to be, on average, lower than one atomic diameter of the particle. It implies that the particles are not arranged in  monolayers along the width of the spacings. It is believed to arise from the packing frustration (a competition between the preferred packing of the fluid particles and the packing constraints imposed by the confining walls <cit.>) among the particles in supercritical phase. The effect of packing frustration may be seen more prominently for relatively smaller spacings. 







Fig.<ref> shows the distribution profile (Fig.<ref>.(a)) and the two-dimensional projections of particle configurations (Fig.<ref>.(b)) for a specific range of confined spacings (20 Å ⩽ H ⩽ 30 Å). Spacings with the ratio of H/σ close to integer values (H = 20 Å, 24 Å and 30 Å) appear to accommodate particles in layers characterized by large amplitude oscillations with the number of peaks scaling linearly with H/σ. The corresponding two-dimensional projections of particle distribution show well-formed periodic structures. When the spacings are not close to the integer values of H/σ (H = 22 Å, 26 Å and 28 Å), the layers, particularly in the central region of the width, exhibit weak oscillations, due to the frustration involved in accommodating the particles. The corresponding two-dimensional projection of particle configuration clearly shows the breakdown of ordered patterns around the central region (near z = 0). It is of interest to note that well-developed layers seem to form when H/σ is precisely an integer for hard spheres under confinement <cit.> for certain packing fractions. However, in our study we observe that not all integer multiple of σ spacings between the walls favour distinct ordering of supercritical fluid particles. While H = 7σ (23.8 Å) shows well-developed layering structure, H = 8σ (27.2 Å) doesn't show prominent layering (Fig.<ref>). After H = 30 Å, the distribution is nearly featureless in the central region around z = 0, corresponding to the bulk average density (Fig.<ref>). Thus H = 30 Å appears to indicate a sharp transition in particle number distributions along z. 





  §.§.§ Translational order parameter and Two-body excess entropy studies : Quantification of ordering normal to the walls

Although the particle distribution normal to the walls before crossing the Frenkel line clearly shows the ordering of the particles perpendicular to the walls due to confinement, there are several other measures that help to quantify the structure. Two such measures are the translational order parameter (τ) and two-body excess entropy (s^(2)), both of which can be determined by the radial distribution function normal to the walls denoted by g_⊥(r). g_⊥(r) can be defined as 

    g_⊥(r) ≡   1/ρ^2V∑_i≠ jδ(r -r_ij)[θ(|x_i-x_j|+δ x/2)
        -θ(|x_i-x_j|-δ x/2)]


where V is the volume, ρ is mass density, r_ij is the distance normal to the walls between molecules i and j, x_i is the x coordinate of the molecule i, and δ(x) is the Dirac δ function. The Heaviside function θ(x) restricts the sum to a pair of particles located in the same slab of thickness δ x. Translational order parameter (τ), which can be used as a tool to probe local density modulations in a system, is defined as <cit.>

    τ = ∫_0^r_c |g_⊥(r)-1|dr

where, r_c is the numerical cut-off of the RDF, along normal plane to the wall. For a completely uncorrelated systems like ideal gas τ is zero as g_r ≡ 1 for such systems. As system becomes ordered and structured, its value becomes relatively large <cit.>. We calculate τ from the g_⊥(r) for both the confined and bulk systems before and after crossing the Frenkel line. Fig.<ref>.(a) shows the translational order parameter (τ) variation for different confined spacings along the Frenkel line of supercritical Argon. For comparison, the corresponding values for bulk Argon in supercritical phase have been computed.


Another estimate of the structure of the confined system is the excess entropy, defined as the difference between the entropy of the probed system and the ideal gas calculated at a same density, temperature combination. We use the two-body excess entropy s^(2), defined <cit.> as

    s^(2) = -2 πρ k_B ∫[g_⊥(r)lng_⊥(r)-g_⊥(r)+1]r^2dr

As g_⊥(r) → 1 for completely uncorrelated and disordered systems, s^(2) vanishes(ideal gas behaviour) and becomes large and negative as an ordered structure starts forming <cit.>. Fig.<ref>.(b) shows the two-body excess entropy as a function of confined spacings for the supercritical Argon along the Frenkel line. 

For small spacings, the values of τ and the scaled s^(2) can be seen to be respectively higher and lower than the corresponding bulk phase values. It is worth noting that τ for bulk phases before and after the Frenkel line are positive, indicating that the supercritical bulk phases are more ordered than the ideal gas phase. Similarly, the bulk values of scaled s^(2) show non-zero negative values , which indicates that some ordering is present even in bulk with respect to the ideal gas structure. The variation of the scaled s^(2) with confinement is arising purely from ordering since the density has been kept constant for all confinements at a chosen P,T point.









In Fig.<ref>.(a) and Fig.<ref>.(b) we show the variation of τ and s^(2) with different confined spacings. We observe oscillating behaviour of both τ and s^(2) under the atomistic boundary confinement for a range of spacings (H ⩽ 30 Å) at 300K. Spacings leading to well developed sharper layers (comparatively lesser packing frustration) gives rise to higher values of τ and lower values of s^(2), in contrast to the spacings where under-developed layers are formed (packing frustration dominates). Fig.<ref> pictorially depicts these features. 

The spacing (H = 30 Å) at 300K, mentioned earlier, shows distinct jumps in both τ and s^(2) indicating significant loss of ordering. For spacings beyond H = 30 Å, τ and s^(2) can be seen to tend to their respective bulk values. In phase "B"(after crossing the Frenkel line) τ and s^(2) indicate mild ordering only for very small spacings. 



  §.§.§ Structural features parallel to the walls : Before and after crossing the Frenkel line

At 300K we have seen that confinement leads to layering over the width of the spacing. To understand the effect of the confinement on the structural behaviour of supercritical fluid parallel to the walls, we examine the parallel components of RDF for each of these layers. The radial distribution function for each of these layers can be evaluated from <cit.>

    g_∥(r) ≡   1/ρ^2V∑_i≠ jδ(r -r_ij)[θ(|z_i-z_j|+δ z/2)
        -θ(|z_i-z_j|-δ z/2)]

where V is the volume, ρ is the mass density, r_ij is the distance parallel to the walls between molecules i and j, z_i is the z coordinate of the molecule i, and δ(x) is the Dirac δ function. The Heaviside function θ(x) restricts the sum to a pair of particles located in the same slab of thickness δ z. We have considered δ z to be same as the width of each layer. We use a uniform bin width and bin number of 80 to calculate g_∥(r) for all the cases. Two well-defined classes of g_∥(r), henceforth labelled "P" and "Q", can be seen from Fig.<ref>.(a) and (b) respectively. g_∥(r) for Class "P" is characterized by the absence of positional shift of the peak positions for different layers. 









However, g_∥(r) for Class "Q" shows a gentle shift in the peak positions for different layers. With respect to the central layer the positional shift gradually increases and reaches a maximum for the layer closest to the wall. We note that class "Q" deals with confinement spacings which are not nearly equal to the integer multiple of H/σ leading to packing frustration. It has been found that g_∥(r) for H ⩾ 30 Å belongs to class "Q".

Confinement with H = 30 Å shows an additional feature although it belongs to class "P": it shows amorphous-like structure formation parallel to the walls (Fig.<ref>). Well-defined coordination spheres upto a length of 20 Å (≈ 6 σ) parallel to the walls, observed in g_∥(r), clearly shows more ordering than a typical "liquidlike"

phase, where local ordering persists only upto a maximum of 2 to 3 coordination spheres.

At 1500K (after crossing the Frenkel line), we don't see distinct structural features parallel to the walls. g_∥(r) exhibits only a single peak, shown in Fig.<ref> reminiscent of gaslike phase. It may be noted that g_∥(r) for the layer closest to the wall is no different from the g_∥(r) of the central region.


  §.§.§ Structural features of supercritical fluid for sufficiently narrow spacings under confinement: Before and After crossing the Frenkel line 

We investigate the structural features of supercritical Argon for extremely tight confinements (3.4 Å ⩽ H ⩽ 10 Å). 

Fig.<ref> shows the distinct layering of the particles in supercritical state at 300K, before crossing the Frenkel line. The number of layers formed due to confinement, scales linearly with the ratio of spacing to the atomic diameter. For example, a spacing of 1 atomic diameter(H=3.4 Å) gives 1 peak, 2 atomic diameters(H=6.8 Å) give 2 peaks and so on. Successive peaks are spaced lower than the atomic diameter of a particle, on average. Thus, despite having prominent layer arrangement of particles as seen from the number distribution profile, it is different from the close-packed structure.

To examine the role of the wall-fluid interaction potential, simulations are carried out with higher interaction strength(ϵ=0.0471ev). Fig.<ref> shows that, while, the general features of layering remain unaltered, there is a shift in the peak location and an enhancement of the peak heights.

Fig.<ref> shows g_∥(r) for each of these layers, exhibits

 structural features corresponding to amorphous phase at 300K. This property has also been seen in the case of g_∥(r) for H = 30 Å. Fig.<ref>.(b) and (c) show almost identical trends in g_∥(r) for different layers, corresponding to spacings of 6.8 Å and 10 Å respectively. 

Increasing the wall-fluid interaction strength leads to enhanced peak-heights in g_∥(r) as can be seen from Fig.<ref>.In Fig.<ref>.(a) a prominent central peak is seen for 3.4 Å spacing at 1500K, after crossing the Frenkel line. For all other spacings at 1500K, we observe weak layering near the wall and almost flat distribution in between. g_∥(r) for all spacings, including H = 3.4 Å, show structural features resembling the gas phase (Fig.<ref>.(b)).




  §.§.§ Evolution of structural features along the isobaric line under strong confinement


Until now we have discussed about the structural features under confinement for two state points of supercritical Argon at a pressure of 5000 bar: One(point A) before crossing the Frenkel line(300K) and the other(point B) after crossing the Frenkel line(1500K). Though these two state points A and B, clearly showed distinct liquidlike and gaslike features respectively, it will be interesting to see how the structural features evolve in the vicinity of the Frenkel line, along the isobaric path at 5000bar. We investigate the structural features of supercritical Argon for various state points close to the Frenkel line under the confinement with atomistic boundaries. We take the case of one narrow confined spacing(H = 10 Å) and closely monitor the transformation of the features of number distribution profiles along z (Fig.<ref>.(a)). 









In the liquidlike phase at 300K, three prominent layers are observed, resembling a highly structured fluid. This structural feature gradually decays as we go towards the gaslike phase crossing the Frenkel line. We observe a systematic decay of the peak-height of the central layer from a well defined value to complete disappearance as we go from 300K to 1500K along the Frenkel line (Fig.<ref>.(a)). This gradual disappearance of the central peak along the Frenkel line(∼ 600K-700K) clearly reconfirms the two-phase heterogeneity of supercritical fluids in confinement. The layers close to the walls, though decreased in heights, seem to exist even after crossing Frenkel line due to strong wall-particle correlation near the walls. g_∥(r) of the central layer also shows gradual transition with multiple coordination spheres(amorphous-like) transform to fewer (2,3) coordination spheres (liquidlike) and ultimately reduces to a single coordination sphere (gaslike) across the Frenkel line (Fig.<ref>.(b))









  §.§.§ Variation of structural features of supercritical fluid with varying rigidity of the atomistic walls:

In this last section we describe our investigation to analyse the role of the rigidity of the walls on the structural properties of supercritical fluid under confinement. We consider the confined spacing of 10 Å of supercritical Argon at a state point before the Frenkel line(300K) for this study. In our model, the wall-atoms are attached to their corresponding lattice sites by harmonic springs with the stiffness coefficient k. We vary the stiffness coefficient(k) of the springs attached to the wall-atoms from a very high value of 5000 ev/Å^2 (rigid, restricting the mean squared displacement (MSD) of the wall-atoms with respect to their lattice sites) to a very low value 0.005 ev/Å^2 (soft, increasing the mean squared displacement of the wall-atoms with respect to their lattice sites). 

The number distribution profiles for different k values are shown in Fig.<ref>.(a) for supercritical Argon at 300K with a confined spacing of 10 Å. For 10 ev/ Å^2 ⩽ k ⩽ 5000 ev/ Å^2, the peak-heights are found to be almost same. Further lowering of the k-value give rise to shorter peaks and for k ⩽ 0.05 ev/ Å^2 the peak heights decrease considerably to give flattened distribution of particles. Further, for these k values the depletion region vanishes.

While the ordering is quite similar for k > 0.5 ev/ Å^2, we observe a nearly 50 % decrement of the negative pair-excess entropy (-s^(2)) for k ⩽ 0.05 ev/ Å^2, which explains the flattening of peaks of the number distribution of supercritical fluid due to reduced ordering (see appendix, Table.<ref>). Decreasing the rigidity of the walls by  increasing the MSD of the wall atoms with respect to their lattice sites (lowering k value) shows a transition from an amorphous-like structure to a liquidlike ordering, parallel to the walls. Fig.<ref>.(b) shows the variation of g_∥(r) for the central layer of H = 10 Å spacing at 300K. Reduction of density of the layers causes this redistribution of the particles parallel to the walls from a highly ordered amorphous phase to a comparatively less ordered liquidlike phase.     



§ SUMMARY AND CONCLUSION

MD simulations have been carried out on bulk and partially confined supercritical LJ fluid to investigate the structural aspects of supercritical fluids across the Frenkel line. The study, done using LAMMPS, considered a system of 10^5  particles, interacting via the Lennard-Jones potential, at P = 5000 bar and temperatures ranging from 240K to 1500K simulating a wide range of densities of Argon.
 

VACF and RDF, evaluated using MD simulations for an isobaric line at 5000 bar over a temperature ranging from 240K to 1500K, confirm the characteristics of liquidlike and gaslike phases across the Frenkel line of supercritical Argon in the bulk. In this process the Frenkel line crossover point was identified to be in the range T ≈ 600-700K. Investigations of the density fluctuations in the bulk reveal that the changes in compressibility are consistence with the liquidlike phase going over to the gaslike phase as the Frenkel line is crossed. 












For the first time, MD simulations of confined supercritical fluid are reported for two P-T state points, one before (P= 5000 bar, T=300K) and one after crossing the Frenkel line of supercritical Argon (P= 5000bar, T=1500K) using both smooth and atomistic walls. At each P-T state point, the confinement spacing ranged from very narrow spacings like 3.4 Å (= 1 atomic diameter) to larger spacing of 70 Å (≈ 21 atomic diameters) while maintaining a constant density corresponding to the chosen P-T state. We further investigate the effect of confinement in the vicinity of the Frenkel line by considering state points across the line. The effect of rigidity of the walls on the structural properties of supercritical fluids is also studied in the context of confinement. 
 
In the "liquidlike" regime (point A, before crossing the Frenkel line), layering of particles perpendicular to the confining walls is pronounced. Analysis of the component of RDF that is parallel to the confining walls of successive layers reveals that the particles arrange themselves in a close-packing formation when smooth walls are imposed as boundaries(see appendix, Fig.<ref>). On the contrary, in the presence of atomistic walls this close-packing formation breaks down due to the appearance of near-wall depletion layers and associated packing frustration in the confinement. Also more ordered patterns are observed under atomistic wall-confined systems. The accommodation of the particles under confinement is governed by spacing between the walls, which, depending on less or more frustration in packing, selectively allows particles to form well-developed layer along z. We have found a spacing(H = 30 Å) at 300K (before crossing Frenkel line), where the distribution of particles exhibit maximum ordering and amorphous-like structure forms parallel to the walls. Immediately after which ordering disappears around z = 0. Rigidity of atomistic walls play a crucial role on the structural properties of the supercritical fluids. We observe a significant loss of ordering of the fluid, both normal and parallel to the walls, on modelling the walls softer by increasing the average MSD of the wall-atoms with respect to its lattice sites. This correlation between the ordering of fluids and the rigidity of the walls, has been confirmed by the two-body excess entropy measurements.
 


Extreme confinements with very narrow widths (spacing ∼ 1, 2 or 3 times the diameter of Argon) show prominent layering which can be enhanced using higher interaction strength between wall and fluid. More interestingly, amorphous-like structural features are confirmed by g_∥(r) along these well-defined layers parallel to the walls. Further, studying state points at the vicinity of the Frenkel line under confinement shows the gradual disappearance of layering across the Frenkel line. In the “gaslike” regime (point B, after crossing the Frenkel line), such ordering is not seen apart from the region very close to the walls.


The structural ordering has been quantified by two-body excess entropy and translational order parameter calculations which clearly indicate the correlation between confinement and ordering in a supercritical fluid. 



The discrete jump in translational order parameter at the spacing (H = 30 Å) confirms the sudden loss of ordering for H > 30 Å (Fig.<ref>.(a)).

While layering is known to occur in normal fluids at high density, it is not typical of liquids to show amorphous-like structures parallel to the walls under confinement at 300K. Close-packing like structural feature, found in supercritical fluids in the liquidlike regime under smooth walls, is also not a  general feature for liquids. Since layering is not significant in the "gaslike" phase of the supercritical fluid, it is unlikely to occur for normal fluids in the "gaslike" regimes as they have much lower densities.
 
The structural aspects of very narrow confined systems further open up possibilities of unusual dynamics parallel to the walls, where "liquidlike" phase of SCF suffers a 





transition to amorphous-like fluid. Also, our findings on the wall-rigidity dependence of the structural features of SCF offer new scopes to elucidate the relation between the supercritical fluid and wall dynamics. The key features of our findings are summarized in a tabular form in Table <ref>.

The presence of ordering in the ‘liquidlike” regime and the absence of it in the “gaslike” regime under confinement may have significant implications for the variation of transport properties across the Frenkel line. This heterogeneity, present both in bulk and confined systems, might be responsible for breakdown of the universal scaling  between structure and dynamics of fluids and stimulates possibilities of having a unique relationship between them. 


§ ACKNOWLEDGEMENT

We acknowledge the help of HPCE,IIT Madras for high performance computing. KG would like to express his gratitude to Department of Science and Technology(DST)(India) for INSPIRE Fellowship for funding.





§ BULK STUDIES OF SUPERCRITICAL FLUID



 §.§ Temperature variation of First and second peak-height of RDF(radial distribution function) and its derivatives in bulk phase of supercritical Argon for all P,T state points chosen for study at P = 5000 bar: Identification of Frenkel line:

Figure <ref> describes the temperature variation of the first
and second RDF peaks in detail.
















§ DATASET FOR NUMBER DENSITY FLUCTUATIONS AND COMPRESSIBILITY IN BULK SUPERCRITICAL ARGON

Table. <ref> shows the dataset of compressibility in bulk
supercritical argon at 5000 bar across the Frenkel line.


§ CONFINEMENT STUDIES OF SUPERCRITICAL FLUID



 §.§  Layering of supercritical fluid normal to the walls before crossing the Frenkel line(300K): Comparison between smooth, purely reflective and atomistic walls

We observe at wider spacings (e.g. H = 70 Å) the number distribution profiles are quite similar with the presence of bulk number distribution around z= 0 for both flat and atomistic walls. Prominent differences can be seen as we go towards smaller spacings (H = 20 Å and 30 Å) (Fig.<ref>).



 §.§  Translational order parameter and Two-body excess entropy studies normal to the purely reflective walls:

The smooth, purely reflective boundaries manifest nearly monotonous trends of τ and s^(2) as a function of confined separations. τ and s^(2) show monotonically decreasing and increasing trends respectively towards the bulk value with increasing spacing (Fig.<ref>) . 



 §.§  Radial distribution function parallel to the smooth, purely reflective walls before crossing Frenkel line:


Smooth, purely reflective boundaries, unlike to the atomistic walls, show only one class of g_∥(r) where positional shift of the coordination spheres is found for all spacings ranging from 20 Å to 70 Å (class Q) (Fig.<ref>). These positional shifts, along with the periodic the mono-layer formation along z, suggest a close packing structure parallel to the flat walls and are depicted in Fig.<ref>. For small spacings the close packing structure is seen to prevail across the entire thickness (Fig.<ref>). 

As the spacing(H) between the walls increases, the close packing arrangement gradually seems to disappear in the central region of the confinement.


 §.§  Radial distribution function parallel to the smooth, purely reflective walls after crossing Frenkel line:

After crossing the Frenkel line, layering becomes insignificant. Fig.<ref> shows one such case at 1500K temperature for 20 Å wall-spacing for smooth walls, where we don't see much change in comparing the g_∥(r) of central region and that of the region close to the walls. 




 §.§  Excess entropy dataset for different k values (stiffness coefficient)of the walls for 10 Å spacing at 300K:

Table <ref> describes the variation of pair-excess entropy as a function of stiffness co-efficients of the walls.































39
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[E.Kiran et al.(2000)E.Kiran, P.G.Debenedetti, and C.J.Peters]Kiran2000
  author author E.Kiran, author P.G.Debenedetti,  and author
  C.J.Peters, title Supercritical
  fluids: Fundamentals and applications,  (publisher Kluwer
  Academic Publishers, address Dordrecht, year
  2000)NoStop
[J.McHardy and S.P.Sawan(1998)]McHardy1998
  author author J.McHardy and  author S.P.Sawan, title
  Supercritical fluid cleaning: Fundamentals, technology and applications,   (publisher Noyes Publications, address New
  Jersey, year 1998)NoStop
[J.Wang et al.(2009)J.Wang,
  C.N.Sheaff, B.Yoon, R.Addleman, and C.Wai]Wang2009
  author author J.Wang, author C.N.Sheaff, author B.Yoon, author R.Addleman,  and author C.Wai, @noop  journal
  journal Chemistry volume 15,  pages 4458 (year 2009)NoStop
[M.E.Fisher and B.Widom(1969)]Fisher1969
  author author M.E.Fisher and  author B.Widom, @noop  journal journal J.Chem.Phys volume
  50, pages 3756 (year 1969)NoStop
[K.Nishikawa and I.Tanaka(1995)]Nishikawa1995
  author author K.Nishikawa and  author I.Tanaka, @noop  journal journal Chem.Phys.Lett volume
  244, pages 149 (year 1995)NoStop
[T.Morita et al.(1997)T.Morita, K.Nishikawa, M.Takematsu,
  H.Iida, and S.Furutaka]Morita1997
  author author T.Morita, author K.Nishikawa, author M.Takematsu, author H.Iida,  and author S.Furutaka, @noop  journal
  journal J. Phys. Chem. B volume 101,  pages 7158 (year 1997)NoStop
[F.Gorelli et al.(2006)F.Gorelli, M.Santoro, T.Scopigno,
  M.Krisch, and G.Ruocco]Gorelli2006
  author author F.Gorelli, author M.Santoro, author T.Scopigno, author M.Krisch,  and author G.Ruocco, @noop  journal
  journal Phys. Rev. Lett volume 97,  pages 245702 (year 2006)NoStop
[G.G.Simeoni et al.(2010)G.G.Simeoni, T.Bryk, F.A.Gorelli,
  M.Krisch, G.Ruocco, M.Santoro, and T.Scopigno]Simeoni2010
  author author G.G.Simeoni, author T.Bryk, author F.A.Gorelli, author M.Krisch, author G.Ruocco, author M.Santoro,  and author T.Scopigno,  @noop  journal journal Nature Physics  volume 6, pages 503–507 (year 2010)NoStop
[D.T.Banuti et al.(2017)D.T.Banuti, M.Raju, and M.Ihme]Banuti2017
  author author D.T.Banuti, author M.Raju,  and author M.Ihme, @noop  journal journal Phys.  Rev. E volume 95, pages 052120
  (year 2017)NoStop
[C.Yang et al.(2015)C.Yang,
  V.V.Brazhkin, M.T.Dove, and K.Trachenko]Yang2015
  author author C.Yang, author V.V.Brazhkin, author M.T.Dove,  and author K.Trachenko,  @noop  journal journal Phys. Rev. E  volume 91, pages 012112 (year 2015)NoStop
[V.V.Brazhkin et al.(2012)V.V.Brazhkin, Yu.D.Fomin, A.G.Lyapin,
  V.N.Ryzhov, and K.Trachenko]Brazhkin12
  author author V.V.Brazhkin,
  author Yu.D.Fomin, author
  A.G.Lyapin, author V.N.Ryzhov,   and author K.Trachenko, @noop 
  journal journal Phys. Rev. E volume 85, pages 031203 (year
  2012)NoStop
[V.V.Brazhkin et al.(2013)V.V.Brazhkin, Yu.D.Fomin, A.G.Lyapin,
  V.N.Ryzhov, E.N.Tsiok, and K.Trachenko]Brazhkin2013
  author author V.V.Brazhkin,
  author Yu.D.Fomin, author
  A.G.Lyapin, author V.N.Ryzhov,
  author E.N.Tsiok,  and author
  K.Trachenko, @noop  journal journal Phys. Rev. Lett volume 111, pages 145901 (year 2013)NoStop
[D.Bolmatov et al.(2013)D.Bolmatov, V.V.Brazhkin, and K.Trachenko]Bolmatovnew2013
  author author D.Bolmatov, author V.V.Brazhkin,  and author
  K.Trachenko, @noop  journal journal Nat. Commun. volume 4, pages 2331 (year 2013)NoStop
[Bolmatov et al.(2013)Bolmatov, V.V.Brazhkin, Yu.D.Fomin,
  V.N.Ryzhov, and K.Trachenko]Dima2013
  author author D. Bolmatov, author V.V.Brazhkin, author Yu.D.Fomin, author V.N.Ryzhov,  and author K.Trachenko,  @noop  journal journal J.Chem.Phys  volume 139, pages 234501 (year 2013)NoStop
[D.Bolmatov et al.(2015)D.Bolmatov, M.Zhernenkov, D.Zav'yalov,
  S.N.Tkachev, A.Cunsolo, and Y.Q.Cai]Bolmatov2015
  author author D.Bolmatov, author M.Zhernenkov, author D.Zav'yalov, author S.N.Tkachev, author A.Cunsolo,  and author Y.Q.Cai, @noop  journal journal Sci.
  Rep. volume 5, pages 15850 (year 2015)NoStop
[Y.D.Fomin et al.(2015)Y.D.Fomin, V.N.Ryzhov, E.N.Tsiok, and  V.V.Brazhkin]Fomin2015
  author author Y.D.Fomin, author V.N.Ryzhov, author E.N.Tsiok,  and author V.V.Brazhkin,  @noop  journal journal Sci. Rep.  volume 5, pages 14234 (year
  2015)NoStop
[Y.D.Fomin et al.(2014)Y.D.Fomin, V.N.Ryzhov, E.N.Tsiok,
  V.V.Brazhkin, and K.Trachenko]Fomin2014
  author author Y.D.Fomin, author V.N.Ryzhov, author E.N.Tsiok, author V.V.Brazhkin,  and  author K.Trachenko, @noop  journal journal Sci. Rep. volume
  4, pages 7194 (year 2014)NoStop
[D.Bolmatov et al.(2014)D.Bolmatov, D.Zav’yalov, Gao, and  Zhernenkov]Bol2014
  author author D.Bolmatov, author D.Zav’yalov, author M. Gao,  and author M. Zhernenkov, @noop  journal
  journal J. Phys. Chem. Lett. volume
  5, pages 2785−2790 (year 2014)NoStop
[Yu.D.Fomin et al.(2015)Yu.D.Fomin, V.N.Ryzhov, E.N.Tsiok, and  V.V.Brazhkin]Fomin2
  author author Yu.D.Fomin, author V.N.Ryzhov, author E.N.Tsiok,  and author V.V.Brazhkin,  @noop  journal journal Phys. Rev. E  volume 91, pages 022111 (year 2015)NoStop
[M.Stubbs(2016)]Stubbs2016
  author author J. M.Stubbs, @noop  journal journal J.
  of Supercritical Fluids volume 108, pages 104 (year 2016)NoStop
[Mittal et al.(2006)Mittal,
  R.Errington, and M.Truskett]Mittal2006
  author author J. Mittal, author J. R.Errington,
   and author T. M.Truskett,  @noop  journal journal Phys. Rev. Lett  volume 96, pages 177804 (year 2006)NoStop
[Mittal et al.(2008)Mittal,
  M.Truskett, and G.Hummer]Mittal2008
  author author J. Mittal, author T. M.Truskett,
   and author G.Hummer, @noop 
  journal journal Phys. Rev. Lett volume 100, pages 145901 (year
  2008)NoStop
[de Sousa et al.(2016)de Sousa, J.J.Sáenz, Scheffold,
  A.García-Martín, and L.S.Froufe-Pérez]Sousa2016
  author author N. de Sousa, author J.J.Sáenz, author
  F. Scheffold, author
  A.García-Martín,  and author L.S.Froufe-Pérez, @noop  journal journal J.Phys.:Condens.Matter volume 28,  pages 135101 (year 2016)NoStop
[Lam(2004)]Lammps
  http://lammps.sandia.gov title LAMMPS Software
  package:http://lammps.sandia.gov (year 2004)NoStop
[NIS()]NIST
  http://webbook.nist.gov/chemistry/ title NIST
  Data Set:http://webbook.nist.gov/chemistry/NoStop
[Stukowski(2010)]OVITO
  author author A. Stukowski, @noop  journal journal
  Modelling Simul. Mater. Sci. Eng. volume 18,  pages 015012 (year 2010)NoStop
[A.J.Markvoort et al.(2005)A.J.Markvoort, P.A.Hilbers, and S.V.Nedea]Markvoort2005
  author author A.J.Markvoort,
  author P.A.Hilbers,  and author
  S.V.Nedea, @noop  journal journal Phys. Rev. E volume 71, pages 066702 (year 2005)NoStop
[M.P.Allen and D.J.Tildesley(1987)]Allen1987
  author author M.P.Allen and  author D.J.Tildesley, title Computer simulation of liquids,  (publisher
  Oxford,Clarendon, address New York, year
  1987)NoStop
[J.P.Hansen and I.R.McDonald(2007)]Hansen2007
  author author J.P.Hansen and  author I.R.McDonald, title
  Theory of simple liquids,  (publisher Elsevier, address Amsterdam, year 2007)NoStop
[R.E.Nettleton and M.S.Green(1958)]Nettleton1958
  author author R.E.Nettleton and  author M.S.Green, @noop  journal journal J.Chem.Phys volume
  29, pages 1365 (year 1958)NoStop
[R.D.Mountain and H.J.Raveche(1971)]Mountain1971
  author author R.D.Mountain and  author H.J.Raveche, @noop  journal journal J.Chem.Phys volume
  55, pages 2250 (year 1971)NoStop
[Baranyai and J.Evans(1989)]Baranyai1989
  author author A. Baranyai and author D. J.Evans, @noop  journal journal
  Phys. Rev. A volume 40, pages 3817
  (year 1989)NoStop
[McQuarrie(1976)]McQuarrie1976
  author author D. McQuarrie, title Statistical mechanics,   (publisher Harper & Row, address New York,  year 1976) Chap. chapter 3NoStop
[Nygård et al.(2014)Nygård, Sarman, and Kjellander]Kim2014
  author author K. Nygård, author S. Sarman,   and author R. Kjellander,  @noop  journal journal J.Chem.Phys  volume 141, pages 094501 (year 2014)NoStop
[Mittal et al.(2007)Mittal,
  R.Errington, and M.Truskett]Mittal2007
  author author J. Mittal, author J. R.Errington,
   and author T. M.Truskett,  @noop  journal journal J.Chem.Phys  volume 126, pages 244708 (year 2007)NoStop
[Pant and K.Ghorai(2016)]Pant2016
  author author S. Pant and author P. K.Ghorai, @noop  journal journal Mol.
  Phys. volume 114, pages 1771
  (year 2016)NoStop
[T.M.Truskett et al.(2000)T.M.Truskett, S.Torquato, and P.G.Debenedetti]Truskett2000
  author author T.M.Truskett,
  author S.Torquato,  and author
  P.G.Debenedetti, @noop  journal
  journal Phys. Rev. E volume 62,  pages 993 (year 2000)NoStop
[Kumar et al.(2005)Kumar,
  V.Buldyrev, W.Starr, Giovambattista, and Stanley]Pradeep2005
  author author P. Kumar, author S. V.Buldyrev,
  author F. W.Starr, author N. Giovambattista,  and author H. Stanley, @noop 
  journal journal Phys. Rev. E volume 72, pages 051503 (year
  2005)NoStop
[Krott and Bordin(2013)]Krott2013
  author author L. B. Krott and author J. R. Bordin, @noop  journal journal
  J.Chem.Phys volume 139, pages 154502
  (year 2013)NoStop

