Preprint
Article

Bifurcation of Solutions of Ornstein-Zernike Equation in the Gas-Liquid Coexistence Region of Simple Fluids

This version is not peer-reviewed.

Submitted:

23 February 2024

Posted:

23 February 2024

You are already at the latest version

A peer-reviewed article of this preprint also exists.

Abstract
One objective of this paper is to investigate bifurcation of solutions of Ornstein-Zernike equation in liquid state theory for the gas-liquid coexistence region of simple fluids. The analysis uses a discrete form of the Ornstein-Zernike equation in real space, together with a closure relation which provides nearly self consistent thermodynamic properties of fluids governed by interaction potentials like the Lennard-Jones potential. Accurate asymptotic forms of correlation functions are incorporated in the real space algorithm so as to mitigate the influence of cutoff length used in the discrete approximations. It is found that there is no spinodal on the vapor side on a sub-critical isotherm. But there are fold bifurcation points on the vapor and liquid sides together with a 'no-real-solution-region' in between. This is followed by a region of negative compressibility extending up to the spinodal on the liquid side. The emerging picture is quite similar to that in the case of the hypernetted-chain-closure. It is also found that complex solutions connect the physical solutions that exist outside the 'no-real-solution-region' as well as in region of negative compressibility. The second objective of the paper, in the light of bifurcation scenario, is to evolve an interpretation of the coupling-parameter expansion of correlation functions pertaining to the region of negative compressibility. It is shown that the expansion provides an accurate description of results of the restricted ensembles, implied in mean field theories and simulation methods. Bifurcation of singlet density to inhomogeneous distributions at a spinodal is also established using a simple equation, which uses the direct correlation function. Thus the paper reestablishes the fact that restriction to a homogeneous singlet density is the root cause of all inconsistencies found in the gas-liquid coexistence region.
Keywords: 
Ornstein-Zernike equation; bifurcation of solutions; fold bifurcation points; liquid-vapor phase transition; coupling parameter expansion; Lennard-Jones potential
Subject: 
Physical Sciences  -   Condensed Matter Physics

1. Introduction

Ornstein-Zernike equation (OZE) and related developments form the main theoretical apparatus in liquid state physics; the other being computer simulation methods [1]. In addition to being an important branch of statistical mechanics, liquid state theory is an essential component in modeling the high-pressure equation of state [2], with applications in several branches of solid state science and high energy-density physics. So it is important to have quantitative knowledge about the physical and mathematical nature of solutions of OZE and the associated closure relations. Hence this topic has been the subject of detailed investigations over the last few decades. While the OZE equation is a relation defining the direct correlation function c ( r ) in terms of the total correlation functions h ( r ) , a closure relation relates them to the inter-particle potential U ( r ) . Common examples of closures include the Mean spherical approximation (MSA), Percus-Yevick (PY) closure and hypernetted chain (HNC) closure. There are also some approaches consisting of adjustable parameters to achieve self consistency, which provide same results via different routes to thermodynamics [3]. There exist also quite accurate, nearly self-consistent and explicit closures for specific potentials like hard sphere, the Lennard-Jones (LJ) potential, etc. [4,5].
A closure relation is solvable (in principle) so as to express c ( r ) in terms of h ( r ) explicitly, although this relation could be highly nonlinear. On using this relation, the OZE reduces to a nonlinear mapping of the form: h = O h for h ( r ) [6], and its solutions are simply the fixed-points of O , which acts on functions defined in a suitable metric space. The specific details of the mapping O will depend on the closure relation and inter-particle potential. Equally important, changes in fixed-point stability and bifurcation of solutions, resulting due to variation of parameters like density ρ and temperature T occurring in O , will also depend on the discrete approximations to O , and the algorithm chosen to obtain solutions. For instance, algorithms may use 3D Fourier transforms for computing convolution integrals or simply employ real-space based methods. This is an important aspect to be considered when multiplicity and bifurcation of solutions to nonlinear eigenvalue problems are investigated [7].
As the correlation function h ( r ) decays for sufficiently large r, it is reasonable to assume that the constraint I = h ( r ) d r < is satisfied in physically relevant situations. In fact, this integral constraint (IC) is a necessary condition for computing all thermodynamic quantifies and structure factor, as fluctuation theorem [1] relates it to reduced compressibility χ as χ = 1 + ρ I . Further, violation of the IC in any part of phase-plane would generate problems in algorithms using 3D Fourier transforms, as the IC is necessary to ensure the existence of the transforms. This is important as many algorithms for solving OZE compute the convolution integral using Fourier transforms [8,9].
Baxter derived an alternate formulation of OZE [10], defined entirely in a finite range 0 r R 0 , where R 0 is a cutoff length introduced such that c ( r ) = 0 for r > R 0 . The vanishing property of c ( r ) is exactly valid within the PY closure when potential U ( r ) = 0 for r > R 0 . This formulation is derived by imposing the IC mentioned above, and is also useful even in cases where c ( r ) is negligible for r > R 0 where R 0 is chosen suitably. Approaches based on either Baxter’s formulation [11] or 3D Fourier transforms [12], have been used extensively to analyze bifurcation of solutions of OZE. In the latter case, a large cutoff length R 1 is introduced for use of discrete variables and FFT techniques. Violation of the IC may not be directly visible with any choice of R 1 , however large. So it is essential to introduce appropriate asymptotic expressions of correlation functions and their transforms, valid for r > R 1 , in the formulation [12]. Furthermore, this ‘OZE technique’ is unavoidable to obtain the correct variation of transforms in the Fourier variable k near its origin [12].
The OZE can also be written in the radial coordinate covering the full half-line 0 r < , and solved recursively by applying Newton’s method [13]. The linear Fredholm integral equations (FIE) of second kind, which result within this approach, are readily solvable using a proper numerical scheme [14]. Of course, a discrete representation of variable r, sufficiently large cutoff length R 1 , and use of asymptotic expressions of functions beyond R 1 are needed in this approach (see §3). Note that numerical schemes for FIE are well developed, and there are also ways of reducing the half-line to finite intervals using appropriate transformations [15]. Adequacy of R 1 can be checked by generating a sequence of solutions (and integral quantities) for an increasing series of values of R 1 . The Converged solutions obtained with this approach need not satisfy the IC mentioned earlier. Analyses of solutions for the three closures, obtained with a finite element method (using piece wise linear functions) to solve the FIE in 0 r R 1 , together with asymptotic conditions h ( r ) = c ( r ) = 0 for r > R 1 , have been done quite early [13]. Sensitivity of solutions to the chosen R 1 is also analyzed in this extensive work. There are differences in the results for some problems derived using this real-space method and Fourier transform method (see §2.1). This is plausible because the solutions obtained by the latter belong to a sub-space of the more general space of functions accounted within the former.
One of the objectives in the present paper is to investigate bifurcation of solutions of OZE combined with the Martynov-Sarkisov (MS) closure for LJ potential (see §2). This closure is chosen as it provides nearly self consistent results for hard sphere and LJ fluids [4]. A general real-space method, using piece wise quadratic functions for numerical integration, together with improved asymptotic conditions [16], in lieu of the simple ones used earlier [13], are employed in the present work (see §5). This mitigates the ill effects of chosen value of R 1 on the solutions obtained. The resulting real-space algorithms is first verified using an analytical solution [13]. Fold bifurcation points (turning points), on the vapor as well as liquid sides on any sub-critical isotherm, are first found using the pseudo arc length (PAL) continuation algorithm [17,18]. A ‘no-real-solution’ region, a non-physical branch of solution, regions of negative compressibility χ , and spinodal points are all obtained with the method. Further, bifurcation of solutions to the complex domain is analyzed [19], and the method is generalized to compute complex solutions in the gas-liquid coexistence region, wherever they exist (see §7.2).
In the light of the bifurcation scenario for OZE and MS closure, questions arise regarding the meaning of ‘solutions’ obtained via the coupling-parameter expansion (CPE) [20]. This approach uses expansion of correlation functions in terms of a parameter introduced in the potential, and provides real solution in the entire gas-liquid coexistence region [21]. The reference solution corresponds to purely repulsive part of the potential. Spinodal points, on vapor and liquid sides, with negative compressibility in between, are obtained on any sub-critical isotherm. The IC is violated in that region (see §2), thereby indicating some inconsistency in the CPE. However, the gas-liquid phase diagrams for Mie and LJ fluids, obtained via the compressibility route employing Maxwell’s constructions, agree with Monte Carlo simulation data very well. It also agrees excellently with results of an independent method (using thermodynamic conditions) based on free energy computed outside spinodal region [22]. It is important to emphasize that the second method does not use any data corresponding to the negative-compressibility-region. It is now shown (see §8.1) that data on isotherms through the spinodal region agree with those of transition-matrix Monte Carlo method [23]. Therefore, it is proposed (in this paper) that results of CPE are similar to those obtained via the restricted ensemble approach [13]. It is well known that the van der Waals equation, which yields negative compressibility region, is also derived from considerations of a restricted ensemble [24]. In addition, analysis of a simple equation for singlet density (see §8.2), involving the direct correlation function, shows that it bifurcates at spinodal points to an in-homogeneous distribution. This finding explains the origin of inconsistency arising out of negative compressibility obtained in any method.
Remainder of the paper is organized as follows. Basic equations and formulas used in the model are given in §2. Newton’s method in real and Fourier spaces, and some known results on bifurcation, are discussed in §3. Bifurcation of solutions (using real-space methods for MS closure ) inside the gas-liquid coexistence region are discussed in §4. Emergence of an infinite number of solutions in the spinodal region due to ‘singularity-induced bifurcation’ is bought out here. Discrete approximations of the FIE, arising in Newton’s method, using Nyström’s method are discussed in §5, while continuation algorithms form the topic of §6. Generalization of Newton’s method for complex solutions are covered in §7. The coupling parameter expansion and its relevance for restricted ensemble approach are given in §8. These sections (§6, §7, §8 ) also provide numerical results for LJ potential within MS closure. Finally, §9 summarizes the main conclusions of the present work.

2. Equations and formulas of the fluid model

The OZE introduces the short ranged direct correlation function c ( r ) via a convolution integral involving the total correlation function h ( r ) . For a one-component system of fluid particles, interacting via spherically symmetric potentials, it is written as
h ( r ) = c ( r ) + ρ c ( | r r | ) h ( r ) d r ,
where ρ is the average number density of particles. This equation needs to be supplemented with a closure relation involving h ( r ) , c ( r ) and the pair-potential U ( r ) . There is no exact analytical expression for the closure relation; however, all the basic closures are generally expressed concisely in terms of a ‘bridge function’ B ( r ) via the expression [1]
c ( r ) = exp [ β U ( r ) + y ( r ) + B ( r ] y ( r ) 1 ,
where β 1 = ( k B T ) , k B is Boltzmann’s constant and T absolute temperature. Further, y ( r ) = h ( r ) c ( r ) denotes the indirect correlation function. So different approximate forms of B ( r ) generate different closures; for instance, the HNC closure is obtained with B ( r ) = 0 . The PY closure follows from B ( r ) = ln [ 1 + y ( r ) ] y ( r ) , which yields c ( r ) = exp [ β U ( r ) ] 1 ( y ( r ) + 1 ) . Note that PY relation also follows from HNC closure by using the linear approximation exp [ y ( r ) ] 1 + y ( r ) . The MSA closure is defined for potentials with a hard core, and it takes the form c ( r ) = β U ( r ) for r > σ (the hard core diameter) together with the hard core condition g ( r ) = h ( r ) + 1 = 0 for r σ .
A bridge function, which generates the Martynov-Sarkisov (MS) closure ( mentioned in §1) is given by [4]
B ( r ) = 1 + 2 ω ( r ) ω ( r ) 1 , ω ( r ) = y ( r ) ρ σ 3 β U a ( r ) .
The MS closure provides accurate thermodynamic properties for hard sphere as well as LJ potential [4]. The full LJ potential, which will be used throughout this paper, is given by
U ( r ) = ϵ [ ( σ / r ) 12 ( σ / r ) 6 ] ,
where ϵ and σ are energy and length parameters, respectively. Note that ρ σ 3 in Eq.(3) is a dimensionless number-density as σ is a length parameter. Further, U a ( r ) in Eq.(3) is the pure attractive part of the potential as per the Weeks-Chandler-Andersen (WCA) prescription [25].
It is important to note that any solution of OZE must satisfy the relation [ h ( r ) c ( r ) ] d r = ρ [ c ( r ) d r ] [ h ( r ) d r ] , which follows by integrating it over the infinite space. Then, the integral condition (IC) I = h ( r ) d r < , introduced in §1, turns out to be necessary to deduce an alternate relation for reduced inverse compressibility, viz. χ 1 = β P / ρ = 1 ρ c ( r ) d r , where P is pressure. Therefore, the IC must be violated in regions in the phase-plane where χ 1 0 . This is important because there can be situations where h ( r ) slowly decays to zero as r , but I does not exist. Methods which assume that I exists (Baxter’s version of OZE or OZE in Fourier space) can not work reliably in such regions. So an algorithm using Fourier transforms, and a finite cutoff length R 1 , may give negative χ 1 in some region, but those results have to be suspected [12].
Detailed analysis of the asymptotic variation of the correlation functions shows that h ( r ) varies for large r as [16]
h ( r ) χ 2 β U ( r ) + χ ξ 2 1 4 π exp [ r / ξ ] r , C 0 = 0 c ( r ) d r , C 2 = 1 3 ! 0 c ( r ) r 2 d r ,
where ξ = ρ C 2 / ( 1 ρ C 0 ) is the correlation length. This asymptotic form is valid for χ > 0 . The correlation length is also expressed as ξ ρ χ , χ 1 = 1 ρ C 0 , and the amplitude A = χ / 4 π ξ 2 is positive. For the HNC closure, as easily derived from Eq.(2), c ( r ) varies as
c ( r ) ( H N C ) β U ( r ) + 1 2 h 2 ( r ) ,
so that it is a short ranged function compared to h ( r ) . For the MS closure (due to the bridge function), the asymptotic form for c ( r ) is found to be
c ( r ) ( M S ) β U ( r ) + 1 6 h 3 ( r ) ,
which decays faster than that in the case of HNC closure. This feature has some consequences in determining the solution in the spinodal region. In the present paper, asymptotic forms either in Eqs.(5) and (7) or via Nyström’s interpolation formula (see §5) are used for the range r > R 1 .
When c ( r ) decays as r q for large r, it is necessary that q > 3 for χ 1 to exist. In these cases h ( r ) will decay as r q / 2 for HNC closure. This decay law do ensure the existence of the integral in Eq.(1), however, the IC will be violated for the range 3 < q 6 . Similar conclusion follows for MS closure in the range 3 < q 9 . While real-space methods would provide solutions in such cases, methods using Fourier transforms may face numerical problems.

2.1. Known results of bifurcation of solutions

Known results on bifurcation of solutions of OZE, on a sub-critical isotherm in the two-phase region, may be classified according to PY, HNC and MS closures. The following brief discussion uses ρ s v and ρ s l to denote densities at spinodal points on vapor and liquid sides, respectively. Similarly, ρ t v and ρ t l are used for the corresponding turning points (fold bifurcation points).
One clear result is from Baxter’s adhesive hard sphere model [26,27], which treats a system of sticky hard spheres within PY closure; the stickiness being modeled by an infinitely deep but vanishingly thin attractive square well. Analytical solutions of Baxter’s version of OZE (which assumes the IC mentioned in §1) show that ρ s v does not exist. In fact, there is a ρ s v for which χ 1 = 0 , however, it lies on a solution-branch violating the IC. Then, the densities between ρ t v and ρ t l cover a ‘no-real-solution’ region. As density is increased further, another non-physical branch (violating the IC ) follows, which terminates at ρ s l . The model has complex solutions within the ‘no-real-solution’ region.
Numerical treatment of Baxter’s OZE and PY closure, applied to a truncated LJ potential [11], yields both ρ t v and ρ t l and ‘no-real-solution’ in between. This work shows that, as in Baxter’s model, χ 1 is finite at ρ t v as well as ρ t l . Methods using Fourier transforms and PY closure, also come to similar conclusions, either by incorporating proper asymptotic forms of correlation functions [12] or using large cutoff length ( 165 σ ) [28]. Neither the former work (for hard-core Yukawa potential) nor the latter (for full LJ potential) report ρ t l values. This could be because χ 1 becomes negative in the range ρ t l ρ ρ s l . The implicit assumptions involved in using Fourier transforms [12] show that this method should be restricted to regions where χ 1 > 0 , where the IC is satisfied. Complex solutions have also been found for LJ potential near ρ t v [28]. Solutions of OZE and PY closure in real space [13], also show similar results, i.e., evidence for ρ t v and ρ s l , but no ρ t l . In addition, the last work also finds density regions (between ρ t v and ρ s l ) with negative compressibility.
For HNC closure, all approaches find ρ t v and ρ t l densities. Locus of these values (for different T), have been found using Fourier transforms and pseudo arc length (PAL) continuation algorithm [29]. Here it is also shown that χ 1 values on the liquid side are very sensitive to cutoff length R 1 - a lower value of R 1 ( 20 σ ) yields ρ s l , but it is absent with R 1 80 σ . Importantly, the densities ρ s v and ρ s l are not obtained at all when proper asymptotic forms are incorporated into Fourier transforms[12]. Similar results together with complex solutions, near ρ t v and ρ t l , are also obtained using Fourier transforms [28]. But results of real space solutions [13] clearly show the occurrence of both pairs of densities ( ρ t v , ρ t l ) and ( ρ s v , ρ s l ) . In addition, regions of negative χ 1 are also obtained [13] - however, ρ s v is on non-physical branch and ρ t l occur in the region of negative χ 1 . Thus, for the HNC closure, there are some differences between these results and those obtained using Fourier transforms. In fact, results from real space methods are closer to those of Baxter’s model.
The occurrence of ρ t v has been shown with several other closures, including the MS closure [30]. Thus it appears that a turning point on the vapor side of a sub-critical isotherm is found with all closures as well as methods. Emergence of complex solutions near ρ t v , and occurrence of ρ s l , have also been shown for MS closure using Fourier transform method. One of the aims in this paper is to investigate the complete bifurcation scenario within MS closure using real space methods.
Finally, note that analytical solutions exist to the hard-core Yukawa potential, within the MSA closure [31,32], which show the emergence of two real and two complex solutions in the entire range of ρ and T in the super-critical region. Out of these, only one real solution approaches the ideal gas limit for small ρ , and vanishes as the interaction strength K 0 . In the sub-critical region, there is a ’no-real-solution-region’ for a range of densities ρ 1 < ρ < ρ 2 , because both solutions turn complex. Within this model, the spinodal lines, which now occur on both the gaseous and liquid branches, cover the ’no-real-solution-region’, and χ 1 is always non-negative (as in Baxter’s model).

3. Newton’s method in real and Fourier spaces

As the indirect correlation function y ( r ) is more appropriate for numerical approximations, first Eq.(1) is rewritten in term of y ( r ) as
y ( r ) = ρ c ( | r r | ) [ c ( r ) + y ( r ) ] d r ,
Newton’s method in real space 0 r < applied to Eqs.(8) and (2), with bridge function B ( r ) specified in Eq.(3, is the approach used in this paper. Such an approach has been applied earlier, however, with PY and HNC closures, for a truncated and shifted form of the LJ potential [13]. Solutions to another well known non-linear equation in fluid physics, known as Born-Green-Yvon-Kirkwood (BGKY) equation, is also analyzed in this reference. This work used the asymptotic conditions h ( r ) = c ( r ) = 0 for r > R 1 , with largest value R 1 = 12 σ . A finite element method, employing piece wise linear functions, was used to solve the equations in 0 r R 1 . A more general algorithm and several improvements, introduced in the present work, are discussed in this section.
Starting with an initial guess solution y ( r ) (for specified ρ and T), and hence c ( r ) from Eqs.(2) and (3), improved approximations c ( r ) + Δ c ( r ) and y ( r ) + Δ y ( r ) are obtained using the first order correction term
Δ c ( r ) = [ h ( r ) + g ( r ) B ( r ) ] Δ y ( r ) ,
where h ( r ) = g ( r ) 1 = y ( r ) + c ( r ) and B ( r ) = 1 / 1 + 2 ω 0 1 is the derivative of B [ y ] with respect to y. The correction term Δ y ( r ) obeys the linear equation
Δ y ( r ) = ρ [ c ( | r r | ) Δ y ( r ) + { 2 c ( | r r | ) + y ( | r r | ) } Δ c ( r ) ] d r s ( r ) ,
where s ( r ) is the error (or residual ) in Eq.(8) with the starting guess solutions y ( r ) and c ( r ) , viz.
s ( r ) = y ( r ) + ρ c ( | r r | ) [ c ( r ) + y ( r ) ] d r .
Substituting Δ c ( r ) from Eq.(9) gives the linear Fredholm integral equation (FIE)
Δ y ( r ) = ρ [ c ( | r r | ) + q ( r ) { 2 c ( | r r | ) + y ( | r r | ) } ] Δ y ( r ) d r + s ( r ) ,
where q ( r ) = [ h ( r ) + g ( r ) B ( r ) ] . First step in Newton’s method consists of solving this FIE, and obtaining an improved solution y ( r ) y ( r ) + Δ y ( r ) . The process is then repeated with new y ( r ) as the guess solution until certain convergence criterion is met. Usually, these iterations are terminated when a suitable norm | | Δ y | | (say, Euclidean) of Δ y ( r ) is less than a prescribed value (say, 10 6 ). This way of implementing Newton’s method is sufficiently general and is applicable with other closure relations; either B ( r ) is to be changed or an appropriate form of q ( r ) is to be substituted. For instance, in the case of HNC closure B ( r ) = 0 while q ( r ) = exp [ β U ( r ) ] 1 for PY closure.
On writing | r r | = r 2 + r 2 2 r r μ , where μ is the cosine of the angle between r and r , the integral over μ in the interval [-1, 1] can be rewritten in r and r . Then the terms in the kernel of Eq.(12) are given by
C ( r , r ) = 1 2 r r | r r | r + r c ( t ) t d t ,
Y ( r , r ) = 1 2 r r | r r | r + r y ( t ) t d t .
With these explicit definitions, the kernel of the FIE reduces to
K ( r , r ) = C ( r , r ) + q ( r ) { 2 C ( r , r ) + Y ( r , r ) } .
Using same notations, the error term is given by
s ( r ) = y ( r ) + 4 π ρ 0 C ( r , r ) { c ( r ) + y ( r ) } r 2 d r .
Now the FIE for Δ y ( r ) is concisely expressed as
Δ y ( r ) = 4 π ρ 0 K ( r , r ) Δ y ( r ) r 2 d r + s ( r ) ,
which can be solved using standard methods [14]. Quadrature formulas to approximate the integral for a specified value of r, and collocation rules form the tools in Nyström’s method to reduce a linear FIE to a matrix equation [15].

3.1. Version in Fourier space

If the integral constraint I < is valid, it is easy to rewrite Eq.(12) in Fourier space. On taking Fourier transforms it reduces to
Δ y ¯ ( k ) = ρ 2 c ¯ ( k ) + y ¯ ( k ) 1 ρ c ¯ ( k ) [ q ( r ) Δ y ( r ) ] ¯ + [ ρ c ¯ ( k ) 2 1 ρ c ¯ ( k ) y ¯ ( k ) ] .
Here k is the Fourier variable, and a ‘bar’ over the functions indicates their 3D Fourier transforms [22]. Note that first term on the right of Eq.(12) has appeared as 1 ρ c ¯ ( k ) in the denominators. The integral operator in this equation is defined in terms of Fourier transform. It is best solved using Krylov space-based methods like CGNR (Conjugate Gradient to Normal equations to minimize Residual) or GMRES (Generalized Minimum Residual), because the needed matrix-vector products in these methods can be computed via two Fourier transform operations [22]. To employ Fast Fourier Transform (FFT) techniques, first the radial coordinate is truncated to 0 r R 1 with a sufficiently large value of R 1 . It is also essential to use asymptotic forms of c ( r ) and y ( r ) , and their Fourier transforms, to develop a robust algorithm [12]. This method would develop numerical problems whenever the IC is violated as is the case close to spinodal points (if they exist), and on non-physical solution branches where compressibility χ is negative.

4. Bifurcation of solutions

Solving the linear FIE given in Eq.(17) and performing subsequent Newton’s iterations would provide a solution for one set of ρ and T, if a suitable starting solution is provided. However, it is necessary to generate solutions for different values of ρ on each isotherm, and on a number of isotherms. Bifurcation or branching of solutions occur at certain values of ρ and T whenever new solutions appear. Bifurcation analysis, and continuation methods for tracing the different branches along ρ , form a well developed topic. Excellent books cover all the mathematical and computational aspects of this fascinating field in nonlinear physics [17,18]. Applications of these methods for OZE and MS closure in real space are discussed in §6; the present section is devoted to a brief discussion of different branching scenarios.
Suppose y ( r ) is a solution of Eqs.(1) and (3) corresponding to density ρ . Then, the change in the solution Δ y ( r ) , due to a specified change Δ ρ in density, satisfies the linear equation
Δ y ( r ) ρ [ c ( | r r | ) + q ( r ) { 2 c ( | r r | ) + y ( | r r | ) } ] Δ y ( r ) d r + Δ ρ c ( | r r | ) [ c ( r ) + y ( r ) ] d r = 0 .
This equation follows exactly via the the same method elaborated in §3 for Newton’s method. For simplicity of notation, terms arising due to explicit dependence of c ( r ) on ρ are omitted. Once again, q ( r ) = [ h ( r ) + g ( r ) B ( r ) ] , which corresponds to the solution at specified ρ . Note that there is no term s ( r ) in the equation above, as occurring in Eq.(12), because y ( r ) satisfies Eqs.(1) and (3). The last equation may be written in operator form as
( I L ) Δ y + Δ ρ F = 0 ,
where ( I L ) is a linear operator and Δ y is the change in solution vector. Further, F is the operator defining the OZE and 0 a null vector. For a sufficiently small Δ ρ , solution of Eq.(20) will yield an approximate solution y a ( r ) = y ( r ) + Δ y ( r ) for the density ρ + Δ ρ . Then, starting with y a ( r ) as initial guess, Newton’s method is applied to obtain a new solution-pair { y n e w , ρ n e w } . In this way, a set of pairs { y j , ρ j } is generated as points defining a solution-curve. However, in the neighborhoods and at values of ρ where ( I L ) is singular (called critical points), this approach is bound to fail. In fact, the singularity of the linear operator ( I L ) signals bifurcation of solutions.
Analysis of bifurcation of solutions is usually done together with labeling the solution-curve using the arc length parameter . Thus, = 0 corresponds to the initial point { y 0 , ρ 0 } where y 0 is known, and Eq.(20) is considered as a relation between Δ y and Δ ρ . For evolving the pair { y , ρ } as functions of , Eq.(20) is modified as
( I L ) y ˙ + ρ ˙ F = 0 .
Here y ˙ and ρ ˙ denote the components of the tangent vector to solution-curve at . The nature of the singularity of the operator ( I L ) , and the relationship of its null-space to F , decide the type of bifurcation that occurs at that ρ . These aspects will be covered in §6, using a discrete approximation to operator L in Eq.(21). However, there is another type of bifurcation which is considered below.

4.1. Singularity induced bifurcation

The operator ( I L ) can be split into two parts, and Eq.(21) rewritten as
( I L 1 L 2 ) y ˙ + ρ ˙ F = 0 .
Here the part ( I L 1 ) y ˙ , for a specified density ρ and corresponding c ( r ) , is defined as
( I L 1 ) y ˙ = y ˙ ( r ) ρ c ( | r r | ) y ˙ ( r ) d r ,
and L 2 denotes the remaining part. Note that L 2 , which contains the function q ( r ) in its kernel, arises from the change in c ( r ) induced (via Δ y ) due to the change Δ ρ . If ( I L 1 ) is non-singular, then Eq.(22) can be rewritten as [ I L 2 ( I L 1 ) 1 ] z ˙ + ρ ˙ F = 0 , where the new vector z ˙ is defined as ( I L 1 ) y ˙ = z ˙ . Either y ˙ or z ˙ can be used for the solution-curve as far as ( I L 1 ) is non-singular. However, at a value of ρ = ρ b where ( I L 1 ) is singular, the matrix [ I L 2 ( I L 1 ) 1 ] will diverge, and the associated changes in solution is termed as ‘singularity induced bifurcation’ [33]. But as ( I L 1 L 2 ) is non-singular at ρ b , the vectors y ˙ and hence z ˙ are well defined, and these will vary smoothly across ρ b . So Newton’s method (via Eq.(12) or (17) ) will yield a solution y ( r ) to the original OZE in Eq.(8) at ρ b . Singularity induced bifurcation is associated with the divergence of the operator (or matrix) involved in one of the routes towards the solution.
Now, note that the operator involved in Eq.(1) or (8) - with the converged function c ( r ) - is the same as ( I L 1 ) . Therefore, y ( r ) + α ϕ ( r ) is also a solution at ρ b , where α is an arbitrary constant and ϕ ( r ) the eigenfunction corresponding to the zero-eigenvalue of ( I L 1 ) . It is easily verified by direct substitution that the singularity occurs at ρ b whenever 1 4 π ρ b c ¯ ( k ) = 0 , where c ¯ ( k ) is the 3D Fourier transform of c ( r ) [22], and ϕ ( r ) = sin ( k r ) / ( k r ) is the eigenfunction corresponding to zero-eigenvalue. A real value of k satisfying these conditions would exist when χ 1 = 1 4 π ρ b c ¯ ( 0 ) 0 [12]. The pertinent point to note is that the solution y ( r ) + α ϕ ( r ) violates the integral condition I = h ( r ) d r < . Thus singularity induced bifurcations generate an infinite number of solutions (for differenr α ), which violates the IC. These will occur for all values of ρ b for which χ 1 0 ; such ρ b values do arise in the spinodal region ( ρ s v ρ ρ s l ) on the isotherm (see Fig 2 ).

5. Discrete approximation to FIE

For developing discrete approximations, the integral over 0 r < in Eq.(17) is truncated at cutoff length R 1 , and a quadrature scheme (like Trapezoidal or Simpson’s rule) is applied to approximate the integral for each value of r. Then a collocation method is used to generate a matrix equation that can be solved for Δ y ( r ) at discrete set of quadrature points in 0 r R 1 . The kernel K ( r , r ) defined in Eq.(15) requires c ( r ) and y ( r ) in an extended interval 0 r 2 R 1 when the solution is sought in 0 r R 1 . This requirement arises whenever the unknown function, in integral equations in 3D, depends only on the radial co-ordinate r. Therefore, suitable asymptotic forms for c ( r ) and y ( r ) are needed in R 1 < r 2 R 1 . The ‘extreme’ asymptotic forms c ( r ) = 0 and y ( r ) = 0 in the last interval would necessitate the choice of a sufficiently large R 1 [13]. This problem, together with the efforts in computing the kernel K ( r , r ) and solving the FIE are, perhaps, the drawbacks of the real-space method. So the pioneering work [13] did not receive proper recognition.
A partial remedy to this problem, implemented in this paper, is to use accurate asymptotic forms c a s ( r ) and h a s ( r ) in R 1 < r 2 R 1 . Using the current solutions y ( r ) and c ( r ) , Eq.(5) may be used to compute these asymptotic forms whenever χ 1 > 0 . An alternate approach using Nyström’s interpolation formula, applicable to more general situations, is discussed at the end of this section. In developing a method for solving the FIE, first it is assumed that Δ y ( r ) = 0 for R 1 < r < . Note that Δ y ( r ) and s ( r ) vanish (within tolerance) in the entire half-line when Newton’s iterations converge. Then, Eq.(17) is reduced to a finite interval as
Δ y ( r ) = 4 π ρ 0 R 1 K ( r , r ) Δ y ( r ) r 2 d r + s ( r ) .
After determining Δ y , and thereafter c + Δ c and y + Δ y , the asymptotic forms to be used in the next Newton’s iteration are updated. Thus, although there is lag by one iteration, the asymptotic forms in R 1 < r 2 R 1 also converge if Newton’s iterations converge. The real-space method implemented in this paper, with an accurate method to solve the FIE in a finite interval [14,15], would provide solutions even when the IC is violated (see §2).
The Nyström’s method is commonly employed to transform a FIE to a matrix problem [14,15]. In the procedure adapted here, the upper limit R 1 is chosen such that R 1 = ( n × m ) δ , where δ = σ / n . Thus, there are a total of ( n × m ) major intervals and N = n × m + 1 solution-points in the discrete problem. The unknowns Δ y i = Δ y ( r i ) , where r i = ( i 1 ) δ for 1 i N , are the dependent variables. It is assumed that ( n × m ) is even and the integral in Eq.(24), for any value of r, can be approximated using piece wise quadratic polynomials. In the interval Δ j = r j 1 r r j + 1 (j even), any function f ( r ) is expressed as
f ( r ) = f j 1 L j L ( r ) + f j L j M ( r ) + f j + 1 L j R ( r ) , L j L ( r ) = ( r j r ) ( r j + 1 r ) / 2 δ 2 , L j M ( r ) = ( r r j 1 ) ( r j + 1 r ) / δ 2 , L j R ( r ) = ( r r j 1 ) ( r r j ) / 2 δ 2 .
Then, introducing certain weight factors W j ( r ) (which depends on r as explained below) to account for the volume element ( r 2 d r ) , Eq.(24) is approximated as
Δ y ( r ) = 4 π ρ j = 1 N W j ( r ) K ( r , r j ) Δ y j + s ( r ) .
For writing this equation, it is assumed that K ( r , r ) Δ y ( r ) is a quadratic polynomial in r in Δ j , as defined in Eq.(25). The collocation rule in Nyström’s method amounts to enforcing this equation at the discrete set of values { r i } , 1 i N , which reduces it to a matrix equation
( I M ) Δ y = s .
where Δ y and s are N dimensional vectors with components Δ y ( r j ) and s ( r j ) . Elements of the N × N dimensional matrix M are defined as
M 1 , j = 4 π ρ W j ( 2 ) K ( r 1 , r j ) , 1 j N , M i , 1 = 4 π ρ W 1 ( 2 ) K ( r i , r 1 ) , 2 i N , M i , j = 4 π ρ W j ( 1 ) K ( r i , r j ) , 2 i , j N .
where δ i , j denotes Kronecker delta. Two types of weight factors have been introduced here for reasons explained below. Note that the first row of W j ( r i ) (for which i = 1 ) is denoted as W j ( 2 ) for 1 j N . The first column of W j ( r i ) is just the transpose of its first row. The remaining elements are denoted as W j ( r i ) = W j ( 1 ) for 2 i , j N . Components of the source (residual) vector s are given by
s 1 = y ( r 1 ) + 4 π ρ j = 1 N W j ( 2 ) c ( r j ) { c ( r j ) + y ( r j ) } , s i = y ( r i ) + 4 π ρ W 1 ( 2 ) c ( r i ) { c ( r 1 ) + y ( r 1 ) } + + 4 π ρ j = 2 N W j ( 1 ) C ( r i , r j ) { c ( r j ) + y ( r j ) } , 2 i N .
Note that all the components s i would approach zero, if Newton’s method converges.
As the integrals C ( r 1 , r ) = C ( r , r 1 ) = c ( r ) , and similarly Y ( r 1 , r ) = Y ( r , r 1 ) = y ( r ) for r 1 = 0 , the two types of weight factors, W j ( k ) ( for k = 1 , 2 ) , introduced above are
W j ( k ) = r j 1 r j + 1 L j M ( r ) r k d r , j = 2 , 4 , , N 1 , W j + 1 ( k ) = r j r j + 1 L j R ( r ) r k d r + r j + 1 r j + 2 L j + 1 L ( r ) r k d r , j = 2 , 4 , , N 2 , W 1 ( k ) = r 1 r 2 L 2 L ( r ) r k d r , W N ( k ) = r N 1 r N L N R ( r ) r k d r .
These definitions account for the volume element in the integral operator exactly, and allow to retain Δ y 1 explicitly as an unknown in the discrete equations. Further, for evaluating the integrals C ( r , r ) and Y ( r , r ) , additional weight factors Q j ( k ) ( f o r k = 1 , 2 , 3 ) are introduced as
Q j 1 ( 1 ) = r j 1 r j L j L ( r ) r d r , Q j ( 1 ) = r j r j + 1 L j M ( r ) r d r , j = 2 , 4 , , 2 ( N 1 ) , Q j ( 2 ) = r j 1 r j L j M ( r ) r d r , Q j + 1 ( 2 ) = r j r j + 1 L j R ( r ) r d r , j = 2 , 4 , , 2 ( N 1 ) , Q j 1 ( 3 ) = r j 1 r j L j R ( r ) r d r , Q j ( 3 ) = r j r j + 1 L j L ( r ) r d r , j = 2 , 4 , , 2 ( N 1 ) .
Note that W j ( k ) δ k + 1 and Q j ( k ) δ 2 , where the proportionality factors are universal real numbers. Furthermore, the recursion relations defined as
V 1 c = 0 , V j c = V j 1 c + Q j 1 ( 1 ) c ( r j 1 ) + Q j ( 2 ) c ( r j ) + Q j 1 ( 3 ) c ( r j + 1 ) , V j + 1 c = V j c + Q j ( 1 ) c ( r j ) + Q j + 1 ( 2 ) c ( r j + 1 ) + Q j ( 3 ) c ( r j 1 ) , V 1 y = 0 , V j y = V j 1 y + Q j 1 ( 1 ) y ( r j 1 ) + Q j ( 2 ) y ( r j ) + Q j 1 ( 3 ) y ( r j + 1 ) , V j + 1 y = V j y + Q j ( 1 ) y ( r j ) + Q j + 1 ( 2 ) y ( r j + 1 ) + Q j ( 3 ) y ( r j 1 ) ,
for j = 2 , 4 , , 2 ( N 1 ) , allow the integrals C ( r i , r j ) and Y ( r i , r j ) to be expressed as
C ( r i , r j ) = 1 2 r i r j { V i + j 1 c V | i j + 1 | c } , Y ( r i , r j ) = 1 2 r i r j { V i + j 1 y V | i j + 1 | y } .
Here, i and j vary over 2 i , j N . Then the discrete approximation to the kernel ( r i , r j 0 ) reduces to
K ( r i , r j ) = C ( r i , r j ) + q ( r j ) { 2 C ( r i , r j ) + Y ( r i , r j ) } , 2 i , j N .
This method for computing the matrix elements, which relies on recursively defined one-dimensional arrays, is much faster in comparison to direct evaluation.
In Nyström’s method, Eq.(26) (known as Nyström’s interpolation formula) is useful to determine Δ y ( r ) , and thereafter Δ c ( r ) = q ( r ) Δ y ( r ) , for arbitrary values of r in the half-line [14]. In this paper, the interpolation formula is employed for the set { r i } , N + 1 i 2 N 1 , to update the asymptotic forms as y a s ( r i ) y a s ( r i ) + Δ y a s ( r i ) and c a s ( r i ) c a s ( r i ) + Δ c a s ( r i ) . This way of updating the asymptotic forms would work as far as Newton’s method provides a solution to Eq.(1), even in situations where the integral condition is violated.

5.1. Checking the algorithm

To check the accuracy of the discrete algorithm, consider the nonlinear integral equation [13]
ψ ( r ) = 4 π n 0 f ( r ) ψ ( r ) r 2 d r 1 2 r r | r r | r + r f ( t ) ψ ( t ) t d t .
It is easily verified, by direct substitution, that it has the exact solution ψ ( r ) = exp ( α r 2 / 2 ) when the function f ( r ) and n are chosen as f ( r ) = exp ( α r 2 / 2 ) and n = ( 2 α / π ) 3 / 2 . Alternatively, the second integral term can be written as f ( | r r | ) ψ ( | r r | ) (see Eq.(13 also), and the solution verified using Fourier transforms. This non-linear equation is similar to the OZE [13], and may be used to evaluate the dependence of numerical error on mesh size δ in discrete approximation. Application of Newton’s method, that follows exactly the same steps leading to Eq.(24), yields the linear equation
Δ ψ ( r ) = 4 π n 0 R 1 Φ ( r , r ) Δ ψ ( r ) r 2 d r + ω ( r ) ,
where the kernel and source terms are given by
Φ ( r , r ) = 2 f ( r ) 1 2 r r | r r | r + r f ( t ) ψ ( t ) t d t ,
ω ( r ) = ψ ( r ) + 4 π n 0 f ( r ) ψ ( r ) r 2 d r 1 2 r r | r r | r + r f ( t ) ψ ( t ) t d t .
Here ψ ( r ) is a guess solution. Its asymptotic form needed in the interval R 1 r 2 R 1 , for computing the integrals in these equations, is taken as the exact solution. The matrix equation obtained in the discrete version of Eq.(36) is expressed as
( I Φ ) Δ Ψ = Ω ,
where the matrix elements and source vector can be computed as done earlier. In fact the matrix elements Φ i , j are defined as in Eqs.(28) with K i , j replaced by Φ ( r i , r j ) . The latter is to be computed as done in the case of Eq.(34). New arrays as defined in Eq.(32) have to be generated in advance to obtain Φ ( r i , r j ) . Similarly, the components of source vector Ω are defined as in Eq.(29). Due to the use of exact asymptotic form of ψ ( r ) , the error in the discrete approximation arises from the use of quadratic functions only.
Results of computations (with cutoff length R 1 = 10 ) for the error e = f ψ (infinity-norm) versus mesh size δ are plotted (symbols) in Figure 1 for α = 1 and 0.5 . The broken lines are fits given by e = 0.25 δ 4 and e = 0.06 δ 4 for the two values of α . As expected, the global error reduction factor ( δ 4 ) of Simpson’s rule is realized in the algorithm.
Figure 1. Euclidean norm of error e in { ψ ( r i ) } versus mesh size δ are shown (symbols) for α = 1 and α = 0.5 . Broken lines are fits given by e = 0.25 δ 4 and e = 0.06 δ 4 for the two cases.
Figure 1. Euclidean norm of error e in { ψ ( r i ) } versus mesh size δ are shown (symbols) for α = 1 and α = 0.5 . Broken lines are fits given by e = 0.25 δ 4 and e = 0.06 δ 4 for the two cases.
Preprints 99693 g001

6. Continuation algorithms

The discrete set of equations and Newton’s method would provide a solution for given ρ and T, if a suitable starting solution is provided. More appropriately, solutions on an isotherm, starting from ρ = 0 to an upper limit ρ m a x are derived using the pseudo-arc length (PAL) continuation algorithm. This generates a sequence of solutions { y j , ρ j } , and hence a solution-curve, within the discrete approximation. The PAL algorithm, employing Fourier transforms to compute convolution integrals is already used to generate the ‘no-real-solution’ boundary lines for LJ and hard core Yukawa (HCY) potentials within the HNC closure [29]. This method is likely to fail on or inside the spinodal region on a sub-critical isotherm because the Fourier transforms of functions do not exist.
This section is devoted to discussing these methods for OZE and the MS closure (Eqs.(1) and (3)) in real-space. The discrete version of Eq.(20) is expressed as
( I M ) Δ y + Δ ρ F = 0 ,
where 0 is a N dimensional null vector. The elements of F are given by (see also Eq.(29))
F 1 = 4 π ρ j = 1 N W j ( 2 ) c ( r j ) { c ( r j ) + y ( r j ) } ; and   for 2 i N F i = 4 π ρ W 1 ( 2 ) c ( r i ) { c ( r 1 ) + y ( r 1 ) } + 4 π ρ j = 2 N W j ( 1 ) C ( r i , r j ) { c ( r j ) + y ( r j ) } .
In the limit Δ ρ 0 , Eq.(40) implies a system of ordinary differential equations of the form I M ( d y / d ρ ) + F = 0 . Then, Euler’s forward difference scheme implied in Eq.(40), or more accurate difference schemes are used to solve the differential equations, with an appropriate initial condition y 0 at ρ 0 = 0 . For a sufficiently small Δ ρ , solution of Eq.(40) will yield an approximation y a = y 0 + Δ y for the density ρ 1 = ρ 0 + Δ ρ . Thereafter, starting with y a as initial guess, Newton’s method is applied to obtain a new solution-pair { y 1 , ρ 1 } . In this way, a set of pairs { y j , ρ j } is generated as points defining the solution-curve. Generating a discrete set of solutions { y j , ρ j } via Eq.(40), and subsequent Newton’s method, is called the first order continuation method. Incidentally, zeroth order continuation amounts to repeated use of Newton’s method with solution from previous value of ρ as the initial guess. First order continuation belongs to the class of ‘predictor-corrector’ methods: solution to Eq.(40) provides the predictor while Newton’s method works as a corrector.

6.1. PAL algorithm

In the neighborhoods and at values of ρ where ( I M ) is singular (called critical points), first order (or zeroth order) method will fail and different continuation algorithms are needed. In these algorithms, the pairs { y j , ρ j } are considered as points on a curve in N + 1 dimensional space, which is labeled by the arc length (see §4). Thus, = 0 corresponds to the known initial point { y 0 , ρ 0 } . Assuming that solution-pairs up to j are already obtained, Eq.(40) is modified as (see Eq.(21))
( I M j ) y ˙ j + ρ ˙ j F j = 0 .
Here y ˙ j and ρ ˙ j denote the components of the tangent vector to solution-curve at j . However, now there are N + 1 unknowns in Eq.(42), but the deficiency is removed by normalizing the tangent vector to unit length as | | y ˙ j | | 2 + ρ ˙ j 2 = 1 , where | | y ˙ j | | denotes the Euclidean norm of y ˙ j . A simple way to impose this normalization is to append Eq.(42) with an additional equation, y ˙ j k = 1 , which assigns value unity to the k t h component of y ˙ j . Thus Eq.(42) is modified as
I M j F j e k t 0 y ˙ j ρ ˙ j = 0 1 ,
where e k is a unit vector with k t h component unity, and superscript t denotes transpose of the column vector. The value of y ˙ j k is finally recomputed by imposing the normalization condition of the tangent vector. The choice of k is somewhat arbitrary; one choice would be such that | y ˙ j 1 k | is largest of all components. This method will not fail even at a fold bifurcation point, as the enlarged matrix in Eq.(43) is non-singular (see §6.2).
To obtain a proper predictor in PAL algorithm, the arc length equation defined as | | Δ y | | 2 + ( Δ ρ ) 2 = Δ 2 is added to Eq.(40) to close the system of N + 1 equations. Here | | Δ y | | is the Euclidean norm of Δ y . Being nonlinear, this relation is not useful directly. So a linear version of the same (hence the name pseudo arc length), expressed as y ˙ t · Δ y + ρ ˙ Δ ρ = Δ , is used in PAL algorithm. This relation follows by replacing ( Δ y / Δ ) and ( Δ ρ / Δ ) with the corresponding derivatives, which is well justified if Δ is sufficiently small. Using their latest values computed via Eq.(43), the system of equations to be solved in PAL algorithm, for predicting the solution from j , is
I M j F j y ˙ j t ρ ˙ j Δ y Δ ρ = 0 Δ
This solution readily provides the predictor to next solution pair as y j + 1 p = y j + Δ y and ρ j + 1 p = ρ j + Δ ρ . Newton’s method is applied as a corrector, starting with y j + 1 p and ρ j + 1 p , to this extended system of N + 1 equation. The corrector iterations are performed by solving the equation
I M j + 1 F j + 1 y ˙ j + 1 t ρ ˙ j + 1 Δ y Δ ρ = s j + 1 Δ
Note that s j + 1 is now present in the source vector as the predicted (as well as subsequent) solutions do not lie on the solution-curve. All the matrix elements in M j + 1 , tangent components y ˙ j + 1 t , ρ ˙ j + 1 and residual vector s j + 1 in Eq.(45) are to be updated, as iterations progress, using the most recent solutions y j + 1 y j + 1 + Δ y and ρ j + 1 ρ j + 1 + Δ ρ at arc length j + 1 . In a more commonly used version of the PAL algorithm, an additional condition that the correction vector ( Δ y , Δ ρ ) is orthogonal to the tangent vector y ˙ j + 1 , ρ ˙ j + 1 is employed. Correction iterations converge faster as these proceed in orthogonal direction to the tangent vector. This makes Δ = 0 in the source vector and Newton’s corrector equation reduces to
I M j + 1 F j + 1 y ˙ j + 1 t ρ ˙ j + 1 Δ y Δ ρ = s j + 1 0

6.2. PAL algorithm and singularities

PAL algorithm outlined above will work in situations when ( I M ) is singular for a particular ρ , but has a simple rank deficiency. Thus its rank is N 1 at the singularity, i.e., one less than the dimension N . In that case its zero-eigenvalue is simple, which means that the null space is one-dimensional. Now, the rank of the augmented matrix [ ( I M ) | F ] can be either N or N 1 . In the first case, ρ is called a fold bifurcation point, where as in the second case it denotes a simple bifurcation point corresponding to either trans-critical or pitchfork bifurcations [18]. For the case of fold bifurcation point (case-1), as the rank of augmented matrix is N , there is only one vector ( y ˙ , ρ ˙ ) that satisfy the equation
( I M ) y ˙ + ρ ˙ F = 0 ,
which is same as the tangent equation. Now, the augmented matrices in Eq.(43) or (44) can not have a zero-eigenvalue, as the same vector ( y ˙ , ρ ˙ ) can not also satisfy the eigen-value equation. In this case the PAL algorithm allows one to move on the solution-curve smoothly. Determinant of ( I M ) will change sign at the fold bifurcation point, and can be monitored to locate the same [18]. In fact, zeros of linear interpolation polynomials of solutions, constructed with successive positive and negative values of determinant, provide good approximations to the fold bifurcation point [18]. The singularity of ( I M ) at the fold bifurcation point, say at ρ f , implies a consistency condition (for a non-trivial y ˙ to exist) that ( ρ ˙ f F ) has to be orthogonal to its right eigenvector - same as the eigenvector of ( I M ) t . However, as the rank of the augmented matrix is N , the vector F can not be expressed as a linear combination of the columns of ( I M ) , i.e., there does not exit a non-trivial vector c defined as F = ( I M ) c . So the consistency condition is satisfied only when the slope ρ ˙ f = 0 . Using Taylor’s expansion, it is shown that a negative value of second derivative ρ ¨ f , with respect to , leads to a fold turning to left at ρ f , while a positive value yields a right-turning fold (see Figure 2).
Because the rank of augmented matrix (which has N + 1 columns ) is N 1 for the case of simple bifurcations (case-2), there are two independent vectors ( u ˙ k , ρ ˙ k ) , for k = 1 , 2 , that satisfy the equation ( I M ) u ˙ k + ρ ˙ k F = 0 . In fact, these two tangent vectors can be determined [17,18], thereby making it possible to trace both branches passing through the bifurcation point. Applications of these ideas are discussed below (see §7.2, §8.2).
The methods described above is now applied to LJ potential with MS closure. Results obtained for the isotherm corresponding to reduced temperature T * = k B T / ϵ = 1.15 are shown in Figure 2. Mesh size δ = 0.05 σ and cutoff length R 1 = 20 σ leading to N = 401 are used in the discrete approximations. Proper asymptotic forms for y ( r ) and c ( r ) , as explained in §5, also are employed. Inverse compressibility χ 1 = β ( P / ρ ) and determinant Δ I M of the Jacobian matrix are plotted versus reduced density ρ * = ρ σ 3 . This isotherm is well in the sub-critical region as the critical temperature is T c * = 1.316 . Fold bifurcation points, corresponding to Δ I M = 0 , occur at ρ * = 0.1339 and 0.3797 on the vapor and liquid sides, respectively. Similarly, spinodal points at which χ 1 = 0 , are found at ρ * = 0.1327 and 0.4584 , on the two sides. The difference between the two types of densities is very small on the vapor side, although the fold bifurcation point is larger than the spinodal point. The solution-branch for ρ * < 0.1339 is nonphysical because χ 1 increases with density. So it may be concluded that the spinodal at ρ * = 0.1327 is not on a physical branch. The fold bifurcation and spinodal are well separated on the liquid side, and χ 1 0 between them. While there is a unique solution branch for densities higher than the spinodal (on the liquid side), two branches are found for densities lower than the same on vapor side. These results are quite similar to those found via real-space approach for HNC closure [13]. However, computations employing Fourier transforms did not find density-regions where χ 1 is negative [28]. Similar results are found for other isotherms corresponding to T * = 1.25 , 1 and 0.9 .
Figure 2. Inverse compressibility χ 1 = β ( P / ρ ) versus ρ * = ρ σ 3 along isotherm for T * = k B T / ϵ = 1.15 , in the vapor side ( curve 1) and liquid side (curve 2) for LJ potential. Determinant Δ I M of the Jacobian matrix on vapor and liquid sides are shown in curves 3 and 4, respectively. While χ 1 vanish at ρ * = 0.1327 and 0.4584 , fold bifurcation points (where determinant vanish) are at ρ * = 0.1339 and 0.3797 (see text in §6.2).
Figure 2. Inverse compressibility χ 1 = β ( P / ρ ) versus ρ * = ρ σ 3 along isotherm for T * = k B T / ϵ = 1.15 , in the vapor side ( curve 1) and liquid side (curve 2) for LJ potential. Determinant Δ I M of the Jacobian matrix on vapor and liquid sides are shown in curves 3 and 4, respectively. While χ 1 vanish at ρ * = 0.1327 and 0.4584 , fold bifurcation points (where determinant vanish) are at ρ * = 0.1339 and 0.3797 (see text in §6.2).
Preprints 99693 g002

7. Generalization for complex solutions

Complex solutions to OZE and PY closure, inside the two-phase region, were originally found for hard spheres with stickiness [26], and LJ potential [11]. While these investigations employed Baxter’s version of OZE, similar results were also found by generalizing the Fourier transform method to the complex domain [34]. More recently, complex solutions are found for LJ potential with HNC [35] as well as several other closures [28]. All these results may be thought of as bifurcation from real to complex branches of solutions [19]. Most interestingly, bifurcation analysis to complex domain connects different branches of real solutions, which appear as disconnected branches if restricted to real valued functions [19]. In this section, the real-space method to OZE and MS closure is generalized to the complex domain.
If c ( r ) and y ( r ) are generalized as c ( r ) = c R ( r ) + ı c I ( r ) and y ( r ) = y R ( r ) + ı y I ( r ) , where the subscripts R and I denote the real and imaginary parts, respectively, then Eq.(8) will yield separate equations for the real and imaginary components. Alternatively, the linear Fredholm integral equation in Eq.(24), can be generalized as
Δ y R ( r ) = 4 π ρ 0 R 1 [ K R ( r , r ) Δ y R ( r ) K I ( r , r ) Δ y I ( r ) ] r 2 d r + s R ( r ) , Δ y I ( r ) = 4 π ρ 0 R 1 [ K I ( r , r ) Δ y R ( r ) + K R ( r , r ) Δ y I ( r ) ] r 2 d r + s I ( r ) ,
where the kernels K R ( r , r ) and K I ( r , r ) are given by
K R ( r , r ) = C R ( r , r ) + q R ( r ) P R ( r , r ) q I ( r ) P I ( r , r ) , K I ( r , r ) = C I ( r , r ) + q I ( r ) P R ( r , r ) + q R ( r ) P I ( r , r ) .
For simplifying equations, P Z ( r , r ) = { 2 C Z ( r , r ) + Y Z ( r , r ) } is introduced, where, C Z ( r , r ) and Y Z ( r , r ) (subscript Z denotes either R or I) are redefinition of Eqs.(13) and (14) employing c Z ( r ) and y Z ( r ) . Therefore, they can be computed as in the case of real solutions. Similarly, q Z ( r ) = [ h Z ( r ) + g Z ( r ) B R ( r ) ] provides the real and imaginary components, where it is assumed that the bridge function is defined with the real part y R ( r ) . Asymptotic forms of these functions in R 1 < r 2 R 1 are needed in computing the kernels given above. Residual terms s R ( r ) and s I ( r ) are found to be
s R ( r ) = y R ( r ) + 4 π ρ 0 R 1 [ C R ( r , r ) p R ( r ) C I ( r , r ) p I ( r ) ] r 2 d r , s I ( r ) = y I ( r ) + 4 π ρ 0 R 1 [ C R ( r , r ) p I ( r ) + C I ( r , r ) p R ( r ) ] r 2 d r ,
where p Z ( r ) = { c Z ( r ) + y Z ( r ) } . Discrete form of Eq.(48) is obtained following the same procedures and formulas used for the case of real solutions. There are no additional complications, although the dimension of the matrix equation to be solved is doubled to 2 N . Nevertheless, some of the formulas are listed in next section for completeness.

7.1. Discrete form for complex solutions

For complex solutions, it is straightforward to generalize Nyström’s interpolation formulas in Eq.(26) as
Δ y R ( r ) = 4 π ρ j = 1 N W j ( r ) [ K R ( r , r j ) Δ y R j K I ( r , r j ) Δ y I j ] + s R ( r ) , Δ y I ( r ) = 4 π ρ j = 1 N W j ( r ) [ K I ( r , r j ) Δ y R j + K R ( r , r j ) Δ y I j ] + s I ( r ) .
Application of the collocation rule readily yields the matrix equation
I M R M I M I I M R Δ y R Δ y I = s R s I .
Elements of the matrices M Z are defined exactly as in Eq.( 28) with K Z ( r i , r j ) (where Z denotes either R or I) in place of K ( r i , r j ) . Computation of matrix elements also proceeds exactly as in the case of real solutions, however, additional arrays, like in { V j c } and { V j y } ( see Eq.( 32) ), need to be generated with real and imaginary components { c Z ( r j ) } and { y Z ( r j ) } . The known vectors s R and s I are also computed using expressions analogous to those in Eq.(29).
Nyström’s interpolation formulas in Eq.(51) are useful to determine Δ y Z a s ( r ) , for arbitrary values of r along the half-line [ 0 , ) , as in the case of real solutions. Thereafter Δ c Z a s ( r ) are readily computed using Δ c Z a s ( r ) = q Z ( r ) Δ y Z a s ( r ) . These are employed for the set of ordinates { r i } in the range N + 1 i 2 N 1 to update the asymptotic forms to y Z a s ( r i ) y Z a s ( r i ) + Δ y Z a s ( r i ) and c Z a s ( r i ) c Z a s ( r i ) + Δ c Z a s ( r i ) . Continuation of solutions for different values of ρ is obtained using the same PAL algorithm, generalizing Eq.(43) to get tangent vector, Eq.(44) for predictor and Eq.(46) for corrector iterations. Of course, the starting point is ρ = 0 and it is necessary to employ 2 N dimensional matrices throughout.

7.2. Complex bifurcation

The method outlined above, is along the well known procedure to look for complex solutions to non-linear equations [19]. The imaginary components of the solutions are automatically generated, as discussed below, on continuing the solution along the arc length - it is found unnecessary to tweak the OZE by adding an imaginary term as done in earlier works [28,35]. For the complex solution, the tangent equation in Eq.(47) still applies, but the matrix M , and vectors y ˙ and F are taken as complex objects (with subscripts R and I denoting real and imaginary parts, respectively). This equation together with its derivative with respect to , are written as
( I M ) y ˙ + ρ ˙ F = 0 , ( I M ) y ¨ M y y ˙ y ˙ ( M ρ F y ) y ˙ ρ ˙ + ρ ¨ F = 0 .
Here M y , F y and M ρ denote derivatives with components of y and ρ , respectively. Note that M y is a complex matrix of dimension ( N × N 2 ) and y ˙ y ˙ , which denotes a direct product, is a N 2 dimensional vector. Expressions for these are easily derived using Eqs.(28) and (41), after generalizing them to complex solutions. As explained earlier (with reference to Eq.(19), the explicit dependence of c ( r ) , and hence of F , on ρ is omitted, which explains the absence of the F ρ ρ ˙ ρ ˙ term. Also, y ¨ and ρ ¨ are second derivatives with arc length . At the fold bifurcation point ρ f , the identities M I = M y I = F I = 0 and ρ ˙ f = 0 hold, although y ˙ and y ¨ can still be complex vectors. So at density ρ f , Eq.(53) reduces to
(54) ( I M R ) y ˙ = 0 , (55) ( I M R ) y ¨ M y R y ˙ y ˙ + ρ ¨ f F R = 0 ,
Now, Eq.(54) shows that y ˙ = θ f where f is the (real) eigenvector corresponding to zero-eigenvalue of ( I M R ) and θ an arbitrary (complex) constant. The corresponding eigenvector of ( I M R ) t is denoted by g . On substituting for y ˙ , and pre-multiplying with g t , Eq.(55) yields
θ 2 g t · M y R f f ρ ¨ f g t · F R = 0 .
This equation consists of two unknowns, however the possibility of choosing multiple values of θ (consistent with the equation) would indicate the emergence of different solution-branches at ρ f . Now, θ has to be either purely real ( θ 1 ) or purely imaginary ( ± ı θ 2 ) . The former case yields ρ ¨ f = θ 1 2 ( A / B ) , where A = [ g t · M y R f f ] and B = ( g t · F R ) . This provides the real-branch at ρ f found earlier in §6.2. The latter case gives ρ ¨ f = θ 2 2 ( A / B ) , which shows the emergence of imaginary branch at ρ f .
Exactly similar analysis is done for the spinodal point starting with the discrete version of the tangent equation in Eq.(22). Its complex version and derivative with are expressible exactly as in Eq.(53), although with minor differences. Now, M denotes the discrete version of operator L 1 , and F also contains an additional term arising out of the operator L 2 in Eq.(22). Again, all the matrices and vectors (except y ˙ ) are real at the spinodal density ρ s , however, ρ ˙ s 0 . Thus, the equations to be analyzed at ρ s are given by
(57) ( I M R ) y ˙ + ρ ˙ s F R = 0 , (58) ( I M R ) y ¨ M y R y ˙ y ˙ ( M R ρ F y R ) y ˙ ρ ˙ s + ρ ¨ s F R = 0 .
Note the presence of extra terms, in comparison to those in Eqs.(54) and (55), because ρ ˙ s 0 . Now, the crucial point is that ( g t · F R ) = 0 , because the solution y ˙ exists, although ( I M R ) is singular at ρ s . Then, the general solution of Eq.(57) is expressed as y ˙ = θ f + v , where θ is a (complex) constant and v is a particular solution. This solution is chosen such that ( g t · v ) = 0 , as any component of v along f is absorbed in θ f . Now, the steps leading to Eq.(56) are repeated to obtain the quadratic equation: A θ 2 + ( B 1 ρ ˙ s B 2 ) θ + ( C 1 ρ ˙ s C 2 ) = 0 . Here A is already defined, B 1 = [ g t · M y R ( f v + v f ) ] and B 2 = [ g t · ( M R ρ F y R ) f ] . The remaining constants are obtained by replacing f by v in A and B 2 , respectively. That is, C 1 = [ g t · M y R v v ] and C 2 = [ g t · ( M R ρ F y R ) v ] . The quadratic equation in θ has two real roots or a pair of complex conjugate roots, depending on the values of ρ ˙ s and other parameters. In addition, it yields ρ ˙ s = ( θ 2 A + B 1 + C 1 ) / ( θ B 2 + C 2 ) . The two real values of θ give ( y ˙ , ρ ˙ s ) on the real solution-branches while the complex conjugate pair generates the complex solution-branches.
Thus, two sets of complex solutions may occur in the two-phase region; one originating at the fold bifurcation point and the other at the spinodal point. The second may be absent in some specific cases, as in the case of Baxter’s model [26].
The algorithm described above is applied to LJ potential. Mesh size δ = 0.05 σ and cutoff length R 1 = 12 σ , which yields N = 241 , are used. Asymptotic forms for real and imaginary parts of y ( r ) and c ( r ) are also used, as explained earlier in §7.1. Computations started at ρ = 0 and the complex solutions are generated automatically by the algorithm. Complex solutions y ( r ) and c ( r ) lead to complex values for thermodynamic quantities. Values of inverse compressibility χ 1 = β ( P / ρ ) against ρ * = ρ σ 3 , along isotherm for T * = k B T / ϵ = 1.15 , are shown in Figure 3. Imaginary part (red curve1) is identically zero outside the spinodal region, while the real part (blue curve 2 ) connects the spinodal points on vapor and liquid sides. Symbols, marked as 3 and 4 in the figure, show χ 1 versus ρ * from purely real solution branches discussed in Figure 2. Importantly, the complex solutions are seen to cover the entire spinodal region, as clearly evident on the liquid side. It is possible that there is another complex branch connecting the fold bifurcation points, however, this was not found in the results. Special procedures to trace a complex branch at the fold on liquid side may be necessary for this. Comparison of χ 1 (on real branches) in these figures also shows that all its variations are retained either with R 1 = 12 or R 1 = 20 .
Figure 3. Complex inverse compressibility χ 1 = β P / ρ versus ρ * = ρ σ 3 along isotherm for T * = k B T / ϵ = 1.15 . Real part ( red curve2) and imaginary part (blue curve1) connect spinodal points on vapor and liquid sides. Symbols (marked as 3 and 4 ) show χ 1 vs ρ * from purely real solution branches (see Figure 2). Note that the complex solutions cover the entire spinodal region (see text in §7.2).
Figure 3. Complex inverse compressibility χ 1 = β P / ρ versus ρ * = ρ σ 3 along isotherm for T * = k B T / ϵ = 1.15 . Real part ( red curve2) and imaginary part (blue curve1) connect spinodal points on vapor and liquid sides. Symbols (marked as 3 and 4 ) show χ 1 vs ρ * from purely real solution branches (see Figure 2). Note that the complex solutions cover the entire spinodal region (see text in §7.2).
Preprints 99693 g003

8. Coupling parameter expansion

The coupling parameter expansion (CPE) starts with splitting the potential as U ( r , λ ) = U 0 ( r ) + λ U a ( r ) , where U 0 ( r ) and U a ( r ) are, respectively, purely repulsive and attractive components, as per WCA prescription [25]. The coupling parameter λ tunes the strength of the attractive potential; λ = 0 corresponds to the reference part while λ = 1 provides the full potential. Now, all correlation functions also are considered as functions of λ . In CPE, these functions are assumed to be analytic around λ = 0 , and expanded as Taylor’s series in λ about the respective functions corresponding to their reference parts [21]. Thus, the function y ( r , λ ) is expanded as
y ( r , λ ) = n = 0 1 n ! λ n y n ( r ) .
Here y n ( r ) denotes the n t h derivative ( 0 n < ) of y ( r , λ ) at λ = 0 ; y 0 ( r ) being the total correlation function for reference potential. Other correlation functions, c ( r , λ ) and h ( r , λ ) , are also represented via such series expansions. The main assumption in CPE is that these expansions converge over the entire range 0 λ 1 , and the contributions from the attractive part of potential are adequately accounted via the derivatives like y n ( r ) . Truncation of the series at zeroth order ( n = 0 ) yields the first order perturbation theory [1]. But the general formulation of CPE can incorporate quite higher order derivatives, however, by taking recourse to the OZE [22].
Using the binomial expansion of derivatives of product of two functions, the n t h derive of Eq.(8) is readily found to be
y n ( r ) = ρ c n ( r ) [ c 0 ( | r r | ) + y 0 ( | r r | ) ] + + c 0 ( | r r | ) [ c n ( r ) + y n ( r ) ] d r + q n ( r ) , q n ( r ) = ρ m = 1 n 1 C n m c n m ( | r r | ) [ c m ( r ) + y m ( r ) ] d r , n 1 ,
where { C n m } are the binomial coefficients. Similarly, the closure relation in Eq.(2) is differentiated n times to obtain
c n ( r ) = [ h 0 ( r ) + g 0 ( r ) B 0 ( r ) ] y n ( r ) + w n ( r ) , n 1 w n ( r ) = g 0 β U A δ n 1 + g 0 B n * + + m = 1 n 1 C m n 1 [ β U A δ n m , 1 + y n m + B n m ] [ y m + c m ] , n 1 ,
where B 0 ( r ) is derivative of B [ y ] computed at y 0 (see Eq.(9)). The derivative of bridge function, with respect to λ , are expressed as B n = B 0 y n + B n * ( n 1 ) , so that B n * contains only derivatives of y ( r ) of order less than n. The derivatives B n * for the MS closure are obtained via recursion relations [21,22]
B n * = B 0 ρ * β U a δ n 1 ( 1 + B 0 ) m = 0 n 2 C m n 1 O m + 1 O n 1 m , n 1 , O n = ( 1 + B 0 ) y n + B n * , n 1 .
The summation term contributes only for n 2 ; so this recursion equation provides B n * explicitly at all orders. Now, substitution of the above expression for c n ( r ) in Eq.(60) gives the linear equation
y n ( r ) = ρ [ c 0 ( | r r | ) + q 0 ( r ) { 2 c 0 ( | r r | ) + y 0 ( | r r | ) } ] y n ( r ) d r + + s n ( r ) , n 1 , s n ( r ) = q n ( r ) + ρ [ 2 c 0 ( | r r | ) + y 0 ( | r r | ) ] w n ( r ) d r , n 1 ,
where q 0 ( r ) = [ h 0 ( r ) + g 0 ( r ) B 0 ( r ) ] . Note that the kernel of the above FIE is similar to that in Eq.(12), and so can be re-expressed as in Eq.( 15). Computation of correlation functions y 0 ( r ) and c 0 ( r ) for the reference system (on any isotherm) is readily done using Newton’s method and first order continuation. No singularities are encountered for the reference system as it is purely repulsive. The FIE is easily converted into a matrix equation like ( I M 0 ) y n = s n using Nyström’s method elaborated in §5. The source vector s n contains only derivatives of order less than n, and the matrix to be inverted is same at all orders.
The derivatives y n ( r ) and c n ( r ) are short ranged functions, just like y 0 ( r ) and c 0 ( r ) , and so their Fourier transforms exist even in the coexistence region. Therefore, it is not necessary to continue in real-space; the convolution integrals in Eq.( 63) are easily computed in Fourier space. Taking Fourier transforms yields
y ¯ n ( k ) = ρ 2 c ¯ 0 ( k ) + y ¯ 0 ( k ) 1 ρ c ¯ 0 ( k ) [ q 0 ( r ) y n ( r ) ¯ + w ¯ n ( k ) ] + q ¯ n ( k ) 1 ρ c ¯ 0 ( k ) .
This equation is identical to that obtained earlier [21,22]. Further, it is quite similar to Eq.(18), which represents Newton’s method in Fourier space. However, the factor ( 1 ρ c ¯ 0 ) in the denominators will never become zero as it refers to the reference system. The operator involved in Eq.(64) is the Fourier transform integral, and so this equation is easily solved via Krylov space methods [21]. For applications, λ is set to 1 in Eq.(59), and the series is terminated at a suitable order. All results obtained below use 7th order expansion [22], thus keeps terms up to n 6 .
Two independent methods for computing thermodynamic quantities, one via the free energy route and the other using the compressibility route, were used to obtain the phase diagram of LJ fluid [22]. Both approaches agree very well each other, and with results derived from simulation data [36]. The former method does not need data in the spinodal region for computing the phase diagram; it imposes equality of pressure and chemical potential of the gaseous and liquid phases. But the latter uses integration of inverse compressibility ( χ 1 = β P / ρ ) along ρ , starting from ρ = 0 . Because it is computed as 1 ρ c ( r ) d r , negative values of χ 1 - signaling slow decay of c ( r ) - are obtained along any isotherm in the sub-critical region, which is discussed next.

8.1. Negative compressibility

Convergence of the CPE was analyzed in detail earlier [22]. There it is shown that the series truncated at 7th order, which retain terms up to n 6 , converge to well defined functions in the super-critical as well as sub-critical regions for LJ fluid. However, the integral constraint I = h ( r ) d r < is violated in regions where χ 1 is negative (see §2). Suppose y e x ( r ) and c e x ( r ) are solutions (presumably converged) obtained via the CPE. Then, the starting Eq.(8) shows that y α ( r ) = y e x ( r ) + α ϕ ( r ) is also a solution where ϕ ( r ) satisfies the homogeneous equation ϕ ( r ) ρ c e x ( | r r | ) ϕ ( r ) d r = 0 , and α an arbitrary constant. In fact, ϕ ( r ) is the eigenvector corresponding to zero-eigenvalue of the integral operator in OZE, with a given function c ( r ) . It was argued in §4 that such functions do exist whenever χ 1 is negative. Equally important, the spatial integrals of ϕ ( r ) and y α ( r ) do not exist. Thus CPE gives rise to an infinite number of solutions, which violate the integral constraint in the spinodal region.
It has been argued [13] that Monte Carlo simulations [37] and mean field theories of fluids [24] can be constructed within a restricted ensemble approach. Here the total volume occupied by the system is divided into a large number of cells, and the particles in each cell is restricted to move within that cell. This restricts the range of allowed statistical (number) fluctuations within the system, thereby making it possible to connect vapor and liquid branches (through the spinodal region) on a sub-critical isotherm. Recent simulations using the transition-matrix Monte Carlo (TMMC) clearly show the emergence of van der Waals loops in pressure versus volume curves in the spinodal region [23]. So it seems appropriate to relate the results of CPE to those of restricted ensemble approach. However, CPE does not fix the range of allowed fluctuations arbitrarily, as is the practice in simulations. The correlation functions of the reference system automatically fix the range of fluctuations - like in other mean field theories. Solutions to OZE incorporate fluctuations at all ranges, although approximately [13]. The fold bifurcation points (on vapor and liquid sides) and ‘no-real-solution regions’ between them, obtained via Newton’s method (see §6), seem to be a consequence of these very long range fluctuations. In fact, real-space method applied to OZE with lower cutoff length ( R 1 = 6 σ ) , although within the PY closure [38], do not show fold bifurcation points; it yields connected isotherms (with van der Waals loops) between vapor and liquid branches.
Figure 4. Comparison of β P σ 3 versus ρ * = ρ σ 3 obtained using CPE (lines) and TMMC simulations (symbols) on different isotherms, labeled with T * = k B T / ϵ . The blue and red lines are obtained via the free energy and compressibility routes, respectively. Both set of results show the van der Waals loops in the spinodal region (see text in §8.1).
Figure 4. Comparison of β P σ 3 versus ρ * = ρ σ 3 obtained using CPE (lines) and TMMC simulations (symbols) on different isotherms, labeled with T * = k B T / ϵ . The blue and red lines are obtained via the free energy and compressibility routes, respectively. Both set of results show the van der Waals loops in the spinodal region (see text in §8.1).
Preprints 99693 g004
To check the accuracy of CPE, reduced pressure β P σ 3 versus reduced density ρ * = ρ σ 3 is compared in Figure 4, for different isotherms labeled with T * = k B T / ϵ , with TMMC simulation data for T * = 0.8 , 1.0 , 1.2 , 1.35 , 1.5 . As critical temperature ( T c * = ( 1.326 ) obtained in CPE [22] compare well with that from simulations ( T c * = 1.316 ) [36], it may be concluded that CPE provides excellent results in the super-critical region ( see isotherms for 1.35 , 1.5 ). In the spinodal region, the agreement is quite good, except for T * = 0.8 . Pressure computed via free energy route (blue curves) agree better with data. The differences between CPE results obtained via free energy (blue curves) and compressibility (red curves) routes are indicative of the lack of exact thermodynamic consistency of the MS closure. These results support the view that CPE provides accurate results within the restricted ensemble approach.

8.2. Non-uniform singlet density

Finally, there is an issue within CPE because χ 1 computed via c ( r ) and g ( r ) will not match in the spinodal region, because the former is negative. As argued below, this is due to the inherent drawback of using a spatially uniform density in that region. It is shown next that the singlet density ρ ( r ) bifurcates to a non-uniform distribution when χ 1 = 0 . A starting point to show the emergence of this bifurcation is the equation for non-uniform singlet density ρ ( r ) in a fluid [39]
ln [ 1 + ψ ( r ) ] = ρ c ( | r r | ) ψ ( r ) d r ,
where ψ ( r ) = [ ρ ( r ) ρ ] / ρ , and c ( r ) is the direct correlation function obtained within CPE or any other mean field theories. Being the correlation function in an isotropic fluid, c ( r ) depends only on the radial coordinate r. Here ρ denotes the mean density used in OZE, CPE, mean field theories, etc. Clearly, uniform density, corresponding to the solution ψ ( r ) = 0 of Eq.(65), exits for all ρ . What is to be checked is the emergence of bifurcation to a solution ψ ( r ) 0 at a spinodal point (see §7.2). For this analysis, ψ ( r ) and ρ are considered as functions of the arc length parameter , used for labeling a solution-curve. The slopes ψ ˙ and ρ ˙ , with respect to arc length , along the solution-curve, satisfy the equation
1 1 + ψ ψ ˙ ρ L ψ ˙ ρ ˙ L + L ρ ψ = 0 ,
where the operator L is defined as L ψ = c ( | r r | ) ψ ( r ) d r . The dependence of c ( r ) on ρ is incorporated in L ρ . Differentiating once more with yields
1 1 + ψ ψ ¨ ρ L ψ ¨ ρ ¨ L + L ρ ψ 1 ( 1 + ψ ) 2 ψ ˙ ψ ˙ 2 ρ ˙ L + L ρ ψ ˙ ρ ˙ L ρ + L ρ ρ ψ = 0 ,
where L ρ ρ accounts for second derivative of c ( r ) with ρ . Suppose ρ b is a density where bifurcation occurs from ψ ( r ) = 0 solution. At that point ( ρ b , 0 ) on the solution-curve, these equations reduce to
ψ ˙ ρ b L ψ ˙ = 0 , I ρ b L ψ ¨ ψ ˙ ψ ˙ + 2 ρ b ˙ L + L ρ b ψ ˙ = 0 .
The first equation shows that ψ ˙ ( r ) = α ϕ ( r ) is a solution where ϕ ( r ) is the eigenvector of ( I ρ b L ) corresponding to zero-eigenvalue, and α an arbitrary constant. By direct substitution it is verified that ϕ ( r ) = exp ( ι k b · r ) is a non-trivial eigenvector, provided 1 ρ b c ¯ ( k b ) = 0 , where c ¯ ( k b ) is the 3D Fourier transform of c ( r ) , and k b the (radial) Fourier variable. For a non-trivial solution ψ ¨ in the second equation, the terms in the square bract must be orthogonal to ϕ ( r ) as ( I ρ b L ) is singular. This condition provides a quadratic equitation for α , viz.
α 2 δ k b , 0 + 2 α ρ b ˙ c ¯ ( k b ) + c ¯ ρ b ( k b ) = 0 ,
where δ k b , 0 denotes Kronecker delta, and c ¯ ρ b ( k b ) is the density-derivative of c ¯ ( k b ) at ρ b . This equation shows that α = 0 is always a solution, which corresponds to uniform density distribution ψ ( r ) . However, there is a non-trivial solution when ρ b = ρ s and k b = 0 , where ρ s is the density corresponding to a spinodal point. This is given by α = 2 ρ s ˙ c ¯ ( 0 ) + c ¯ ρ s ( 0 ) . Note that c ¯ ( 0 ) = ρ s 1 and c ¯ ρ s ( 0 ) = ρ s 2 ρ s 1 χ ρ s 1 where the subscript ρ s denotes density-derivative at ρ s . Thus the singlet density bifurcates to a non-uniform distribution at the spinodal points. Although the linear part ( I ρ b L ) is singular for k b 0 , i.e. inside the spinodal region, the bifurcation occurs only at the spinodal points.
For illustrative purposes, consider ρ close to ρ s . Then ψ ( r ) < < 1 , and so ln [ 1 + ψ ( r ) ] in Eq.(65) may be replaced with ψ ( r ) . Further, assuming that ψ ( r ) varies slowly with r , so that ψ ( r ) is approximated by Taylor’s series, Eq.(65) reduces to
3 ρ C 2 2 ψ [ 1 ρ C 0 ] ψ ( r ) = 0 ,
where the constants C 0 and C 2 are the same as in Eq.(5). When χ 1 = [ 1 ρ C 0 ] is negative, i.e., for ρ within the spinodal region, Eq.(70), together with (say) periodic boundary conditions, will yield oscillatory solutions. Numerical treatment of Eq.(65) will be necessary to establish the occurence of non-trivial ψ ( r ) in the entire spinodal region.

9. Summary

One of the objectives in this paper is to investigate bifurcation of solutions of OZE and MS closure, applied to full LJ potential, in the gas-liquid co-existence region. This closure yields nearly self consistent results in most part of the phase plane except for low sub-critical temperatures. The method adopted is based on a discrete approximation of OZE in real space. The error reduction factor of algorithm scales as δ 4 , where δ is the mesh size. This approach is chosen in lieu of the more conventional approach in Fourier space because Fourier transform of OZE does not exist when χ 1 < 0 . Accurate asymptotic forms of correlation functions, which are essential for the real-space method, are incorporated so that the results are much less sensitive to cutoff length ( R 1 ) used in the method.
Main findings on bifurcation of solutions are as follows: (i) Spinodal point on the vapor side does not exist on any sub-critical isotherm, although there is a density on a nonphysical solution-branch where inverse compressibility vanishes. (ii) Fold bifurcation points exist on vapor and liquid sides of the isotherm, with a ‘no-real-solution’ region in between. (iii) Spinodal point exists on the liquid branch at a density larger than that of the fold bifurcation point on that side. (iv) In addition, solution-branches (violating the integral condition) with negative compressibility are found on vapor and liquid sides. (v) Complex solutions continuously connect the physical solution-branches (see Figure 3). (vi) It is shown that there are an infinite number of solutions when χ 1 is negative, which may also be seen as ‘singularity-induced’ bifurcation points. First four observations are quite similar to that found earlier, using real-space methods, for HNC closure [13].
The second objective in the paper is to relate results of CPE, in gas-liquid co-existence region, with those of restricted ensembles. This possibility is corroborated by displaying excellent agreement of pressure-density isotherms with simulation data, where the results show emergence of van der Waals loops in the sub-critical region. Analysis of a simple equation for singlet density distribution shows that it bifurcates to a spatially in-homogeneous distribution at the spinodal points. Perhaps, this is the root cause of all inconsistencies (negative compressibility, van der Waals loops, ‘no-real-solution’ region, etc.) in approaches to OZE in the co-existence region.

References

  1. J. P. Hansen, I. R. McDonald, Theory of Simple Liquids (Forth Edition) Academic Press, Cambridge. 2013.
  2. S. V. G. Menon, Bishnupriya Nayak, An equation of state for metals at high temperature and pressure in compressed and expanded volume regions Condensed Matter (MDPI) 4, 2019 , 71.
  3. C. Caccamo, Integral equation theory description of phase equilibria in classical fluids Physics Reports 274, 1996, 1.
  4. G. N. Sarkisov, Structure of simple fluid in the vicinity of the critical point: approximate integral equation theory of liquids J. Chem. Phys. 119, 2003, 373.
  5. N. Choudhury, S. K. Ghosh, Integral equation theory of Lennard-Jones fluids: A modified Verlet bridge function approach J. Chem. Phys. 116, 2002, 8517.
  6. G. Malescio, P. V. Giaquinta, Y. Rosenfeld, Iterative solutions of integral equations and structural stability of fluids Phys. Rev. E 57, 1998 , R-3723.
  7. C. A. Stuart, G. Vuillaume, A discrete nonlinear eigenvalue problem with many spurious branches of solutions Z. Anal. ihre Anwend 20, 2001, 183.
  8. G. Zerah, An efficient Newton’s method for numerical solution of fluid integral equations J. Comp. Phys. 61, 1985, 280.
  9. L. Belloni, A hypernetted chain study of highly asymmetrical poly-electrolytes Chem. Phys. 99, 1985, 43.
  10. R. O. Baxter, Method of solution of the Percus-Yevick, Hypernetted-chain or Similar eqautions Phys. Rev. 154, 1967 , 170.
  11. R. O. Watts, Percus-Yevick equation applied to a Lennard-Jones fluid J. Chem. Phys. 48, 1968, 50.
  12. L. Belloni, Inability of the hypernetted chain integral equation to exhibit a spinodal line J. Chem. Phys. 98, 1993, 8080.
  13. J. Kerins, L. E. Scriven, H. Ted Davis, Correlation functions in sub-critical fluid Adv. Chem. Phys. LXV, 1986 , 215.
  14. K. E. Atkinson, The numerical solution of integral equations on the second kind Cambridge University Press, Cambridge. 2009.
  15. K. E. Atkinson, L. F. Shampine, Algorithm 876: Solving Fredholm integral equations of the second kind in Matlab ACM Trans. Math. Softw. 34, 2008, 1.
  16. R. J. F. Leote de Carvalho, R. Evans, D. C. Hoyle, J. R. Henderson, The decay of the pair correlation function in simple fluids: long - versus short-ranged potentials J. Phys.: Condens. Matter 6, 1994, 9275.
  17. H. B. Keller, Lectures on numerical methods in bifurcation problems (Tata Institute of Fundamental Research) Springer, Berlin, Heidelberg, New York, Tokyo. 1986.
  18. R. Seydel, Practical Bifurcation and Stability Analysis (Third Edition) Springer, New York, London. 2010.
  19. M. E. Henderson, H. B. Keller, Complex bifurcation from real paths SIAM J. Appl. Math. 50, 1990, 460.
  20. A. Sai Venkata Ramana, S. V. G. Menon, Coupling-parameter expansion in thermodynamic perturbation theory Phys. Rev. E 87, 2013, 022101.
  21. S. V. G. Menon, Scaling of phase diagram and critical point parameters in liquid-vapor phase transition of metallic fluids Condensed Matter (MDPI) 6, 2021 , 6.
  22. S. V. G. Menon, Convergence of coupling-parameter expansion-based solutions to Ornstein-Zernike equation in liquid state theory Condensed Matter (MDPI) 6, 2021 , 29.
  23. V. K. Shen, J. R. Errington, Meta-stability and instability in the Lennard-Jones fluid investigated by transition-matrix Monte Carlo J. Phys. Chem. B 108, 2004, 19595.
  24. N. G. Van Kampen, Condensation of a classical gas with long-range attraction Phys. Rev. 135, 1964 , A362.
  25. J. D. Weeks, D. Chandler, H. C. Andersen, Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids J. Chem. Phys. 54, 1971, 5237.
  26. R. O. Baxter, Percus-Yevick equation of hard spheres with surface adhesion J. Chem. Phys. 49, 1968 , 2770.
  27. Menon, S. V. G.; Manohar, C.; Srinivasa Rao, K. A new interpretation of the sticky hard sphere model J. Chem. Phys. 1991, 95, 9188. [Google Scholar] [CrossRef]
  28. G. N. Sarkisov, E. Lomba, The gas-liquid phase transition singularities in the framework of the liquid-state integral equation formalism J. Chem. Phys. 122, 2005, 214504.
  29. A. T. Peplow, R. E. Beardmore, F. Bresme, Algorithms for the computation of solutions of Ornstein-Zernike Equation Phys. Rev. E 74, 2006, 046705.
  30. D. A. Tikhonov, G. N. Sarkisov, Singularities of solution of the Ornstein-Zernike Equation within the gas-liquid transition region Russian J. Phys. Chem. 74, 2000, 470.
  31. Waisman, E. The radial distribution function for a fluid of hard spheres at high densities Mean spherical integral equation approach. Mol. Phys. 1973, 25, 45. [Google Scholar] [CrossRef]
  32. Cummings, P. T.; Smith, E. R. Liquid-gas transition for hard spheres with attractive Yukawa tail interactions. Chem. Phys. 1979, 42, 341. [Google Scholar] [CrossRef]
  33. V. Venkatasubramanian, Singularity induced bifurcation and the van der pol Oscillator IEEE Trans. Circuits and Systems-I 41, 1994, 765.
  34. P. A. Monson, P. T. Cummings, Solution of the Percus-Yevick equation in the coexistence region of a simple fluid Int. J. Thermo. Phys. 6, 1985, 573.
  35. E. Lomba, J. L. López-Martin, On the solutions of the hypernetted chain equation inside the gas-liquid coexistence region J. stat. Phys. 80, 1995, 825.
  36. J. K. Johnson, J. A. Zollweg, K. E. Gubbins, The Lennard-Jones equation of state revisited Mol. Phys. 78, 1993, 591.
  37. J. P. Hansen, L. Verlet, Phase transitions of Lennard-Jones system Phys. Rev. 184, 1969, 151.
  38. L. Mier y Teran, A. H. Falls, L. E. Scriven, H. T. Davis, Structure and thermodynamics of stable, meta-stable, and unstable fluid - Efficient numerical method for solving the equations of equilibrium fluid physics Proc. Eighth Symp. on Thermo-physical Properties 8, 1982 , 45.
  39. R. Lovett, On the stability of a fluid toward solid formation J. Phys. Chem. 66, 1977, 66, 1225.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

Altmetrics

Downloads

224

Views

210

Comments

0

Subscription

Notify me about updates to this article or when a peer-reviewed version is published.

Email

Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2025 MDPI (Basel, Switzerland) unless otherwise stated