Limits of the R-matrix analysis of Mayer cluster coefficients: the spectrum-shift regime, reality of the matrix elements, and error-controlled critical activities

For purely repulsive (hard-core) lattice and continuum gases, the Mayer cluster coefficients can be encoded in a real symmetric tridiagonal R matrix whose elements converge to two constants; the edges of its spectrum fix the two termination points of the fluid equation of state in the activity variable, and the behavior of the first eigenvector component at those edges yields exact critical exponents for the density singularity (Rotman and Eisenberg, Phys. Rev. E 80, 031126 (2009)), generalizing the earlier 1/n^2 treatment of Eisenberg and Baram. The same paper identified precisely where this exact analysis fails, and these failures remain open. First, when the fitted correction amplitudes are so strong that the radicand of the exponent formula becomes negative, the matrix spectrum edge shifts away from its constant-matrix value, and neither the termination activity nor the critical exponent in this regime is known (realized in the triangular-lattice gas with first- and second-neighbor exclusion, where the R-matrix fit gives a termination activity of 9.35 against the Monte Carlo value 5.794). Second, the reality of the R-matrix elements, on which the whole construction rests, is observed in every model computed to date but remains unproven for general purely repulsive models. Third, the extraction of the critical activity from finite-order fitted parameters is highly sensitive to small errors in the asymptotic constants, and no error-controlled methodology exists. Of the source's open items, only the square-lattice N4 (exclusion of all sites within the first four neighbor shells) universality-class dispute has since been settled (in favor of the Ising assignment, with the discovery of a second transition), leaving these three as the surviving open core: extend the R-matrix treatment to the spectrum-shift regime, prove or refute the reality of the construction, and quantify the uncertainty of the extracted critical activities.

Background

Hard-core lattice gases, in which each occupied site excludes a fixed set of neighboring sites so that every Boltzmann weight is 0 or 1 and any ordering transition is driven purely by entropy, are standard laboratories for entropy-driven phase transitions. Writing the grand-canonical thermodynamics in terms of the activity (fugacity) $z$, the fluid-phase density is given by the Mayer cluster series $\rho(z)=\sum_{n\ge 1}n b_n z^n$, where $b_n$ are the Mayer cluster integrals; for lattice models these can be obtained exactly to finite order by transfer-matrix enumeration (exact counting by matrix products over lattice rows), and for the exactly solved hard-hexagon model (nearest-neighbor exclusion on the triangular lattice; critical activity $z_c=(11+5\sqrt{5})/2$, density exponent $2/3$) the exact equation of state yields hundreds of coefficients (Baxter 1980; Joyce 1988). For purely repulsive interactions the radius of convergence of this series is set by a nonphysical singularity on the negative real $z$-axis, whose exponent is universal ($\sigma=1/6$ in two dimensions; Baram and Luban 1987 and later confirmations), and this nearby branch point defeats standard series-analysis tools such as ratio methods, Dlog-Pade (Padé approximation of the logarithmic derivative of the series), and differential-approximant methods, which then fail to detect the physical ordering transition at all.

Baram, Eisenberg and co-workers showed that the cluster integrals can be rearranged into a real symmetric tridiagonal matrix, the R matrix, defined by $(R^n)_{11}=(-1)^n(n+1)b_{n+1}$ for all $n\ge 1$. The density then takes the spectral form $\rho(z)=\sum_{\lambda}\psi_1(\lambda)^2/(z^{-1}+\lambda)$, where the sum runs over the eigenvalues of $R$ and $\psi_1(\lambda)$ is the first component of the corresponding normalized eigenvector; the reciprocals of the eigenvalues are the Yang-Lee zeros, i.e. the zeros of the grand-canonical partition function in the complex activity plane. Singularities of $\rho$ occur where $-z^{-1}$ meets an edge of the R-matrix spectrum, and the local behavior of $\psi_1$ at an edge fixes the exponent of the density singularity at that termination point, defined by $\rho(z_c)-\rho(z)\propto(z_c-z)^{\sigma'}$ at the physical termination point and analogously with $\sigma$ at the nonphysical one. For every repulsive model studied the matrix elements converge to two constants, $R_{n,n}\to B$ and $R_{n,n+1}\to A$, which yields termination points $z_t=(2A-B)^{-1}$ and $-z_0=-(2A+B)^{-1}$, with square-root exponents; corrections to the constants modify the exponents. When the corrections decay as $1/n^2$, a continuum (Bessel-type) limit of the eigenvalue equation gives closed exponents (Eisenberg and Baram 2007). Rotman and Eisenberg (2009) treated the general asymptotic form

$R_{n,n}=B+\frac{b}{n^2}+\frac{b'\cos(qn)}{n},\qquad R_{n,n+1}=A+\frac{a}{(n+1/2)^2}+\frac{a'\cos(q(n+1/2))}{n+1/2},$

by a real-space decimation (renormalization) argument, obtaining the closed exponent formula

$\sigma=\frac{1}{2}\sqrt{1-\frac{4(2a+b)}{A}-\frac{[2a'\cos(q/2)+b']^2}{[1-\cos(q)]A^2}},\qquad \sigma'=\frac{1}{2}\sqrt{1-\frac{4(2a-b)}{A}-\frac{[2a'\cos(q/2)-b']^2}{[1-\cos(q)]A^2}},$

which reduces to the earlier result when the oscillatory amplitudes vanish. The formula was validated against the exactly solved hard hexagons model ($\sigma'=0.6662$ versus the exact $2/3$; nonphysical $\sigma=0.1655$ versus the universal $1/6$) and against synthetic matrices with prescribed asymptotics. Applied to the square-lattice gases in which each particle excludes all sites within the first $k$ neighbor shells (the N$_k$ models, studied by Monte Carlo by Fernandes, Arenzon and Levin 2007), it performed well for N5 (no oscillations, $1/n^3$ corrections being irrelevant in the renormalization sense) but exposed its own limits, which are the subject of this problem: for the triangular-lattice gas with first- and second-neighbor exclusion (believed from Monte Carlo to lie in the four-state Potts universality class, a standard exactly solved spin-model family; critical chemical potential $\mu_c=1.75682(2)$, i.e. critical activity $z_t=e^{\mu_c}=5.794$; Zhang and Deng 2008), the fitted parameters give $4(2a-b)/A\simeq 2.7>1$, so $\sigma'$ is formally imaginary, the spectrum edge shifts from $2A-B$, and the predicted $z_t=(2A-B)^{-1}=9.35$ is far from the Monte Carlo value; for the square-lattice N4 model the method gave $\sigma'=0.28(6)$, contradicting the Ising assignment of the simulations (a contradiction since resolved in favor of the simulations by Nath and Rajesh 2014); and the construction itself rests on the unproven assumption that the R-matrix elements are real for every purely repulsive model.

Problem Statement

Extend the exact R-matrix analysis of the Mayer cluster series beyond its established range of validity, so that the extraction of termination points and critical exponents from cluster coefficients works, and is provably applicable, for the models on which it currently fails. Three linked results are required, corresponding to the source paper's stated limits of usage of the exponent formula:

(1) The spectrum-shift regime. Consider R matrices whose elements approach their limiting constants with the asymptotic corrections

$R_{n,n}=B+\frac{b}{n^2}+\frac{b'\cos(qn)}{n},\qquad R_{n,n+1}=A+\frac{a}{(n+1/2)^2}+\frac{a'\cos(q(n+1/2))}{n+1/2},$

where $A$ and $B$ are the limiting constants, $a$, $b$, $a'$, $b'$ the correction amplitudes, and $q$ the oscillation wavenumber. When the correction amplitudes are so strong that the radicand of the exponent formula is negative — for the physical edge, $4(2a-b)/A+[2a'\cos(q/2)-b']^2/([1-\cos(q)]A^2)>1$, so that the formal $\sigma'$ is imaginary — the constant-matrix spectrum edges $2A+B$ and $2A-B$ no longer bound the spectrum: the spectrum edge shifts, the termination points are no longer $(2A\pm B)^{-1}$, and the exponent formula no longer applies. Determine, for this strong-correction regime, the location of the shifted spectrum edge, the resulting physical termination activity $z_t$, and the exponent $\sigma'$ of the density singularity at it. The triangular-lattice gas with exclusion of first and second neighbors is the physical instance (R-matrix fit $z_t=9.35$ versus Monte Carlo $z_t=5.794$), but the question concerns the general asymptotic class of matrices, not only this model.

(2) Reality of the R matrix. Prove that the recursively constructed R-matrix elements are real for every purely repulsive model (hard-core lattice gases and purely repulsive continuum gases, the classes studied in the source), as observed in every model computed to date, or exhibit an explicit counterexample.

(3) Error-controlled extraction. Develop a robust procedure for determining the physical termination activity $z_t$ — equal to $(2A-B)^{-1}$ where the constant-matrix edges bound the spectrum and given by the shifted value of (1) in the strong-correction regime — and the exponents, from finite-order fitted parameters, with quantified error bars, given that small errors in $A$ and $B$ are strongly amplified in $z_t$ and that the R-matrix construction is itself exponentially sensitive to errors in the cluster integrals.

These three results jointly complete the method: (1) extends the exact treatment to the regime where it demonstrably fails, (2) secures the hypothesis on which the entire construction rests, and (3) makes its numerical output trustworthy on the finite cluster data any real model provides. The existing exponent formula was itself established by a combination of analytic derivation and numerical validation, recovering the exact hard-hexagon benchmarks ($z_c=(11+5\sqrt{5})/2$, $\sigma'=2/3$) and documenting the termination-point discrepancy of the triangular first- and second-neighbor exclusion model; since the hard hexagons lies inside the validity range of the existing formula, the spectrum-shift-regime triangular-lattice model provides a genuinely new test of any extension.

The square-lattice N4 universality dispute (exclusion of the first four neighbor shells) is out of scope: the question whether that transition is Ising-like, contrary to the R-matrix prediction $\sigma'=0.28(6)$, has been settled by later Monte Carlo work (Nath and Rajesh 2014: the transition is Ising, and the model has a second, Ashkin-Teller-type transition, in the universality family of coupled Ising spin pairs). The model remains a natural consistency check, since a correct solution of (1)-(3) should account for why the R-matrix analysis failed there; diagnosing that failure is not a separate deliverable.

References

1. Z. Rotman and E. Eisenberg, 'Critical exponents from cluster coefficients', Phys. Rev. E 80, 031126 (2009). DOI: 10.1103/PhysRevE.80.031126; arXiv:0907.4037. https://doi.org/10.1103/PhysRevE.80.031126
2. E. Eisenberg and A. Baram, 'Estimation of critical exponents from cluster coefficients and the application of this estimation to hard spheres', Proc. Natl. Acad. Sci. USA 104, 5755-5758 (2007). DOI: 10.1073/pnas.0700778104. https://doi.org/10.1073/pnas.0700778104
3. R. J. Baxter, 'Hard hexagons: exact solution', J. Phys. A: Math. Gen. 13, L61-L70 (1980). DOI: 10.1088/0305-4470/13/3/007. https://doi.org/10.1088/0305-4470/13/3/007
4. G. S. Joyce, 'On the hard-hexagon model and the theory of modular functions', Phil. Trans. R. Soc. Lond. A 325, 643-702 (1988). DOI: 10.1098/rsta.1988.0077. https://doi.org/10.1098/rsta.1988.0077
5. A. Baram and M. Luban, 'Universality of the cluster integrals of repulsive systems', Phys. Rev. A 36, 760-765 (1987). DOI: 10.1103/PhysRevA.36.760. https://doi.org/10.1103/PhysRevA.36.760
6. H. C. M. Fernandes, J. J. Arenzon and Y. Levin, 'Monte Carlo simulations of two-dimensional hard core lattice gases', J. Chem. Phys. 126, 114508 (2007). DOI: 10.1063/1.2539141. https://doi.org/10.1063/1.2539141
7. W. Zhang and Y. Deng, 'Monte Carlo study of the triangular lattice gas with first- and second-neighbor exclusions', Phys. Rev. E 78, 031103 (2008). DOI: 10.1103/PhysRevE.78.031103; arXiv:0803.3334. https://doi.org/10.1103/PhysRevE.78.031103
8. T. Nath and R. Rajesh, 'Multiple phase transitions in extended hard-core lattice gas models in two dimensions', Phys. Rev. E 90, 012120 (2014). DOI: 10.1103/PhysRevE.90.012120; arXiv:1404.6902. https://doi.org/10.1103/PhysRevE.90.012120
9. Z. Rotman and E. Eisenberg, 'Ideal glass transition in a simple two-dimensional lattice model', Phys. Rev. E 80, 060104(R) (2009). DOI: 10.1103/PhysRevE.80.060104; arXiv:0908.2595. https://doi.org/10.1103/PhysRevE.80.060104
10. Z. Rotman and E. Eisenberg, 'Direct measurements of the dynamical correlation length indicate its divergence at an athermal glass transition', Phys. Rev. Lett. 105, 225503 (2010). DOI: 10.1103/PhysRevLett.105.225503; arXiv:1008.3688. https://doi.org/10.1103/PhysRevLett.105.225503
11. A. A. A. Jaleel, D. Mandal and R. Rajesh, 'Hard core lattice gas with third next-nearest neighbor exclusion on triangular lattice: One or two phase transitions?', J. Chem. Phys. 155, 224101 (2021). DOI: 10.1063/5.0066098. https://doi.org/10.1063/5.0066098
12. A. A. A. Jaleel, D. Mandal, J. E. Thomas and R. Rajesh, 'Freezing phase transition in hard-core lattice gases on the triangular lattice with exclusion up to seventh next-nearest neighbor', Phys. Rev. E 106, 044136 (2022). DOI: 10.1103/PhysRevE.106.044136; arXiv:2206.04985. https://doi.org/10.1103/PhysRevE.106.044136
13. S. Darjani, J. Koplik, S. Banerjee and V. Pauchard, 'Liquid-hexatic-solid phase transition of a hard-core lattice gas with third neighbor exclusion', J. Chem. Phys. 151, 104702 (2019). DOI: 10.1063/1.5123231; arXiv:1908.05555. https://doi.org/10.1063/1.5123231