Measuring pressure in equilibrium and nonequilibrium lattice-gas models

Pressure is a key thermodynamic parameter encoding the phase behavior of macroscopic systems in thermal equilibrium.¹ As first realized by Tammann long ago,^2,3 for systems with many degrees of freedom, the pressure behavior may not be as a simple as a mechanical interpretation would suggest. This is essentially due to the entropic contribution coming from the integral (or the sum) over the whole configuration space, which is at the root of several unusual and counter-intuitive phenomena, such as inverse melting.^4–6 There is also a fundamental interest in extending the notion of pressure to active matter, a crucial step toward the formulation of a nonequilibrium statistical mechanics.⁷ It is therefore of paramount importance to rely on flexible and viable algorithms that are capable of measuring the pressure directly, with no reference to a special state as in the standard thermodynamic integration method. This is all the more true for (athermal) lattice systems where Monte Carlo algorithms are far less developed than their molecular dynamics counterpart.^8,9 Previously proposed methods generally sacrifice the usual periodic boundary as they involve hard walls and volume fluctuations due to the removal or addition of a layer next to the wall.¹⁰ We consider here an elegant proposed some time ago by Dickman in the context of lattice polymers,^11,12 which bears a strong resemblance with the ghost particle insertion trick suggested by Widom for computing the chemical potential.¹³ The Dickman method consists of three simple yet important observations. First, for a lattice gas with volume V and Helmholtz free energy F, the pressure P = −∂F/∂V is written as

where Z is the canonical partition function. Then, one observes that the partition function of a system with volume V is equivalent to that of a modified system with volume V + ΔV in which one has added an infinite repulsive potential acting on the sites contained in ΔV (preventing particles from entering ΔV). Thus, if we introduce a repulsive potential of finite magnitude U_λ = −k_BT ln λ acting on the sites in ΔV (in such a way that λ = 0 corresponds to U = ∞, and λ = 1 corresponds to U = 0) and denote with Z_λ the partition function of this modified system, we get Z(V) ≡ Z₀(V + ΔV) and Z(V + ΔV) ≡ Z₁(V + ΔV), whence the partition function ratio in Eq. (1) becomes

To transform this into something useful for the design of a Monte Carlo algorithm, one finally needs to remark that the partition function Z_λ for the modified system can be written as

where the summation is taken over the possible configurations of ΔN particles in ΔV and Ω is the total number of configurations with ΔN particles in ΔV and N + ΔN in V + ΔV. Thus, the average number ⟨ΔN⟩ of particles in ΔV can be obtained as

This leads to the desired expression for the pressure

The practical implementation of Eq. (5) in a Monte Carlo algorithm obviously requires the interval of λ values be discretized and numerical simulations carried out for different finite values of the repulsive potential. In addition, the sites of ΔV on which the potential acts needs to be finite. To the best of our knowledge, this is typically chosen to be as small as 1, and the single pressure-probe site is generally embedded in a wall of codimension 1 confining the system from one side. In this form, it has been observed that the algorithm becomes problematic in the presence of a strong external bias, as in a lattice gas under gravity.¹⁴ When dealing with pressures as tiny as those generated by the Casimir effect in a nonequilibrium driven system,¹⁵ the algorithm may require excessively long execution times, which is all the more troublesome for cooperative driven dynamics. Moreover, a naive extension of the algorithm allowing probing the pressure at different locations, simultaneously, leads to unreliable results in nonequilibrium conditions.

We show in this paper that these drawbacks can be actually overcome and that an efficient formulation of the algorithm can be consistently applied to a variety of physical situations of growing complexity, from a simple sedimentation problem to fluctuation-induced forces in a driven diffusive system. First of all, one should notice that some care is required when dealing with systems in which ⟨ΔN⟩ ≈ ΔV. In these cases, there exists a limitation on the maximum value of pressure P_max that can be measured, as Eq. (4) diverges for λ → 0, while P_max ≈ −log λ_min, (where λ_min is the minimum value of the λ discretization). Evidently, to correctly quantify the pressure when ⟨ΔN⟩ ≈ ΔV, one has to choose a suitable value of λ_min. Second, it is important to emphasize that the presence of a confining wall is not a necessary ingredient for computing pressure in a Monte Carlo simulation. If for some reason the spatial translational invariance of system geometry and interactions has to be preserved, one can use periodic boundary conditions (unless, of course, one is interested in the wall effects, as in wetting phenomena, for example). Third, the size of ΔV can be relatively large, ΔV ≫ 1, provided that it is sub-extensive, ΔV ≪ V. This improves data statistics and becomes particularly relevant when the system under study displays an inhomogeneous density profile. In this case, multiple pressure probes distributed uniformly over the whole system can give access to the entire pressure profile with a single simulation run. Finally, once the wall is suppressed, it is important to specify how the (ghost) sites of ΔV are added to the original system. We do so by adding an extra nearest neighboring site to a certain number, ΔV, of sites of the original system. In this way, the local geometric structure of the original system is preserved, and its intrinsic dynamics will be only weakly perturbed by the presence of these extra ghost sites (for a schematic representation, see Fig. 1). This point is particularly important as statistical mechanical models of interest generally involve static or dynamic spatially extended interactions.

In the following, the above points are implemented and carefully examined in systems for which the pressure has a purely entropic content, i.e., athermal systems. Comparison with theory is made whenever analytical results are partially or fully available. Monte Carlo simulations are carried out in the grand canonical ensemble, in which the system of volume V (or one of its parts) is coupled to a reservoir R at chemical potential μ = log[ρ_R/(1 − ρ_R)], where ρ_R is the reservoir density. Unless otherwise mentioned, the pressure-probe ghost sites are equally spaced in both the transverse and longitudinal directions of the system. At every Monte Carlo time step, a lattice bond in the coupled system (ΔV + V + R) is randomly selected, and depending on the pair of sites involved, one of the three possibilities is considered: hole-particle exchange within the volume V, hole-particle exchange between V and the reservoir R, and hole-particle exchange between V and the ghost volume ΔV. It is important to emphasize that the dynamics of the particle reservoir R is not explicitly simulated. Rather, every time the randomly selected bond involves a site of R, its occupation state is randomly established by extracting a number with probability ρ_R. The results we present were obtained by using the cutoff λ_min = 0.001 and discretizing λ ∈ [0, 1] into 24 values equally spaced over three contiguous subintervals [0.01, 0.04], [0.04, 0.2], and [0.2, 1], with steps Δλ = 0.01, 0.02, and 0.1, respectively. Compared to a simple uniform discretization of [0, 1], this makes the numerical evaluation of Eq. (5) more effective because the contribution of smaller values of λ, which give a large contribution to the integral, is better weighted.

To begin with, we discuss what is perhaps the simplest equilibrium system with a non-trivial inhomogeneous density profile, the hard-core lattice gas under gravity. This provides a natural test bed for the algorithm as the system thermodynamics is simple enough to allow for a detailed comparison with the exactly known results. For a system of height h_top enclosed on the bottom and with the top in contact with a reservoir at density ρ_top, the pressure at height h is given by

where ρ(h) is the density profile defined as

with mg being the gravity force. A further advantage of this problem is that the average number ⟨ΔN⟩ of particles in ΔV, Eq. (4), can be computed easily,¹⁴ 

providing a further detailed check of the algorithm. In the Monte Carlo, simulation we assume that the boundary condition in the direction transverse to the gravity force is periodic and the ghost sites are uniformly distributed along the system diagonal. The square lattice is tilted by 45°, and particle hopping is ruled by the standard Metropolis transition probability, p = min{1, exp(−βmgΔh)}, where Δh = ±1 is the height difference due to a single hop. In Fig. 2, we show the results for the density profile P(h) for two different values of the gravity force, when the system top is kept in contact with a particle reservoir at density ρ_top = 0.01. As we can see, the analytic results are perfectly recovered.¹⁶ In this form, the algorithm can be applied to more complex lattice models of colloidal sedimentation and granular compaction, like those of Refs. 17 and 18.

Next, we turn to a more interesting system of interacting particles, that is, the lattice gas with nearest neighbor exclusion, where particles are forbidden to occupy the same or neighboring sites. This model has a long tradition in the study of critical phenomena as a schematic lattice model of hard-sphere fluids and is closely related to the classical NP-hard combinatorial optimization problem of finding a minimum vertex cover on a random graph.¹⁹ Numerical and theoretical studies for various finite dimensional lattices^20–27 have shown that there exists a critical density ρ_c (for a square lattice ρ_c = 0.367 74 …) at which the model exhibits a continuous phase transition to an ordered state with exponents that belong to the Ising universality class.²⁵ Above ρ_c, particles tend to occupy preferentially one of the two sublattices in which the original lattice can be decomposed (in such a way that the nearest-neighbors of any sites of one sublattice belong to the other sublattice and vice-versa). The ground state, therefore, corresponds to a fully occupied sublattice, and the maximum packing density is ρ_max = 1/2. In Fig. 3, we show the parametric plot of the equation of state P(ρ) obtained in a grand-canonical simulation by slowly increasing the chemical potential. In agreement with the well-established numerical and theoretical results, the equation of state displays an inflection point, corresponding to a divergent compressibility, near the critical density at which sublattice ordering occurs. The full line below the critical density is a one-parameter fit obtained from the approach in Ref. 27, while the full line at higher densities is the non-interacting entropic pressure −log(1 − ρ/ρ_max) near the closest packing density ρ_max = 1/2.

The previously considered model can be interestingly generalized to situations in which every particle cannot have more than ℓ particles as nearest neighbors (thus, the nearest neighbor exclusion models is recovered for ℓ = 0). This class of models was introduced in the attempt of providing a microscopic finite-dimensional realization of the thermodynamic scenario for the elusive glass transition.²⁸ Subsequent works have shown that in finite dimension, a freezing transition toward an ordered ground state occurs.^30,31 However, some more sophisticated variants appear to be quite stable against crystallization and display properties of fragile glass-forming liquids.^32,33 For the purpose of testing the present algorithm in a geometric structure other than the usual Euclidean lattices, we focus here on the monodisperse lattice glass model on the Bethe lattice. This is also known as the regular random graph, which is a locally tree-like geometric structure with a fixed connectivity k + 1, and with long loops of order ln N. This latter feature, an analog of generic boundary conditions, prevents unphysical surface effects and gives rise to frustration precluding crystalline orderings. The resulting glassy behavior can be studied analytically with the tools of disordered systems theory^28,29 and has been extended to the quantum domain as well.³⁴ In Fig. 4, we show a parametric plot of the equation of state obtained in a grand canonical Monte Carlo simulation for the Bethe lattice glass (with connectivity k + 1 mimicking the square and cubic lattice) for two values of ℓ. Consistent with the analytical calculation,²⁸ no sign of crystallization is observed in the equation of state and, for large chemical potential, the pressure tends to diverge at a density very near the closest packing limit, μ → ∞, obtained with the cavity method of disordered systems theory.²⁹ 

As a final more stringent test of the present algorithm, we now consider a two-dimensional system driven into a nonequilibrium steady state by two reservoirs at different densities, ρ₀ and ρ₁, located at its edges x = 0 and x = L, respectively. For such driven diffusive systems, it is known that there exists long-range correlations,³⁵ even for purely hard-core interactions.³⁶ A striking consequence is the appearance of Casimir-like forces in a finite geometry.¹⁵ To the leading order in ρ₁ − ρ₀ and for a simple symmetric exclusion process³⁷ on a square lattice of size L × d, this fluctuation-induced pressure, Π(x), is given by¹⁵ 

where ρ(x) = ρ₀ + (ρ₁ − ρ₀)x/L is the system density profile, and we set k_BT = 1. For the lattice system with only infinite hard-core repulsion, the thermodynamic limit of the average density of particle, ⟨ρ_ΔV⟩, in the ghost volume ΔV is exactly known, see Eq. (8). This means that the Casimir pressure, Π_MC(x), in a Monte Carlo simulation with a finite ΔV can be measured as

where it is understood that the integral has to be properly discretized. In Fig. 5, we compare the Monte Carlo results for Π_MC(x) at reservoirs densities ρ₀ = 0.1 and ρ₁ = 0.7 and different system sizes, with the theoretical result for Π(x), Eq. (9), in the rescaled form $Π̃(x)=LdΠ(x)$⁠. Data points of each curve were simultaneously obtained for a dynamical evolution long 10⁹ Monte Carlo sweeps and averaged over a sample of 10–20 elements (statistical errors are not visible because of the order of symbol size). The agreement is excellent and fully confirms the efficiency and functionality of our algorithm.

In conclusion, we have developed an algorithm based on the Dickman method which can considerably improve the accuracy of Monte Carlo calculations of pressure in lattice models, and is especially suited for systems with inhomogeneous density profiles, both in equilibrium and nonequilibrium steady states. The algorithm is quite general and holds the key to many other statistical mechanics applications.

I gratefully acknowledge Daniel Martin and José Luis Iguain for valuable suggestions in the early stage of this work. I also thank Avi Aminov, Jef Arenzon, and Yariv Kafri for kindly clarifying some aspects of their work.^15,27 Numerical simulations were performed by using ICTP computer facilities.

The data that support the findings of this study are available from the corresponding author upon reasonable request.