format: text

text:
# Implementation of the Hierarchical Reference Theory for simple one-component fluids

## I Introduction

The Hierarchical Reference Theory of fluids (hrt) pioneered by Parola and Reatto  has been found well capable  of describing structural and thermodynamic properties of fluids even in the vicinity of the critical point and for subcritical temperatures, yielding rigorously flat isotherms in the coexistence region (thus eliminating the need for Maxwell constructions) and non-classical values for the critical exponents . Still, adoption by a significant part of the liquid physics community of this renormalization group (rg) theoretical approach to the integral equation description of fluids has largely been lacking so far. While this may be partially attributed to hrt’s inherent difficulties and rather high computational cost, lack of an easy to use yet flexible, well-documented implementation of hrt may also have played a rôle. To fill this gap, we have written a program***Available on the world wide web from http://purl.oclc.org/NET/ar-hrt-1/. suited as a general framework for the exploration and application of hrt to simple one-component fluids with hard sphere reference systems with various combinations of physical systems, approximations, and solution algorithms. Within the natural limitations of the method, it has proved well applicable to a variety of model systems including the hard-core Yukawa (present contribution) and hard-core multi-step potentials while most attention has been devoted to the square well fluid .

Of course, ours is not the first implementation of hrt for simple one-component fluids: indeed, there has been a series of earlier programs  by the authors of the theory and their collaborators, but it was the one used in , henceforth referred to as the “original” implementation, that was a vital step in demonstrating the viability of hrt for continuous systems below the critical temperature; though never published or formally released, it has been circulating among interested physicists for quite some time. Our new software, on the other hand, differs from its precursors in many respects: adoption of a meta language in our version, programming style and documentation-to-code ratio may be most obvious, number and nature of hard-coded limitations, important details of the numerical procedure and a possible speed gain through generation of customized code might be less apparent. Most importantly, though, the original implementation’s structure makes experimentation with different combinations of approximations, partial differential equation (pde) solving algorithms, parameter settings and physical potentials rather cumbersome; in contrast, the fully modular approach adoption of a meta language allowed us to take seems far better suited to a more general survey of hrt’s numerous attractive features. In addition to the necessary flexibility of our software, great care has been taken to ensure the numerical soundness of every step in the calculation and hence of the results produced, so that the generation of numerical errors necessarily arising from finite-precision arithmetic is uniformly spread over all of the problem’s domain. To this end we introduce one central parameter, $`ϵ_\mathrm{\#}`$, characterizing the maximum relative error introduced at any step; together with a number of criteria relying on $`ϵ_\mathrm{\#}`$, this parameter governs virtually all of the numerics. Any deviation from this strategy is made explicit, as are all the other approximations entering the calculation. Ultimately, our goal was to provide the liquid physics community with a general and versatile yet numerically reliable tool for the systematic exploration and assessment of hrt and of the effects introduced by different approximations.

This paper is meant to serve a twofold purpose: to present the software we have written and its capabilities, and to provide its prospective users with some rudimentary documentation. To this end, after a brief presentation of the standard formulation of hrt for one-component fluids and some of the theory’s properties as far as they concern our implementation (section II), we first give a general outline of our program in section III, only touching upon the meta-language it has been written in. Due to our implementation’s fully modular design it is only natural to then proceed by a discussion of the more important of its building blocks and the various approximations they implement (sub-sections III B through III E). Presentation and critical assessment of the kind of results that can be attained with these and concluding remarks (section IV) are followed by two appendices dealing with some technical aspects of the formulation used.

## II The Theory

While a much more detailed account with additional references can be found in the review article , in what follows we want to limit ourselves to only a rough sketch of hrt; in doing so we are going to stress several aspects — some of them hardly discussed in the literature — relevant to our implementation of hrt, but only to the extent necessary for the discussion thereof. No prior knowledge of hrt is assumed.

The basic ingredient of hrt, already present in its precursor , is the gradual transition from a reference potential $`v^{\mathrm{ref}}(r)`$ to the full potential $`v(r)=v^{\mathrm{ref}}(r)+w(r)`$ describing the interaction between pairs of particles of a fluid, with any one of the intermediate potentials serving as a reference system with respect to which the properties of a successor potential are calculated \[superscripts always indicate the system a quantity refers to\]. In our work we have restricted ourselves to the case of a spherically symmetric pure two-body interaction, and we have taken advantage of the additional simplifications possible by identifying the reference system with a pure hard sphere system, $`v^{\mathrm{ref}}=v^{\mathrm{hs}}`$, as can always be achieved via the well-known Weeks-Chandler-Andersen scheme . Note however that this restriction to a hard sphere reference system is not present in the original implementation of hrt ; on the other hand, our program’s fully modular framework is flexible enough to accommodate any of these extensions should the need arise.

Other than , hrt achieves the transition from $`v^{\mathrm{ref}}`$ to $`v`$ in infinitesimally small steps. Inspired by momentum-space renormalization group (rg) theory, a cut-off wavenumber $`Q`$ varying from infinity to zero is introduced, and for every $`Q`$ the potential $`v^{(Q)}=v^{\mathrm{ref}}+w^{(Q)}`$ is defined such that Fourier components $`k<Q`$ of the perturbational part $`w^{(Q)}`$ of the $`Q`$-potential $`v^{(Q)}`$ are strongly suppressed whereas those for $`k>Q`$ coincide with those of the original potential $`w`$. Consequently, the reference system and the fully interacting system are recovered in the limits $`Q\mathrm{}`$ and $`Q0`$, respectively:

$$\begin{array}{ccc}\hfill v^{(\mathrm{})}& =& v^{\mathrm{ref}},\hfill \\ \hfill v^{(0)}& =& v.\hfill \end{array}$$
(1)
The rôle of the $`Q`$-potential just introduced becomes clear when we consider a functional expansion in $`\stackrel{~}{w}^{(Q)}`$ of thermodynamic and structural properties of the system with pair interaction $`v^{(Q)}`$ \[a tilde always denotes the Fourier transform\]: as $`\stackrel{~}{w}^{(Q)}(k<Q)`$ is small, the integrals in the expansion are effectively truncated for $`k<Q`$, in keeping with the rg picture.

In principle, the precise manner in which the potential is cut off should not matter, and one can easily conceive of many different ways of doing so. On the other hand, for such a procedure to be usable it must not introduce instabilities when truncating the hrt hierarchy, which is usually done at the two-particle level. Apart from approaches valid only for special types of potentials, we are aware of only two cut-off procedures suitable at least for attractive potentials (the standard formulation of hrt can easily be shown to become unstable for $`\stackrel{~}{w}(0)>0`$, as implicitly stated already in ); in our work we opted for the prescription presented in the review article  which seems to have been used almost exclusively so far  rather than the smooth cut-off formulation of , the latter being numerically cumbersome and predicting non-universal critical exponents. Thus we define the $`Q`$-potential $`v^{(Q)}=v^{\mathrm{ref}}+w^{(Q)}`$ by

$$\stackrel{~}{w}^{(Q)}(k)=\{\begin{array}{ccc}\stackrel{~}{w}(k)& :& k>Q\\ 0& :& k<Q;\end{array}$$
(2)
in $`r`$-space, $`w^{(Q)}`$ differs from $`w`$ by the addition of a convolution integral, viz.

$$w^{(Q)}(r)=w(r)\frac{1}{\pi r}_0^{\mathrm{}}\left(\frac{\mathrm{sin}Q(r^{}r)}{r^{}r}\frac{\mathrm{sin}Q(r^{}+r)}{r^{}+r}\right)r^{}w(r^{})dr^{}.$$

Obviously this $`Q`$-potential is a rather artificial function in $`r`$-space hardly resembling the full potential except in the limits of eq. (1); furthermore, the range in $`r`$-space over which $`v^{(Q)}`$ has to be considered is bound to be much larger than that of the original potential — a property immediately carrying over to related quantities, the direct correlation functions in particular. As an immediate consequence, numerical Fourier transformations involving the $`Q`$-potential or any of the correlation functions for the $`Q`$-system are computationally expensive and must be treated with extreme care; in fact, they should be avoided if possible at all, with obvious repercussions for the implementation of the core condition (v. i.).

In the transition from $`v^{\mathrm{ref}}`$ to $`v`$ or, equivalently, from $`Q=\mathrm{}`$ to $`Q=0`$, hrt treats the $`Q`$-system as reference for the properties of the system with infinitesimally lower cut-off $`Q\mathrm{d}Q`$; re-summation of terms in the resulting expansions in $`\mathrm{d}Q`$ and identification of quantities with a well-defined limit for $`Q0`$ finally yields the hrt equations : for every density $`\varrho `$ there is a formally exact hierarchy of coupled integro-differential equations involving a suitably modified free energy $`𝒜^{(Q)}`$ defined as

$$\frac{\beta 𝒜^{(Q)}}{V}=\frac{\beta A^{(Q)}}{V}\frac{\varrho ^2}{2}\left(\stackrel{~}{\varphi }(0)\stackrel{~}{\varphi }^{(Q)}(0)\right)+\frac{\varrho }{2}\left(\varphi (0)\varphi ^{(Q)}(0)\right)$$
(3)
($`\varphi =\beta w`$, $`\beta =1/k_BT`$; analogously, $`\varphi ^{(Q)}=\beta w^{(Q)}`$), a modified two-particle direct correlation function

$$𝒞^{(Q)}(r)=c_2^{(Q)}(r)+\varphi (r)\varphi ^{(Q)}(r),$$
(4)
and all higher order correlation functions $`c_n^{(Q)}(r)`$, $`n>2`$. The additional terms introduced in $`𝒜^{(Q)}`$ and $`𝒞^{(Q)}`$ explicitly take into account a discontinuity at $`Q=0`$ present in the unmodified free energy $`A^{(Q)}`$ and direct pair correlation function $`c_2^{(Q)}`$; as is apparent from eqs. (2), (3) and (4), modified and unmodified quantities coincide for $`Q=0`$. — Furthermore it should be noted that it is customary to include the ideal gas terms in the definition of the $`c_n^{(Q)}`$: for the two-particle case this is an additional term of $`\delta (\stackrel{}{r})/\varrho `$ so that the Ornstein-Zernike  equation takes the form

$$\stackrel{~}{c}_2=(1/\varrho )\varrho \stackrel{~}{h}\stackrel{~}{c}_2,$$
(5)
where $`h(r)`$ is the usual two-particle total correlation function; for the higher order correlation functions cf. .

The full hierarchy the derivation of which we just touched upon yields expressions for the total derivatives with respect to $`Q`$ for $`𝒞^{(Q)}`$ and all the $`c_n^{(Q)}`$, $`n>2`$; for the evolution of $`𝒜^{(Q)}`$ we have the particularly simple relation

$$\frac{\mathrm{d}}{\mathrm{d}Q}\left(\frac{\beta 𝒜^{(Q)}}{V}\right)=\frac{Q^2}{4\pi ^2}\mathrm{ln}\left(1\frac{\stackrel{~}{\varphi }(Q)}{\stackrel{~}{𝒞}^{(Q)}(Q)}\right).$$
(6)
As $`\mathrm{d}c_n^{(Q)}/\mathrm{d}Q`$, $`n2`$, involves $`𝒜^{(Q)}`$, $`𝒞^{(Q)}`$ and all higher order correlation functions up to $`c_{n+2}^{(Q)}`$ (so that, in particular, $`\mathrm{d}𝒞^{(Q)}/\mathrm{d}Q`$ depends on $`c_3^{(Q)}`$ and $`c_4^{(Q)}`$ via eq. (4)), the equations never decouple and we have to introduce some kind of closure. In doing so, it is usually desirable to retain thermodynamic consistency as embodied in the sum rule

$$\stackrel{~}{𝒞}^{(Q)}(0)=\frac{^2}{\varrho ^2}\left(\frac{\beta 𝒜^{(Q)}}{V}\right)$$
(7)
rigorously true for the exact solution of the full hierarchy; the derivatives with respect to $`\varrho `$ present in eq. (7) then mandate the transition from equations at fixed $`\varrho `$ to a pde in the $`(Q,\varrho )`$-plane with boundary conditions supplied at two densities, $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$. In addition, we need to retain the core condition, viz. $`g(r)=0`$ for $`r<\sigma `$ where $`g(r)=h(r)+1`$ is the pair distribution function; indeed it is one of hrt’s main advantages to conserve information on all length scales, ranging from the hard-sphere diameter $`\sigma (\varrho )`$ of the reference system at density $`\varrho `$ and the associated core condition up to the cut-off wave length $`1/Q`$, in the limit $`Q0`$ allowing criticality to arise from fluctuations of arbitrarily large wave length.

As noted above, the long-ranged nature of $`w^{(Q)}`$ and the correlation functions due to the cutting-off of eq. (2) is a strong argument in favor of any closure allowing an approximate implementation of the core-condition without the need for costly Fourier transforms. This is a likely reason for the up to now seemingly exclusive use of a closure in the spirit of the Lowest-Order $`\gamma `$-ordered Approximation (loga, ) or the equivalent Optimized Random-Phase Approximation (orpa, ) despite this closure’s known deficiencies : with the argument $`\varrho `$ silently to be added in earlier equations when used within the context of the pde, we make the ansatz

$$\begin{array}{cc}\hfill \stackrel{~}{𝒞}^{(Q)}(k,\varrho )& =\stackrel{~}{\varphi }(k,\varrho )+\gamma _0^{(Q)}(\varrho )\stackrel{~}{u}_0(k,\varrho )+\stackrel{~}{𝒦}^{(Q)}(k,\varrho ),\hfill \\ \hfill \stackrel{~}{𝒦}^{(Q)}(k,\varrho )& =\stackrel{~}{𝒢}^{(Q)}(k,\varrho )+\stackrel{~}{c}_2^{\mathrm{ref}}(k,\varrho ),\hfill \\ \hfill \stackrel{~}{𝒢}^{(Q)}(k,\varrho )& =_{n=1}^{\mathrm{}}\gamma _n^{(Q)}(\varrho )\stackrel{~}{u}_n(k,\varrho ),\hfill \end{array}$$
(8)
thereby introducing a set of $`Q`$-independent basis functions $`u_n`$ and corresponding expansion coefficients $`\gamma _n^{(Q)}`$. Here, $`u_0(r,\varrho )`$ is chosen proportional to $`w(r,\varrho )`$ (which has the undesirable effect that, from eqs. (8) and (4), the correlation function $`c_2^{(Q)}`$ of the $`Q`$-system depends upon the full potential $`w`$ rather than $`w^{(Q)}`$, as would be appropriate) and normalized so that $`\stackrel{~}{u}_0(0,\varrho )=1`$; the $`u_n(r,\varrho )`$, $`n1`$, on the other hand, are taken to form a basis for a suitable function space over $`[0,\sigma (\varrho )]`$. With these provisions, the problem of implementing both core condition and thermodynamic consistency reduces to that of an appropriate choice of the expansion coefficients $`\gamma _n^{(Q)}(\varrho )`$, $`n0`$, for every point in the $`(Q,\varrho )`$-plane. With the short-hand notations

$$\alpha ^{(Q)}(\varrho )=\frac{^3}{Q\varrho ^2}\left(\frac{\beta 𝒜^{(Q)}}{V}\right)$$

and, for an arbitrary function $`\psi (k,\varrho )`$,

$$\widehat{}^{(Q)}[\psi (k,\varrho ),\varrho ]=_^3\frac{\mathrm{d}^3k}{(2\pi )^3}\frac{\psi (k,\varrho )}{\left(\stackrel{~}{c}_2^{(Q)}(k,\varrho )\right)^2},$$
(9)
the condition (7) for thermodynamic consistency is easily re-written as

$$\frac{\gamma _0^{(Q)}(\varrho )}{Q}=\alpha ^{(Q)}(\varrho )\underset{n=1}{\overset{\mathrm{}}{}}\frac{\gamma _n^{(Q)}(\varrho )}{Q}\stackrel{~}{u}_n(0,\varrho ),$$
(10)
and following  the core condition can be shown to be equivalent to

$$\underset{n=0}{\overset{\mathrm{}}{}}\widehat{}^{(Q)}[\stackrel{~}{u}_j(k,\varrho )\stackrel{~}{u}_n(k,\varrho ),\varrho ]\frac{\gamma _n^{(Q)}(\varrho )}{Q}=\frac{Q^2}{2\pi ^2}\frac{\stackrel{~}{\varphi }(Q,\varrho )\stackrel{~}{u}_j(Q,\varrho )}{\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\left(\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\stackrel{~}{\varphi }(Q,\varrho )\right)},j1.$$

The latter can be combined with eq. (10) to the more explicit

$$\begin{array}{c}_{n=1}^{\mathrm{}}\widehat{}^{(Q)}[\stackrel{~}{u}_j(k,\varrho )\left(\stackrel{~}{u}_n(k,\varrho )\stackrel{~}{u}_0(k,\varrho )\stackrel{~}{u}_n(0,\varrho )\right),\varrho ]\frac{\gamma _n^{(Q)}(\varrho )}{Q}\\ =\alpha ^{(Q)}(\varrho )\widehat{}^{(Q)}[\stackrel{~}{u}_j(k,\varrho )\stackrel{~}{u}_0(k,\varrho ),\varrho ]+\frac{Q^2}{2\pi ^2}\frac{\stackrel{~}{\varphi }(Q,\varrho )\stackrel{~}{u}_j(Q,\varrho )}{\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\left(\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\stackrel{~}{\varphi }(Q,\varrho )\right)},j1;\end{array}$$
(11)
as both the sum rule and the core condition must hold for the reference system, the above evolution equations for the $`\gamma _n^{(Q)}`$ are readily supplemented with the initial conditions $`\gamma _n^{\mathrm{ref}}=0`$, $`n0`$.

As eq. (11) stands, it is no more amenable to direct implementation than the previous formulation; not only must this infinite-dimensional matrix equation be truncated to a finite number $`1+N_{\mathrm{cc}}`$ of basis functions, but even then the resulting integrals need to be evaluated at every $`Q`$ and $`\varrho `$ — a tedious process no less demanding than the Fourier transformations this approach is meant to replace. What makes this closure manageable, however, is the observation that the discontinuity in $`\widehat{}^{(Q)}`$’s integrand due to the appearance of $`\stackrel{~}{c}_2^{(Q)}`$ instead of the continuous $`\stackrel{~}{𝒞}^{(Q)}`$ leads to a term in $`\widehat{}^{(Q)}[\psi (k,\varrho ),\varrho ]/Q`$ made up of functions evaluated at $`k=Q`$ alone; following , only this single term is retained, leading to the approximation

$$\frac{}{Q}\widehat{}^{(Q)}[\psi (k,\varrho ),\varrho ]=\psi (Q,\varrho )\frac{Q^2}{2\pi ^2}\frac{2\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\stackrel{~}{\varphi }(Q,\varrho )\left(\stackrel{~}{\varphi }(Q,\varrho )\right)^2}{\left(\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\right)^2\left(\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\stackrel{~}{\varphi }(Q,\varrho )\right)^2}$$
(12)
leaving out the non-local contribution $`2_{n=0}^{\mathrm{}}\widehat{}^{(Q)}[\psi (k,\varrho )\stackrel{~}{u}_n(k,\varrho )/\stackrel{~}{c}_2^{(Q)}(k,\varrho ),\varrho ](\gamma _n^{(Q)}(\varrho )/Q)`$; with the above approximation, the task of evaluating one of the integrals of eq. (11) reduces to only an initial integration for the reference system followed by the solution of an ordinary differential equation (ode) coupled to the hrt-pde as well as analogous odes for all the other integrals of the $`\widehat{}`$-type.

Of course, to fully specify the mathematical problem, the pde must be amended by both initial and boundary conditions; while the former take the simple form of vanishing expansion coefficients (v. s.), the latter also impose some constraint on the $`\gamma _n^{(Q)}`$. However, as long as we retain the core condition, such an additional constraint is already sufficient to determine the expansion coefficient $`\gamma _0^{(Q)}`$; unless the $`\gamma _0^{(Q)}`$ so found exactly reproduces eq. (10), thermodynamic consistency can no longer be imposed without introducing mathematical inconsistencies. By the same token, eq. (11), derived by incorporating the sum rule (7) into the core condition, is no longer valid but must be changed to

$$\begin{array}{cc}\hfill _{n=1}^{\mathrm{}}\widehat{}^{(Q)}[\stackrel{~}{u}_j(k,\varrho )\stackrel{~}{u}_n(k,\varrho ),\varrho ]\frac{\gamma _n^{(Q)}(\varrho )}{Q}=& \widehat{}^{(Q)}[\stackrel{~}{u}_j(k,\varrho )\stackrel{~}{u}_0(k,\varrho ),\varrho ]\frac{\gamma _0^{(Q)}(\varrho )}{Q}\hfill \\ & +\frac{Q^2}{2\pi ^2}\frac{\stackrel{~}{\varphi }(Q,\varrho )\stackrel{~}{u}_j(Q,\varrho )}{\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\left(\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )\stackrel{~}{\varphi }(Q,\varrho )\right)},j1.\hfill \end{array}$$
(13)
to reflect the transition from eq. (7) to said constraint determining the $`\gamma _0^{(Q)}`$ appearing on the above equation’s right hand side; furthermore, elimination of thermodynamic consistency obviously means decoupling the pde to a set of odes at fixed density. This is exactly what is needed at the boundaries $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$ of the density interval considered where the $`\varrho `$-derivatives defining $`\alpha ^{(Q)}(\varrho )`$ cannot be evaluated: While specialization to $`\varrho _{\mathrm{min}}=0`$ uniquely determines the solution there by virtue of the divergence of the ideal gas term $`1/\varrho `$ in $`\stackrel{~}{c}_2^{\mathrm{ref}}`$ (cf. appendix A), some prescription for finding $`\gamma _0^{(Q)}`$, accompanied by the necessary modification of the truncated matrix equation according to eq. (13), must be imposed at $`\varrho _{\mathrm{max}}`$; as long as we retain eq. (11) in the pde’s domain’s interior and numerically evaluate $`\alpha ^{(Q)}(\varrho )`$ there, it is natural to use the loga/orpa condition $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})=0`$  at the boundary, which is sufficient to also determine $`𝒜^{(Q)}(\varrho _{\mathrm{max}})/Q`$ (or, equivalently, appendix A’s $`f(Q,\varrho _{\mathrm{max}})`$) from the $`\gamma _n^{(Q)}(\varrho _{\mathrm{max}})`$, $`n1`$, alone.

Unfortunately, it turns out that a scheme retaining $`\alpha ^{(Q)}(\varrho )`$ in eq. (11) for $`\varrho _{\mathrm{min}}<\varrho <\varrho _{\mathrm{max}}`$ presents significant numerical problems for all but extremely high temperatures, precluding reaching $`Q=Q_0`$ at least for the potentials that we have looked at. This is where the so-called “decoupling assumption” comes into play: based upon the different ranges of $`u_0(r)`$ and $`u_n(r)`$, $`n1`$, the authors of  argue that $`\alpha ^{(Q)}(\varrho )=0`$ might be a good approximation, thus eliminating the $`\widehat{}`$-integral on the right hand side of eq. (11); it turns out that this change, invariably adopted in all later publications, is often sufficient to allow generating a solution all the way to $`Q=Q_0`$. From our previous discussion of conditions imposed on the $`\gamma _n^{(Q)}`$ it should be obvious that this decoupling assumption is incompatible not only with the loga/orpa condition $`\gamma _0^{(Q)}(\varrho )=0`$ retained in the original implementation for $`\varrho =\varrho _{\mathrm{max}}`$ but also with thermodynamic consistency (eq. (7)) altogether; thus we are left with only a few possibilities: we may either retain logical consistency by using the decoupling assumption $`\alpha ^{(Q)}(\varrho )=0`$ as a closure for the hrt-equations, reducing the pde to a set of odes in $`Q`$ only; or we may prefer to retain the structure of a pde so as to make use of thermodynamic consistency at least to a certain degree; yet another possibility is to maintain both mathematical and thermodynamic consistency by not implementing the core condition at all. The original implementation’s approach relying on three mutually incompatible concepts, viz. the loga/orpa condition at $`\varrho _{\mathrm{max}}`$, decoupling, and the compressibility sum rule, seems particularly unattractive; at least one should use the decoupling assumption as a boundary condition at high density instead.

Retaining thermodynamic consistency in the form of eq. (10) as well as, in an approximate, way, the core condition via the truncated eq. (11) or (13) together with the approximation (12), we thus arrive at a set of equations implementing hrt with the loga/orpa-like closure (8) on the two-particle level well suited for numerical processing. While these expressions lend themselves to discretization in a straightforward way, it is computationally much more convenient to cast the pde in a form superficially resembling a quasi-linear one  so that an implicit finite-difference scheme requires only the inversion of a tridiagonal matrix. The re-writing we adopted — detailed in appendix A, very similar to the one of  — results in the introduction of an auxiliary function $`f(Q,\varrho )`$ so that the pde implied by eqs. (6) and (7) can be written in the form

$$\frac{}{Q}f(Q,\varrho )=d_{00}[f,Q,\varrho ]+d_{01}[f,Q,\varrho ]\frac{}{\varrho }f(Q,\varrho )+d_{02}[f,Q,\varrho ]\frac{^2}{\varrho ^2}f(Q,\varrho ).$$
(14)
Now it is easy to demonstrate the pde’s stiffness for subcritical temperatures: recalling the definitions of appendix A, we find that the inverse isothermal compressibility $`1/\kappa _T`$ of the system with potential $`v`$ can be written as $`1/\kappa _T=\varrho ^2\stackrel{~}{w}(0,\varrho )/\overline{\epsilon }(0,\varrho )`$, where $`\overline{\epsilon }(Q,\varrho )+1=\epsilon (Q,\varrho )\mathrm{exp}(f(Q,\varrho )\stackrel{~}{u}_0^2(Q,\varrho ))`$. For subcritical temperatures there is a density interval $`[\varrho _v,\varrho _l]`$ where $`1/\kappa _T=0`$ (implying diverging $`\overline{\epsilon }`$ and hence $`f`$); here, $`\varrho _v`$ and $`\varrho _l`$ are the densities of the coexisting vapor and liquid, respectively (recall that hrt yields rigorously flat isotherms within the coexistence region, with binodal and spinodal coinciding in three dimensions ). By construction, however, the limit $`Q0`$ is a continuous one (cf. eqs. (3) and (4)) so that $`f`$, $`\epsilon `$ and related quantities must be large already well before $`Q=0`$ is reached; at the same time, the rg mechanism introduced by the definition of $`w^{(Q)}`$ via eq. (2) precludes any divergence at non-vanishing $`Q`$. Considering the region of the $`(Q,\varrho )`$-plane where $`f`$ and $`\epsilon `$ are large, we easily find from the explicit expressions of appendix A that the $`d_{0i}`$, and hence $`f(Q,\varrho )/Q`$, are of order $`\epsilon ^1`$; restricting ourselves to a specific density $`\varrho `$ and sufficiently small $`Q`$ (so that $`\stackrel{~}{u}_0(Q,\varrho )=1`$; extension of the argument to a larger $`Q`$-range is cumbersome but straightforward) we can write

$$\frac{\mathrm{d}f(Q,\varrho )}{\mathrm{d}Q}=e^{f(Q,\varrho )}d_0(Q),$$

where $`d_0`$ is now of order unity. Inspection of the solution of this ode immediately shows that the average of $`d_0(Q)`$ over the interval $`[Q_1,Q_2]`$, $`0<Q_1<Q_2`$, is rigorously bounded from above by $`\mathrm{exp}(f(Q_2,\varrho ))/(Q_2Q_1)`$ or else there were a singularity in that $`Q`$-interval; translating back to $`f`$ we see that, while $`|f(Q,\varrho )/Q|`$ is still of order $`\epsilon `$, $`f`$ must be a rapidly oscillating function of $`Q`$ (with a period of order $`\epsilon ^1`$), the average slope of which is much smaller, viz. of order $`1/Q`$.

It should be noted that this stiffness is not an artifact of the re-formulation of the pde as summarized in appendix A but is manifest just the same when directly solving the pde for the modified free energy $`𝒜^{(Q)}(\varrho )`$ rather than that for the auxiliary $`f(Q,\varrho )`$ . The above argument relies only on some general properties of hrt in the current formulation as applied to one-component fluids: the divergence of the isothermal compressibility in the coexistence region (the reproduction of which is one of hrt’s main achievements), continuity of the limit $`Q0`$, and the suppression of divergences for $`Q>0`$ as a result of the rg mechanism implemented via the truncated potential $`v^{(Q)}`$. The essential additional ingredient, viz. the behavior of the ratio of the $`Q`$\- and the $`\varrho `$-derivatives as the divergence in the compressibility builds up, while obvious in the formulation via $`f`$, is not easily seen in terms of $`𝒜^{(Q)}(\varrho )`$; this, however, comes as no surprise since $`f`$ is essentially the free energy’s derivative with respect to $`Q`$ so that we would have to reason about third- and second-order derivatives rather than first- and second-order ones if we were to repeat the arguments without resorting to the re-writing of appendix A.

## III Overview of the program

With the theory and approximations outlined in the previous section, we are now in a position to undertake the task of implementing hrt for one-component fluids in fully standards conforming Fortran-90; the only non-standard feature we make use of is the availability of the special values NaN and $`\pm \mathrm{𝙸𝚗𝚏}`$ for numerically undefined values and signed overflows, respectively, as defined in the floating-point standard IEEE 754. These requirements should not pose a serious restriction for our program’s possible users: after all, Fortran-90 compilers have been available for a wide range of platforms for several years, and the desired floating-point behavior can usually be requested — albeit at a small performance penalty — via compiler switches. While our implementation is more appropriately described as a collection of mutually compatible building blocks rather than as a monolithic program so that the details of the numerical procedure are best left to these parts, for the combination of different selections to work all versions of all the modular constituents must adhere to a common model of the computation:

Most obviously, we have to make the transition from the pde’s domain, viz. the infinite strip $`[0,\mathrm{})\times [\varrho _{\mathrm{min}},\varrho _{\mathrm{max}}]`$, to a discrete mesh defined by a finite number of points in the $`(Q,\varrho )`$-plane. Evidently, the placement of these “nodes”, as we shall call them, is of utmost importance for the quality of the discretization so that it is only natural to define $`ϵ_\mathrm{\#}`$, the central parameter governing all of the numerics, in terms of the properties of this mesh: the coarser a mesh we chose, the larger $`ϵ_\mathrm{\#}`$ will be. — In principle, the locations of the nodes, the data structures of which are organized in linked lists, can be chosen freely; in particular, the cut-offs $`Q`$ of all the systems in such a list of nodes are not taken to necessarily coincide, even though this is usually the case except for a low-density boundary at $`\varrho =0`$. As for the densities of the nodes, implementation of the core condition via the truncated eq. (11) and eq. (12) makes anything but constant (though not necessarily equispaced) density values impractical; if the grid is to be refined for low $`Q`$, additional nodes must be inserted at the same densities in all the node lists in the calculation. — After initialization of the nodes’ data structures, solution of the pde proceeds by applying a (possibly iterated) predictor-corrector scheme to generate an approximate solution for the nodes most advanced towards $`Q=0`$ from the information available through the node lists at higher $`Q`$; in the interest of the code’s simplicity, the number of such node lists has been fixed to exactly three, thus facilitating determination of appropriate step sizes $`\mathrm{\Delta }Q`$.

As mentioned before, the goals set for our implementation of hrt necessitate a fully modular architecture of our program; and while we did not want to forgo the well-known advantages of Fortran-90 for the numerical work, experience with prior versions of our code taught us that the kind of flexibility we need cannot be accommodated within the rather rigid framework Fortran-90’s modules with their one-way flow of information provide. Instead we opted for a simple meta-languageAvailable on the world wide web from http://purl.oclc.org/NET/arfg/. for self-configuring construction of code customized to the chosen combination of approximations and the physical system at hand, at the same time enhancing readability and maintainability of the source and encouraging modularization; for a more detailed discussion of this approach and the numerous technical advantages it affords we refer the reader to .

### A Main parts

Our software can best be understood as a collection of mutually compatible and freely exchangeable building blocks corresponding to the underlying physical and mathematical notions; the resulting natural organization of the code cleanly separating conceptually unrelated approximations is a direct consequence of our adoption of a meta-language and the use of automatic code generation techniques. The implementation’s modular constituents, henceforth dubbed “main parts”, must, however, be clearly distinguished from Fortran-90’s modules: in general, there is no simple mapping from main parts to modules, and every main part may give rise to any number of modules, incorporating all the information available within the code base.

In the following sub-sections we take a closer look at some of the main parts, their physical meaning, the algorithms and approximations implemented, and at some of the information they make available to the other parts; we will, however, exclude from this discussion the program’s infrastructure, e. g. the implementation of logging, of reading and parsing of options files, handling of node lists and the definition of a versatile, lossless and storage-efficient albeit platform dependent file format for the results at $`Q=Q_0`$. In a similar vein, we only mention the assortment of accompanying tools for reading these files and dumping their content in human-readable or Mathematica-usable form, for locating the critical point or calculating phase diagrams. — Thus only main parts “potential”, “reference”, “ansatz” and “solver” remain to be discussed:

### B Properties of the potential

First and foremost, we obviously have to provide information on the fluid’s potential $`v=v^{\mathrm{ref}}+w`$ and its properties: this is the purpose of the main part labeled “potential”. Just as the full potential is a sum of a reference part $`v^{\mathrm{ref}}`$ and a perturbational part $`w`$, the functions and parameters to be provided by this main part fall into two distinct categories, pertaining to either $`v^{\mathrm{ref}}`$ or $`w`$; in addition, as the temperature enters the calculation only as a pre-factor to $`w`$, viz. via $`\varphi =\beta w`$, the inverse temperature $`\beta `$ is also defined here.

As far as the reference system is concerned, restriction to hard spheres (for the rationale cf. section II) means that only a function returning the hard sphere diameter $`\sigma (\varrho )`$ and a flag indicating any deviation of $`\sigma (\varrho )`$ from the unit of length need to be made available. — A similar parameter pertaining to the perturbational part $`w`$ of the potential, viz. a flag indicating any density-dependence of $`w`$, also plays an important rôle in many parts of the program as substantial simplifications and, in many cases, significant speed-ups by caching previous results are possible whenever $`\stackrel{~}{\varphi }(Q,\varrho )`$ only depends on $`Q`$. In addition, at every cut-off $`Q`$ the program must have access to the Fourier transforms $`\stackrel{~}{w}(Q,\varrho )`$ and $`\stackrel{~}{\varphi }(Q,\varrho )`$ as well as the derivatives $`\stackrel{~}{\varphi }(Q,\varrho )/Q`$ and $`^n\left(\stackrel{~}{\varphi }(Q,\varrho )\right)^m/\varrho ^n`$, whereas powers of the volume integral, $`\stackrel{~}{\varphi }(0,\varrho )^n`$, and their derivatives $`^n\left(\stackrel{~}{\varphi }(0,\varrho )\right)^m/\varrho ^n`$ obviously do not depend on $`Q`$ (here, $`m`$ and $`n`$ are appropriate integers known during code construction).

For the benefit of the pde-solving algorithm, this main part also has to set a parameter $`\lambda _{[v]}`$ related to the maximum relative curvature of the second $`Q`$-derivative of $`\stackrel{~}{\varphi }(Q,\varrho )`$, defined in such a way as to coincide with $`\lambda `$ for the square well potential $`v^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}`$. Indeed, most of our efforts so far  have concentrated on this particularly simple type of potential given by

$$w^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}(r)=\{\begin{array}{ccc}ϵ& :& r<\lambda \sigma \\ 0& :& r>\lambda \sigma ,\end{array}$$

where the hard sphere diameter $`\sigma `$ of the reference part and the strength $`ϵ`$ of the potential are usually chosen as units of length and energy, respectively. — Another type of potential we have implemented is a generalized step potential, i. e. a succession of stretches of constant $`w(r)`$; more specifically, the perturbational part $`w(r)`$ of an $`n`$-step potential of this type is just a sum of square-well potentials, $`w=_{i=1}^nw^{\mathrm{sw}[ϵ_i,\lambda _i,\sigma ]}`$, $`\lambda _i<\lambda _{i+1}`$, $`1i<n`$; again, we have only considered the $`\varrho `$-independent case, while all the potential’s parameters should be assumed functions of $`\varrho `$ when modeling a specific physical system. — We have also implemented a $`\varrho `$-independent hard-core Yukawa potential $`v^{\mathrm{hcy}}`$,

$$w^{\mathrm{hcy}[ϵ_0,ϵ,z,\sigma ]}(r)=\{\begin{array}{ccc}ϵ_0& :& r<\sigma \\ ϵ\frac{\sigma }{r}e^{z(r\sigma )}& :& r>\sigma ,\end{array}$$

where the parameter $`ϵ_0`$, the value of $`w`$ inside the core, defaults to $`ϵ=w^{\mathrm{hcy}}(\sigma +)`$, which again is usually chosen to coincide with the unit of energy; any mismatch between $`ϵ_0`$ and $`ϵ`$ dominates $`\stackrel{~}{w}^{\mathrm{hcy}}(k)`$ for large $`k`$ and is found to render unstable at least the numerics.

Regarding the stability of the pde, recall from section II that an attractive potential (so that $`\stackrel{~}{w}(0,\varrho )<0`$) is a necessary though not sufficient condition for the stability of the pde.

### C Hard sphere reference system

Due to the specialization of $`v^{\mathrm{ref}}`$ to hard spheres, the reference system enters the expressions of section II only through the direct correlation function $`c_2^{\mathrm{ref}}`$, implementation of which is the task set for main part “reference”: drawing upon the information from main part “potential”, only some initialization code and functions for the evaluation of $`\stackrel{~}{c}_2^{\mathrm{ref}}(Q,\varrho )`$ and $`\stackrel{~}{c}_2^{\mathrm{ref}}(Q,\varrho )/Q`$ have to be exported.

In our program we have so far included two different versions implementing the Percus-Yevick  approximation and the Henderson-Grundke  description; note that all results reported here have been obtained using the Henderson-Grundke $`c_2^{\mathrm{ref}}`$: in a theory relying on internal consistency conditions like eq. (7) as heavily as hrt, the thermodynamic inconsistency present in the Percus-Yevick solution seems particularly undesirable.

### D Discretization, boundary conditions, and other approximations

Main part “ansatz”, where all the approximations on the physical and mathematical level are combined to jointly define a reasonable numerical model of hrt, is at the very core of the pde-solving machinery: for the potential the perturbational and reference parts of which have been described in the two previous sub-sections, the hrt-pde is discretized and solved according to a given set of approximations and on the mesh defined by the node-lists served by main part “solver” (sub-section III E). More precisely, “ansatz” provides a set of facilities in the form of subroutines with standardized interfaces implementing the various stages of the computation, viz. initialization of the node lists at $`Q=Q_{\mathrm{}}`$ and solution of the pde according to a predictor-corrector full approximation scheme. Note, however, that the code must accommodate the possibilities of both iterating the corrector step (which may allow reaching the numerical quality indicated by $`ϵ_\mathrm{\#}`$ with somewhat larger step sizes, thus speeding up the calculation) and of discarding part of the solution should $`ϵ_\mathrm{\#}`$-based criteria not be met; to aid “solver” in these decisions, care has to be taken to detect and signal numerical anomalies. Once a step’s results have been accepted, “ansatz” may perform additional manipulations of the data structures; most importantly, the re-scaling of all quantities affected by exponentiation of $`f`$ necessary whenever $`f`$ is large (cf. our discussion of the pde’s stiffness in section II) is adjusted only when the last corrector’s result has been accepted.

Due to the eminent rôle of the consistency condition eq. (7) in constructing a closure to the underlying ode (6), the pde (14) for $`f(Q,\varrho )`$ is of first order in $`Q`$ and of second order in $`\varrho `$; assuming the lowest possible number of nodes in the discretization (extension to higher order is straightforward), we need at least a $`2\times 3`$ set of nodes. According to the general model of the computation presented in section III, however, we instead keep a third node list in order to allow monitoring of second $`Q`$-derivatives, so that we use a discretization on the $`3\times 3`$ grid schematically presented in fig. 1, including information available via that additional node list. Locally, the discretization is derived from an expansion about the midpoint of the nodes labeled $`(22)`$ and $`(32)`$ in the schematic 1, evaluating the second $`\varrho `$-derivative along the line of constant $`Q`$ through this point (thin horizontal line in fig. 1) by estimating the data at the intersection with the lines of constant density by interpolants defined from node triples $`(i1)`$ and $`(i3)`$, respectively; the resulting finite difference approximation is applied to every set of three adjacent node-triples, substituting suitable boundary conditions at $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$.

As indicated in fig. 1, $`Q`$ is not necessarily constant along a given node list, whereas the stability of the numerical scheme may impose certain geometrical constraints regarding the possible locations of the nodes, e. g. for ensuring that the Courant-Friedrichs-Lewy criterion  is met or for maintaining convexity of the remaining integration region; a suitable representation of these constraints is exported and must be taken into account by main part “solver”. If the latter decides to insert nodes at intermediate densities, the code for initializing the inserted data structures and for interpolating appropriate quantities is negotiated between the main parts, depending upon the order of the interpolation formulæ available. A further consequence of having non-constant $`Q`$ is that some parts of the density range may reach $`QQ_0`$ earlier than others; in this case, the corresponding nodes are locked, preventing further modification, and all of the converged nodes except those necessary for providing a boundary condition for the remaining density interval are removed from the node lists available to main part “ansatz”.

In addition to the discretization of the hrt-pde (14) discussed so far, the implementation of the core condition along the lines of section II and appendix A is also of interest. Relegating discussion of the choice of appropriate basis function $`u_n`$, $`1nN_{\mathrm{cc}}`$, to appendix B, we only point out the extremely slow convergence of the $`\widehat{}^{(Q)}`$-integrals (9) that have to be evaluated at $`Q=Q_{\mathrm{}}`$; furthermore, as the integrand is temperature dependent for $`k>Q_{\mathrm{}}`$, these integrals have to be evaluated for every isotherm — a problem that might be sidestepped by adopting the original implementation’s strategy of consistently using the results for $`Q\mathrm{}`$ rather than those valid at $`Q_{\mathrm{}}`$ for initialization even though such an approach introduces a discontinuity at $`Q=Q_{\mathrm{}}`$. Also, with the usual choice of $`Q_{\mathrm{}}10^2/\sigma `$, integration merely up to $`k=Q_{\mathrm{}}`$ can hardly be deemed sufficient; an appropriate upper integration limit can instead be found by comparing the integrand’s asymptotic behavior with $`ϵ_\mathrm{\#}`$. — Main part “ansatz” also has to identify quantities suitable both for monitoring convergence of the full approximation scheme and for choosing appropriate step sizes $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$, and to make available code fragments for the inspection of nodes in various stages of the computation as well as a description of the boundary conditions at $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$ including mandatory settings for either of these parameters if necessary; in particular, most implementations require $`\varrho _{\mathrm{min}}=0`$ in order to be able to use the divergence of the ideal-gas term in $`\stackrel{~}{c}_2^{\mathrm{ref}}`$ as a $`Q`$-independent boundary condition for $`f`$ (cf. appendix A).

It is this main part that defines the formulation of hrt and the set of approximations used in the calculations; of the numerous versions of this main part we produced only a few are to be mentioned here: both the re-implementation of the original program’s approximations for the core and boundary conditions and the approach combining the pde with $`\alpha ^{(Q)}=0`$ at all densities including $`\varrho _{\mathrm{max}}`$, while mathematically inconsistent, retain thermodynamic consistency at least in some approximate way (cf. our discussion of the decoupling assumption in section II); in addition to these, the two possible approaches at least mathematically meaningful, viz. the odes directly following from decoupling and the pde resigning on the core condition for the benefit of the compressibility sum rule (7) (with the loga/orpa prescription $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})=0`$ as high density boundary condition) will also be used in section IV’s presentation of the results our software yields.

### E Criteria for positioning of nodes

If main part “ansatz” is to provide a discretization on whatever mesh is handed to it, it is the task set for “solver”, the last of the main parts to be discussed here, to define this very mesh and to keep track of the numerical solution’s quality. Based primarily upon the value of $`ϵ_\mathrm{\#}`$ and respecting any restrictions exported by “ansatz”, step sizes $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$ have to be chosen and checked for compatibility with the solution generated, iterating or discarding steps if certain criteria are not met; whenever “ansatz” signals an exception the last step is discarded, accepting the data in the node list corresponding to labels $`(2i)`$ in fig. 1 as the best approximation to the solution for $`Q0`$. At the same time, care has to be taken to locate and identify any problems in the solution, i. e. parts of the $`(Q,\varrho )`$-plane where the solution found does not appear smooth on the scales set by the step sizes, the most basic assumption underlying any finite-difference calculation; whenever this assumption no longer holds, the algorithm will react by locally reducing $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$, inserting node triples (cf. sub-section III D) in order to achieve the latter. — Once we find any nodes already holding the final results for their respective densities they must be taken care of as discussed in sub-section III D; integration of the pde is ended when there is a node with $`QQ_0`$ for every density in the calculation, or when “ansatz” requests an end either because an error condition has occurred (v. s.) or because the current node list is sufficiently close to $`Q=Q_0`$ already. — As noted in section III, the intimate link between this main part’s task and the numerical quality of the solution generated makes it natural to here define $`ϵ_\mathrm{\#}`$, the central parameter governing the numerics, and it is this part of the program that relies upon $`ϵ_\mathrm{\#}`$ and the associated criteria the most; other main parts use $`ϵ_\mathrm{\#}`$ for little more than for switching between full analytic expressions and asymptotic expansions (for a slightly atypical example of which cf. appendix B).

Of the two implementations of this main part, one has been written in the hopes of being able to avoid the problematic region of large $`f(Q,\varrho )`$ altogether, as is indeed possible for some similar pdes. This implementation makes full use of $`ϵ_\mathrm{\#}`$, relying on numerous criteria to control the calculation; in the following discussion the notation $`p_y^{[x]}`$ refers to customization parameters that should usually be taken as real numbers of order unity. A case in point is the choice of the density grid: even though this is not necessary, we decided to always start with an equispaced set of $`\varrho `$-values ranging from $`\varrho _{\mathrm{min}}`$ to $`\varrho _{\mathrm{max}}`$; the number $`N_\varrho +1`$ of such density values is related to $`ϵ_\mathrm{\#}`$ by

$$N_\varrho =\frac{(\varrho _{\mathrm{max}}\varrho _{\mathrm{min}})^2}{ϵ_\mathrm{\#}p_{N_\varrho }^{[\varrho ]}},$$

reflecting the importance of second $`\varrho `$-derivatives in the solution of the pde as well as the static nature of this set of densities due to the $`\widehat{}`$-approach to the core-condition. — Once $`ϵ_\mathrm{\#}`$ has been determined, the system is ready to start determining appropriate step sizes $`\mathrm{\Delta }Q`$; in particular the assumption that the potential $`v(r,\varrho )`$ introduces length scales only in the range from $`\sigma (\varrho )`$ to $`\lambda _{[v]}(\varrho )\sigma (\varrho )`$, where $`\lambda _{[v]}`$ is related to the fourth derivative of $`\stackrel{~}{v}(k,\varrho )`$ with respect to $`k`$ (cf. sub-section III B), places an upper bound on the admissible step sizes, viz.

$$\mathrm{\Delta }Q\frac{\sqrt{12ϵ_\mathrm{\#}p_{\mathrm{\Delta }Q_{\mathrm{max}}}^{[\mathrm{\Delta }Q]}}}{\lambda _{[v]}\sigma (\varrho )}.$$

On the other hand, for a finite difference scheme to be meaningful at least a certain number of bits must remain significant in evaluating the differences, which implies a lower bound on $`\mathrm{\Delta }Q`$ proportional to $`Q`$, and the solution has to be smooth on the scales defined by the mesh, which also rules out abrupt changes in the step sizes; consequently, the ratio of two consecutive $`\mathrm{\Delta }Q`$ steps at the same density is restricted to lie between $`p_{\mathrm{ratio}}^{[\mathrm{\Delta }Q]}`$ and $`1/p_{\mathrm{ratio}}^{[\mathrm{\Delta }Q]}`$. In a similar vein, considering smoothness in the $`\varrho `$-direction we have to postulate that $`(Q_{(22)}+Q_{(32)})/2`$ is greater than either of $`Q_{(31)}`$ and $`Q_{(33)}`$, where the labels coincide with those of the nodes of fig. 1; this condition, unlike the other rules mentioned so far, does not limit the step sizes $`\mathrm{\Delta }Q`$ at any density $`\varrho `$ but rather determines whether $`\mathrm{\Delta }\varrho `$ should be reduced by the insertion of nodes at an additional density. But the most important criteria for choosing $`\mathrm{\Delta }Q`$ come from monitoring the solution generated: for every monitored quantity $`x`$ we make sure that

$$\begin{array}{c}\sqrt{\frac{1}{x_Q}\left|\frac{^2x}{Q^2}\right|}\mathrm{\Delta }Q\sqrt{ϵ_\mathrm{\#}p_x^{[\mathrm{\Delta }Q]}},\\ x_Q=\mathrm{max}_{k>Q}|x(k,\varrho )|;\end{array}$$

the quantities taken for $`x`$ are, of course, defined by “ansatz”, and a usual choice is $`x(Q,\varrho )\{f(Q,\varrho ),1/\stackrel{~}{𝒦}^{(Q)}(Q,\varrho )\}`$. A different set of quantities $`y`$, also chosen by “ansatz” and usually comprising just $`y(Q,\varrho )=f(Q,\varrho )`$, is used to monitor the convergence of the predictor-corrector scheme and to determine whether or not the corrector should be iterated: denoting the absolute difference of consecutive approximations of $`y`$ divided by $`y_Q`$ by $`\mathrm{\Delta }y`$, iterations are performed until $`\mathrm{\Delta }y<ϵ_\mathrm{\#}p_y^{[\mathrm{conv}]}`$, and the ratio of two consecutive step sizes is bounded from above by $`(ϵ_\mathrm{\#}p_y^{[\mathrm{\Delta }Q]})^{1/\left(2+p_{N_{\mathrm{it}}}^{[\mathrm{\Delta }Q]}\right)}/\sqrt{\mathrm{\Delta }^{(1)}y}`$, where $`\mathrm{\Delta }^{(1)}y`$ is $`\mathrm{\Delta }y`$ evaluated after the first corrector step. According to simple heuristic arguments regarding the convergence of corrector iterations and ignoring the effect of other criteria, an average of $`p_{N_{\mathrm{it}}}^{[\mathrm{\Delta }Q]}`$ calls of the corrector can be expected to solve the difference equations to within $`ϵ_\mathrm{\#}`$, and a setting of $`p_{N_{\mathrm{it}}}^{[\mathrm{\Delta }Q]}>1`$ may significantly speed up the calculation by allowing larger steps to be taken without loss of accuracy. — After finding and tentatively using a candidate $`\mathrm{\Delta }Q`$, we still have to check that the assumptions leading to that particular choice for $`\mathrm{\Delta }Q`$ actually hold; to this end we re-evaluate all the criteria with the obvious exception of the one involving $`\mathrm{\Delta }^{(1)}y`$ after the predictor and discard the step unless a slightly smaller step size, viz. $`\mathrm{\Delta }Qp_{\mathrm{discard}}^{[\mathrm{\Delta }Q]}`$ ($`p_{\mathrm{discard}}^{[\mathrm{\Delta }Q]}<1`$), passes the tests. If no step size can be found satisfying all the constraints, the calculation is terminated.

While the above set of prescriptions for finding suitable node locations has proved indispensible in understanding the behavior of the pde’s solution, the oscillatory nature of $`f(Q,\varrho )`$ invariably linked to the build-up of the isothermal compressibility’s divergence for sub-critical temperatures prevent its use for $`\beta >\beta _c`$: considering even the modest value $`f10^3`$, $`\sigma \mathrm{\Delta }Q`$ would have to be smaller than $`e^{10^3}10^{430}`$, which is obviously completely useless for any practical implementation. Thus, even though it means loosing control over the level of accuracy in the solution, we have also implemented a version of this main part with predetermined step sizes that just happen to often be sufficient for reaching $`Q=Q_0`$ even well below the critical temperature while reproducing the overflow necessary for $`\kappa _T`$’s divergence in a density interval the edges of which may then be identified with the coexisting phases’ densities $`\varrho _v`$ and $`\varrho _l`$. Recalling the behavior of $`f`$ wherever it is large we obviously have to drastically reduce $`\mathrm{\Delta }Q`$ as we approach $`Q_0`$; for this we use the very prescription introduced by the authors of  and evidently underlying all later published hrt calculations.

One last aspect of this main part to be mentioned regards the choice of $`Q_{\mathrm{}}`$ in both implementations: As the only reasonable initial condition for the core condition assumes that the structure at $`Q_{\mathrm{}}`$ is basically the same as that for $`Q=\mathrm{}`$ (so that $`c_2^{(Q_{\mathrm{}})}=c_2^{(\mathrm{})}=c_2^{\mathrm{ref}}`$ or, equivalently, $`\gamma _n^{(Q_{\mathrm{}})}=\gamma _n^{(\mathrm{})}=0`$, $`n0`$) and the same set of parameters must also be used for the initialization of nodes at $`Q=Q_{\mathrm{}}+|\mathrm{\Delta }Q|`$ (the nodes labled $`(1i)`$ in fig. 1 in the first step), it is preferable to have $`\gamma _n^{(Q)}(\varrho )/Q=0`$ at $`Q=Q_{\mathrm{}}`$; from eq. (11) one immediately concludes that this is equivalent to $`\stackrel{~}{w}(Q_{\mathrm{}})=0`$ whenever using the decoupling assumption. It is left to main part “ansatz” to decide whether $`Q_{\mathrm{}}`$ should be determined in this way, thereby necessarily introducing $`\varrho `$-dependent $`Q_{\mathrm{}}`$ when dealing with a $`\varrho `$-dependent potential.

## IV Discussion and conclusion

In the preceding sections we had to introduce a number of approximations, some of which may seem rather less justified; their respective importance for and bearing on our program’s predictions of structure and thermodynamics of simple liquids now remains to be assessed. As an exact solution with which to compare numerical results is lacking, for hrt as implemented by our software package to be considered a reliable tool well applicable to realistic physical potentials along the lines outlined so far it is necessary to demonstrate the limited effects variations in the numerical recipe have and to compare the results obtained with those available by other means for certain potentials. For this contribution we turn to the hard-core Yukawa fluid with $`z=1.8/\sigma `$ for illustration, a system that has also been considered in a recent study  comparing various thermodynamically consistent approaches including hrt to short-ranged potentials.

Among the approximations just mentioned, an important one is the necessary truncation of eq. (11) to a finite number $`N_{\mathrm{cc}}+1`$ of basis functions $`u_n`$ and expansion coefficients $`\gamma _n^{(Q)}`$, $`n=0,\mathrm{},N_{\mathrm{cc}}`$; in fact, not only must $`N_{\mathrm{cc}}`$ be finite, it should also be rather small if evaluation of the slowly-convergent $`\widehat{}^{(Q)}`$-integrals at $`Q=Q_{\mathrm{}}`$ is not to dominate program execution time. An obvious test for the minimum number of basis functions to keep is to look at the $`N_{\mathrm{cc}}`$-dependence of the phase behavior predicted as summarized in table I, where the inverse critical temperature $`\beta _c`$, the critical density $`\varrho _c`$, and the coexisting densities at $`\beta =0.9/ϵ`$ are recorded; the latter temperature is sufficiently far away from the critical point so that the differences in $`\varrho _v`$ and $`\varrho _l`$ are not merely to be attributed to the differences in $`\beta _c`$ while maintaining sufficient separation of the binodal from the boundaries at $`\varrho _{\mathrm{min}}=0`$ and $`\varrho _{\mathrm{max}}`$ (v. i.). From this we find inclusion of the core condition to be of vital importance in determining the fluid’s phase behavior, while non-negligible variation of the results remains even for the highest $`N_{\mathrm{cc}}`$-values considered; however, the amount of variation especially in $`\beta _c`$ drops significantly as soon as $`N_{\mathrm{cc}}`$ exceeds 5. — The real test for the use of the truncated eq. (11) and the ad hoc approximation (12) (easily demonstrated to be inadequate at least for certain potentials ), is the pair distribution function $`g^{(Q_0)}(r)`$ as obtained from the final values of the expansion coefficients $`\gamma _n^{(Q_0)}`$ by performing the inverse Fourier transformation implied by the Ornstein-Zernike equation (5): Employing the odes following from the consistent application of the decoupling assumption in order to isolate the approximations pertaining to the implementation of the core condition from other effects (v. i.), the $`g^{(Q_0)}(r)`$ so obtained generally takes on rather large values for small $`r`$ while remaining within a few percent of the contact value $`g^{(Q_0)}(\sigma +)`$ for larger $`r`$ up to $`\sigma `$; while increasing $`N_{\mathrm{cc}}`$ usually does not considerably reduce the magnitude of $`g^{(Q_0)}(r)`$ for $`r`$ close to $`0`$, it instead extends the $`r`$-range of rather small $`g^{(Q_0)}(r)`$ to ever smaller $`r`$. At high density there is no substantial improvement in $`g^{(Q_0)}(r,\varrho )`$, $`r<\sigma (\varrho )`$, for $`N_{\mathrm{cc}}>5`$, nor for $`N_{\mathrm{cc}}>7`$ at low density, a finding corroborated by direct inspection of the final values of the expansion coefficients $`\gamma _n^{(Q_0)}`$ . Accordingly, $`N_{\mathrm{cc}}`$ should probably be chosen no less than 7 (corresponding to a $`6^{\mathrm{th}}`$ order polynomial in $`r`$ for $`𝒞^{(Q)}(r)`$ inside the core), whereas $`N_{\mathrm{cc}}=5`$ may still be sufficiently accurate for some applications; for $`N_{\mathrm{cc}}<5`$, on the other hand, we cannot expect significantly better results than those from consistently solving the pde without implementing the core condition (cf. sub-section III D), which, of course, runs much faster and at least does not rely on inconsistent approximations. But note that the core condition is always but poorly met, irrespective of $`N_{\mathrm{cc}}`$, whenever $`f(Q_0,\varrho )`$ is large (corresponding to the coexistence region or the critical point’s vicinity in implementations relying on a pde). On the other hand, when solving the pde without implementing the core condition at all, $`g^{(Q_0)}(r)`$, $`r<\sigma `$, can, of course, become arbitrarily large; e. g. for $`\beta =0.7/ϵ`$ and $`\varrho =0.9/\sigma ^3`$, $`g^{(Q_0)}(r)=3.26`$ inside the core while the contact value is $`g^{(Q_0)}(\sigma +)=+1.91`$. All in all, while systematic shortcomings in the pair distribution function $`g^{(Q_0)}(r)`$ itself cannot be avoided in an implementation relying on eqs. (11) and (12) with finite $`N_{\mathrm{cc}}`$, we needs must keep the core condition in the calculation due to its bearing on the phase behavior predicted; this is somewhat at variance with earlier findings  indicating only a modest influence of the core condition upon the results, a finding expressly referred to in .

On the other hand, as pointed out already in section II, retaining the core condition is possible only when adopting the decoupling assumption, viz. $`\alpha ^{(Q)}(\varrho )=0`$, an assumption that itself suffices to decouple the pde to a set of odes (which is why we can use $`\alpha ^{(Q)}(\varrho _{\mathrm{max}})=0`$ as a boundary condition). But if such a procedure is to be considered harmless, the results of the consistent application of this approximation to the closure (8) should not differ much from those of the pde applying the decoupling assumption only in order to get rid of the second $`\widehat{}`$-integral in eq. (11). However, the calculations summarized in fig. 2 clearly show that the two approaches yield very different results so that we cannot rule out a non-negligible effect on the structural and thermodynamic properties predicted: Most importantly, the odes cannot reproduce well-defined phase boundaries, and they even yield slightly negative inverse compressibility $`1/\kappa _T`$ in what would otherwise be the coexistence region. Preserving the structure of the pde, so that thermodynamic consistency is at least partly implemented by the pde’s coefficients $`d_{0i}`$ as defined in eq. (A1), seems sufficient to remedy these deficiencies; at any rate, we have to accept the decoupling assumption as indispensable for the implementation of the core condition.

Let us now turn to a brief discussion of several other aspects of hrt in its present formulation and the potential problems they may present; as most of these effects are much more prominent for the square-well and multistep potentials defined in sub-section III B, we relegate more thorough treatment to . — Clearly, the importance of retaining the structure of a pde mandates closer examination of the properties of the set of densities, especially since the terms corresponding to the finite difference approximation to the operator $`(^2/\varrho ^2)`$ have been found the primary limiting factor for the quality of the numerical solution (cf. the definition of $`ϵ_\mathrm{\#}`$ in sub-section III E); while obvious for subcritical temperatures at low $`Q`$ where the near-discontinuity of $`f`$ at the phase boundary betrays the smoothness assumptions underlying finite difference schemes, this is true even for rather small $`f(Q,\varrho )`$. For our hard-core Yukawa system, however, the results’ stability with respect to a variation of the density grid or the location and nature of the high density boundary (cf. table II) is rather satisfactory as long as $`\varrho _{\mathrm{max}}\varrho _l`$ exceeds several $`\mathrm{\Delta }\varrho `$: just as expected from fig. 2, the ode used at the boundary forces the corresponding density to lie outside the coexistence region, which carries over to nearby densities by virtue of the pde’s discretization. Despite the identical phase behavior found for different boundary conditions as evidenced by table II, the question of which condition to impose at $`\varrho _{\mathrm{max}}`$ is far from irrelevant; indeed, when imposing the loga/orpa-condition $`\gamma _0^{(Q)}=0`$ at $`\varrho _{\mathrm{max}}`$ without making use of decoupling at all, inspection of the solution generated close to $`\varrho _{\mathrm{max}}`$ clearly shows the need to replace the original implementation’s approach by an ansatz less inconsistent. — Similarly, the stiffness of the pde for subcritical temperatures shortly touched upon in section II is easily detected by employing the criteria discussed in sub-section III E; still, for $`v^{\mathrm{hcy}}`$ with $`z=1.8/\sigma `$ and $`N_{\mathrm{cc}}=7`$, using pre-determined step sizes and resigning on any control of the predictor-corrector scheme’s convergence we obtain an approximate solution at $`Q_0`$ sufficiently stable outside the coexistence region; thus, while the solution generated in the region of large $`f(Q,\varrho )`$ necessarily differs from the true solution in a fundamental way, the influence of which on the data produced outside the coexistence region is not assessed easily, stiffness here appears less pressing a concern than for other systems .

In addition to the internal consistency of the results, we also have to make contact with data available by other means; as this paper’s purpose is presentation of our software and discussion of some general aspects of hrt rather than a comprehensive study of the hard-core Yukawa system, we here restrict ourselves to the data of table I of  containing the predictions of various thermodynamically self-consistent liquid state theories including hrt in the original implementation, which are found to yield $`T_c`$ ranging from $`1.193ϵ/k_B`$ to $`1.219ϵ/k_B`$ ($`\beta _c=1/k_BT_c`$ between $`0.820/ϵ`$ and $`0.838/ϵ`$), as well as the Monte Carlo (mc) result $`T_c=1.212(2)ϵ/k_B`$ ($`\beta _c=0.825(2)/ϵ`$). As is apparent from our table I, our implementation’s predictions for $`1N_{\mathrm{cc}}4`$ fall precisely into this range and, for $`N_{\mathrm{cc}}\{3,4\}`$, are well compatible with simulation results; that very same $`N_{\mathrm{cc}}`$-range, on the other hand, we have seen (v. s.) is characterized by gross violation of the core condition due to an insufficient number of basis functions retained in the truncated eq. (11). As we further increase $`N_{\mathrm{cc}}`$ so that the core condition is obeyed to a certain extent (v. s.), however, $`\beta _c`$ drops dramatically to values far outside the range quoted in ; while the trend of increasing $`\beta _c`$ evident for $`N_{\mathrm{cc}}7`$ indicates that hrt might match the mc predictions for $`N_{\mathrm{cc}}15`$, we have not performed these cpu intensive calculations.

With this we conclude our superficial sketch of the software we have written and the appraisal of the results it typically produces as illustrated for the hard-core Yukawa potential $`v^{\mathrm{hcy}}`$ with $`z=1.8/\sigma `$: it should be apparent that hrt in its current formulation, while presenting substantial difficulties discussed here as well as in , is well capable of predicting structural and thermodynamic properties of simple one-component fluids; at the same time, the computational difficulties mentioned and the approximations introduced to render the numerics tractable cannot always be shown to be harmless so that great care has to be exercised in the interpretation of isolated results. The fully modular design of our program and the high degree of flexibility brought about by the adoption of a meta language and code construction techniques are key factors in facilitating additional evaluations, providing a means of separating different approximations’ effects; at the same time, our implementation makes for a versatile tool for the systematic exploration of hrt for one-component fluids in its present or alternative formulations.

## V Acknowledgments

ar thanks G. Stell (Stony Brook, N.Y.) and G. Pastore (Trieste) for stimulating discussions, as well as D. Pini (now in Milan) for sharing his experiences with HRT and the code of the original HRT program. The authors gratefully acknowledge support from Österreichischer Forschungsfonds under project number P13062-TPH.

## A Re-formulation in not-quite quasi-linear form

As noted in section II, it is advantageous to replace the highly non-linear pde in the modified free energy $`𝒜^{(Q)}`$ implied by eqs. (6), (7), and the core condition, with a formulation akin to a quasi-linear one. In complete analogy to , we define the auxiliary function $`f(Q,\varrho )`$ via

$$\mathrm{ln}\left(1\frac{\stackrel{~}{\varphi }(Q,\varrho )}{\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )}\right)=f(Q,\varrho )\stackrel{~}{u}_0^2(Q,\varrho )\frac{\stackrel{~}{\varphi }(Q,\varrho )}{\stackrel{~}{𝒦}^{(Q)}(Q,\varrho )};$$

in the ideal gas limit, i. e. for $`\varrho 0`$, we immediately find from the divergence of the terms $`1/\varrho `$ in $`\stackrel{~}{c}_2^{\mathrm{ref}}`$ and $`\stackrel{~}{𝒞}^{(Q)}`$ that $`f(Q,0)`$ = $`f(Q,0)/\varrho =0`$, which is a convenient boundary condition most implementations rely on. In the above definition we have taken advantage of some freedom regarding the $`f`$-term to reduce the number of floating-point multiplications; as $`u_0`$ is usually chosen to be strictly proportional to $`\varphi `$, our $`f`$ differs only by a constant factor from the choice of  which can be recovered by replacing the factor $`\stackrel{~}{u}_0^2(Q,\varrho )`$ by $`\stackrel{~}{\varphi }^2(Q,\varrho )`$; note, however, that the expressions we give here remain valid for a slightly more general choice of $`u_0`$, allowing the proportionality of the basis function and the perturbational part of the potential to hold only up to order $`\stackrel{~}{\varphi }^2`$ — a freedom that might be exploited to use more appropriate basis functions, giving $`𝒞^{(Q)}`$ a larger range in $`r`$-space. Also we should point out that the condition of non-singular coefficients in the pde (14), the very reason for the restriction on the relation between $`u_0`$ and $`\varphi `$ just mentioned as well as for the introduction of the term involving $`\stackrel{~}{\varphi }/\stackrel{~}{𝒦}^{(Q)}`$, only fixes a minimum exponent for the $`\stackrel{~}{u}_0`$-factor with which to multiply $`f`$ .

Inserting the above definition for $`f`$ into the relevant equations of section II and eliminating the expansion coefficient $`\gamma _0^{(Q)}`$ via the consistency condition (7), we can easily re-cast the original pde in the form of eq. (14). Dropping the obvious arguments and with the shorthand notations

$$\begin{array}{c}\epsilon =1\frac{\stackrel{~}{\varphi }}{\stackrel{~}{𝒞}^{(Q)}}=e^{f\stackrel{~}{u}_0^2+x_\varphi }\overline{\epsilon }=\epsilon 1\\ x_\varphi =\stackrel{~}{\varphi }/\stackrel{~}{𝒦}^{(Q)}\stackrel{~}{\varphi }_0=\stackrel{~}{\varphi }(0,\varrho )\stackrel{~}{𝒢}_0=\stackrel{~}{𝒢}^{(Q)}(0,\varrho )\end{array}$$

the coefficients of eq. (14) can be written as

$$\begin{array}{cc}\hfill d_{00}=& +\frac{\stackrel{~}{\varphi }}{Q}\left(\frac{\stackrel{~}{\varphi }_0^2}{\stackrel{~}{𝒦}^{(Q)}\stackrel{~}{\varphi }^2}\frac{\stackrel{~}{𝒦}^{(Q)}\overline{\epsilon }^2\stackrel{~}{\varphi }_0^2}{\epsilon \stackrel{~}{\varphi }^4}\frac{2f}{\stackrel{~}{\varphi }}\right)\hfill \\ & +\frac{\stackrel{~}{𝒦}^{(Q)}}{Q}\left(\frac{\overline{\epsilon }^2\stackrel{~}{\varphi }_0^2}{\epsilon \stackrel{~}{\varphi }^3}\frac{\stackrel{~}{\varphi }_0^2}{\left(\stackrel{~}{𝒦}^{(Q)}\right)^2\stackrel{~}{\varphi }}\right)\hfill \\ & \frac{^2\stackrel{~}{u}_0^2}{\varrho ^2}\frac{Q^2\overline{\epsilon }^2f\stackrel{~}{\varphi }_0}{4\pi ^2\epsilon \stackrel{~}{\varphi }^2}\hfill \\ & \frac{^2x_\varphi }{\varrho ^2}\frac{Q^2\overline{\epsilon }^2\stackrel{~}{\varphi }_0}{4\pi ^2\epsilon \stackrel{~}{\varphi }^2}\hfill \\ & \frac{\stackrel{~}{𝒢}_0^{(Q)}}{Q}\frac{\overline{\epsilon }^2\stackrel{~}{\varphi }_0}{\epsilon \stackrel{~}{\varphi }^2}\hfill \\ \hfill d_{01}=& \frac{\stackrel{~}{u}_0}{\varrho }\frac{Q^2\overline{\epsilon }^2}{\pi ^2\epsilon \stackrel{~}{\varphi }}\hfill \\ \hfill d_{02}=& \frac{Q^2\overline{\epsilon }^2}{4\pi ^2\epsilon \stackrel{~}{\varphi }_0}.\hfill \end{array}$$
(A1)
When evaluating these expressions care has to be taken whenever $`\stackrel{~}{\varphi }(Q,\varrho )`$ is close to zero as all terms of orders $`1/\stackrel{~}{\varphi }^2(Q,\varrho )`$ and $`1/\stackrel{~}{\varphi }(Q,\varrho )`$ must cancel; noting that $`(\overline{\epsilon }^2\stackrel{~}{\varphi }_0^2/\epsilon \stackrel{~}{\varphi }^3)(\stackrel{~}{\varphi }_0^2/(\stackrel{~}{𝒦}^{(Q)})^2\stackrel{~}{\varphi })`$, the coefficient of $`\stackrel{~}{𝒦}^{(Q)}/Q`$ in $`d_{00}`$, can be written in terms of the $`(\varphi /Q)`$-coefficient, in our numerical work the calculation of the terms affected by the cancellation for small $`\stackrel{~}{\varphi }(Q,\varrho )`$ proceeds via application of a fifth-order Taylor expansion of $`(\stackrel{~}{\varphi }_0^2/\stackrel{~}{𝒦}^{(Q)}\stackrel{~}{\varphi }^2)(\stackrel{~}{𝒦}^{(Q)}\overline{\epsilon }^2\stackrel{~}{\varphi }_0^2/\epsilon \stackrel{~}{\varphi }^4)(2f/\stackrel{~}{\varphi })`$; one should note that, even though the criterion for switching between the full analytic expressions and said expansion depends on $`ϵ_\mathrm{\#}`$, the order of this expansion (and a similar one for $`\overline{\epsilon }/\stackrel{~}{\varphi }`$) is not increased for very small values of $`ϵ_\mathrm{\#}`$, which is one of the few hard-coded limitations of our program.

From the given expressions for the $`d_{0i}`$ two more aspects are obvious: (i) there is substantial further simplification for a $`\varrho `$-independent potential (as $`u_0`$ depends upon the density only through the perturbational part $`\varphi `$ of the potential, and the basis functions entering $`𝒢^{(Q)}`$ only through the reference system’s hard-core radius $`\sigma (\varrho )`$), and (ii) the pde does not fall into the class of quasi-linear pdes due to the presence of the second $`\varrho `$-derivative of $`x_\varphi `$ in the coefficient $`d_{00}`$.

## B Basis functions for the core condition

While the basis function $`u_0(r,\varrho )w(r,\varrho )`$ in the closure (8) is fixed by the potential used, in principle there is ample freedom in choosing the $`u_n(r)`$, $`n1`$. Of course, when truncating eq. (11) after $`1+N_{\mathrm{cc}}`$ basis functions it is natural to ask for the set $`\{u_1(r),\mathrm{},u_{N_{\mathrm{cc}}}(r)\}`$ to span the space of polynomials of order up to $`N_{\mathrm{cc}}1`$, so that $`u_n(r)`$ will generally be a polynomial or order $`n1`$ in $`r`$; but whereas different polynomials of this type do not alter the function space, their choice has implications for the numerical properties of the matrix equations implementing the core condition as well as, to a certain extent, for the convergence of the $`\widehat{}`$-integrals to be evaluated at $`Q=Q_{\mathrm{}}`$ (cf. sub-section III D). Other than the original implementation that relied on an affine transformation of the Legendre polynomials, we found it convenient to choose $`u_n(r,\varrho )`$ simply proportional to a power of $`r`$ and normalized to $`\stackrel{~}{u}_n(0,\varrho )=1`$. Thus, dropping the obvious argument $`\varrho `$, we have

$$u_n(r)=\{\begin{array}{ccc}\frac{n+2}{4\pi \sigma ^3}\left(\frac{r}{\sigma }\right)^{n1}& :& r<\sigma \\ 0& :& r>\sigma ,\end{array}$$

the Fourier transform of which is

$$\stackrel{~}{u}_n(k)=\frac{n+2}{(\sigma k)^{n+2}}\left[n!\mathrm{cos}\frac{n\pi }{2}\underset{j=0}{\overset{n}{}}\frac{n!}{(nj)!}(\sigma k)^{nj}\mathrm{cos}\left(\sigma k+\frac{j\pi }{2}\right)\right],$$

an expression used for $`\sigma k>n`$ only for numerical reasons. For smaller $`k`$ we rely on the expansion

$$\stackrel{~}{u}_n(k)=(n+2)\underset{j=0}{\overset{\mathrm{}}{}}\frac{(1)^j}{(n+2j+2)(2j+1)!}(\sigma k)^{2j},$$

truncating the series after $`N_n`$ terms, where $`N_n`$ is the smallest number such that

$$\frac{n+2}{n+2N_n+4}\frac{n^{2N_n+2}}{(2N_n+3)!}ϵ_\mathrm{\#}p_{N_n}^{[u_n]}$$

with a customization factor $`p_{N_n}^{[u_n]}`$ of order unity.


text:
# A Nonlinear Realized Approach of SU(2) Chiral Symmetry Spontaneous Breaking and Properties of Nuclear Matter

## Abstract

A nonlinear realization of SU(2) chiral symmetry spontaneous breaking approach is developed in the composite operator formalism. A Lagrangian including quark, gluon and Goldstone boson degrees of freedom of the chiral quark model is obtained from the QCD Lagrangian. A way to link the chiral symmetry spontaneous breaking formalism at hadron level and that at quark level is predicted. too. The application to nuclear matter shows that the approach is quite successful in describing the properties of nuclear matter and the quark condensate in it.

PACS numbers: 11.30.Rd, 11.30.Qc, 12.39.Fe, 21.65.+f

It has been known that chiral symmetry and its spontaneous breaking is the key points in understanding many features of the nature, such as the generation of some particles in strong interaction physics, the chiral property of many amino acids, and so on. Moreover, it has been widely used in guiding the design and synthesis of chemical compounds to control and improve their functions. In low energy strong interaction physics, Goldstone bosons appear as the chiral symmetry is spontaneously broken. To carry out the constraint by the appearance of Goldstone bosons, several realization formalisms of the chiral symmetry spontaneous breaking ($`\chi `$SB) have been developed . Furthermore, Considering the effective degrees of freedom to describing hadron structure, we know that several formalisms have been developed. For example, the one with constituent quarks and Goldstone bosons has been shown to be successful in describing some properties of hadrons and hadron spectroscopy, even nucleon-nucleon interactions. The one including constituent quarks and gluons are very powerful to describe hadron properties, too. Moreover, Manohar and Georgi proposed that, besides the constituent quarks, both the Goldstone bosons and gluons are necessary to describe hadron structure. Such a scheme has also been widely used to investigate hadron spectroscopy and hadron interactions(see, for example, Refs. ). However, which effective degrees of freedom have more sophisticated QCD foundation is still a significant topic in hadron physics and low energy QCD. In the sprit of $`\chi `$SB and the composite operator scheme of the QCD, we propose, in this letter, a new approach to realize the $`\chi `$SB nonlinearly not only in formalism but also in practical calculation and obtain an effective Lagrangian including constituent quarks, Goldstone bosons and gluons. Then some light is shed on the effective degrees of freedom of hadron structure. We propose also a link between the formalism of $`\chi `$SB at hadron level and that at quark level. As an application, we take the approach to describe the properties of nuclear matter. By implementing the Hellmann-Feynman theorem, we evaluate the quark condensate in nuclear matter in the present formalism, too.

At first, we discuss the linear infinitesimal transformation in general

$$\varphi _n^{}(x)=\varphi _n+iϵ\underset{m}{}t_{nm}^\alpha \varphi _m(x)$$
(1)
where $`t^\alpha `$ is a generator of the symmetry group of the Lagrangian and $`\varphi _n`$ is a spin-zero boson field or spin-zero composite operators. The quantum effective action is defined by

$$\mathrm{\Gamma }[\varphi ]d^3x\underset{n}{}\varphi ^nJ_{\varphi _n}(x)+W[J_\varphi ]$$
(2)
This quantum effective action has the same symmetric property as the Lagrangian. It gives that

$$\underset{n,m}{}\frac{\delta \mathrm{\Gamma }[\varphi ]}{\delta \varphi _n}t_{nm}^\alpha \varphi _m(x)d^4x=0$$
(3)
Introducing an effective potential $`V(\varphi )`$ as $`\mathrm{\Gamma }[\varphi ]=vV[\varphi ]`$, where $`v`$ is the space-time volume, one can rewrite Eq.(3) as

$$\underset{n,m}{}\frac{V(\varphi )}{\varphi _n}t_{nm}^\alpha \varphi _m=0,$$
(4)
where $`\varphi `$ is independent of $`x`$. Taking the differentiation with respect to $`\varphi _l`$, we obtain

$$\underset{n}{}\frac{V(\varphi )}{\varphi _n}t_{nl}^\alpha +\underset{n,m}{}\frac{^2V(\varphi )}{\varphi _n\varphi _l}t_{mn}^\alpha \varphi _m=0$$
(5)
The symmetry spontaneous breaking appears when $`\frac{V(\overline{\varphi })}{\overline{\varphi }_m}=0`$ and $`\overline{\varphi }_m0`$ where $`\overline{\varphi }_m=0|\varphi _m|0`$. In accordance with the condition of the symmetry spontaneous breaking, Eq.(5) can be written as

$$\underset{n,m}{}\left(\frac{^2V(\varphi )}{\varphi _n\varphi _l}\right)_{\varphi =\overline{\varphi }}t_{nm}^\alpha \overline{\varphi }_m=0.$$
(6)
It is apparent that, if the symmetry is broken, $`_mt_{nm}^\alpha \overline{\varphi }_m`$ should not identically equate to zero. The massless eigenvectors of the mass matrix $`\frac{^2V(\varphi )}{\varphi _n\varphi _l}`$ span then a linear space in which the Goldstone bosons lie and the remainder has to be perpendicular to it. This constraint can lead various forms of the realization of the $`\chi `$SB $`SU_L(2)\times SU_R(2)SU(2)`$. In this paper we shall discuss the nonlinear realization of $`SU_L(2)\times SU_R(2)SU(2)`$ on the quark level at first. Then we propose a link between the $`\chi `$SB at quark level and that at hadron level.

Constructing four composite operators in terms of the current quark fields $`u`$ and $`d`$ as

$$\psi _4=\overline{q}q,\psi _i=i\overline{q}\gamma _5\tau ^iq$$
(7)
where $`\tau ^i`$ is the Pauli matrice, we can verify that the transformation of $`\psi _n`$(n=1,2,3,4) under the chiral symmetry transformation of the quarks $`e^{i\gamma _5ϵ^i\tau ^i}`$ can be written as

$$\delta \psi _4=2ϵ^i\psi _i,\delta \psi _i=2ϵ^i\psi _4$$
(8)
where $`ϵ^i`$ are the infinitesimal parameters.

Analyses in QCD sum rules and calculations in lattice QCD have shown that the vacuum expectation value of quarks (quark condensate) is not zero, i.e., $`0\left|\overline{q}(x)q(x)\right|00`$. It indicates that the chiral symmetry $`SU_L(2)\times SU_R(2)`$ is spontaneously broken and chiral condensate appears. As a consequence, the $`\psi _n`$ has to be separated into two parts: one contains Goldstone bosons and the other does not. This can be done by writing

$$\psi _n=\underset{m}{}(e^{2i\xi ^iT^{i4}})_{nm}\stackrel{~}{\psi }_m$$
(9)
where $`T^{\alpha \beta }`$ are the generators of the four dimensional rotation group. To make the $`\stackrel{~}{\psi }_n`$ do not contain Goldstone bosons, the $`\stackrel{~}{\psi }_n`$ can be taken as $`\stackrel{~}{\psi }_n=(\stackrel{~}{\psi }_4,0,0,0)`$. In the infinitesimal form, we have then

$$\delta \psi _i=2\xi ^i\stackrel{~}{\psi }_4$$
(10)
Comparing Eq.(10) with Eq.(8), we obtain the relation between the constituent quarks and the current quarks as

$$\stackrel{~}{q}=e^{i\gamma _5\tau ^i\xi ^i}q$$
(11)
In order to simplify the factor $`e^{i\gamma _5\xi ^i\tau ^i}`$ we introduce the Glodstone bosons as

$$\xi ^i=(\mathrm{tan}^1\left|\stackrel{}{\zeta }\right|)\frac{\zeta ^i}{\left|\stackrel{}{\zeta }\right|}$$
(12)
Using Eq.(12), we can obtain the nonlinear realization of the $`\chi `$SB as

$$e^{i\gamma _5(\mathrm{tan}^1\left|\stackrel{}{\zeta }\right|)\stackrel{}{\tau }\stackrel{}{\zeta }/|\stackrel{}{\zeta }|}=\frac{1i\gamma _5\stackrel{}{\tau }\stackrel{}{\zeta }}{\sqrt{1+\stackrel{}{\zeta }^2}}.$$
(13)
Since the theory leaves the SU(2) invariant, the $`\stackrel{~}{q}`$ constructs a representation of SU(2). Then it is necessary to have

$$e^{i\gamma _5\stackrel{}{ϵ}\stackrel{}{\tau }}e^{i\gamma _5(\mathrm{tan}^1|\stackrel{}{\zeta }|)\stackrel{}{\tau }\stackrel{}{\zeta }/|\stackrel{}{\zeta }|}=e^{i\gamma _5(\mathrm{tan}^1|\stackrel{}{\zeta }^{}|)\stackrel{}{\tau }\stackrel{}{\zeta }^{}/|\stackrel{}{\zeta }^{}|}e^{i\stackrel{}{\tau }\stackrel{}{\theta }}$$
(14)
After some algebraic calculation we obtain the relation between the Goldstone bosons and the constituent quarks as

$$\begin{array}{c}\text{ }\delta \stackrel{~}{q}=i\stackrel{}{\tau }(\stackrel{}{ϵ}\times \stackrel{}{\zeta }),(15a)\text{ }\hfill \\ \text{ }\stackrel{}{\theta }=\stackrel{}{ϵ}\times \stackrel{}{\zeta },(15b)\text{ }\hfill \\ \text{ }\delta \stackrel{}{\zeta }=\stackrel{}{ϵ}(1\stackrel{}{\zeta }^2)+2\stackrel{}{\zeta }(\stackrel{}{ϵ}\stackrel{}{\zeta }),(15c)\text{ }\hfill \end{array}$$

Under these transformations, the QCD Lagrangian can be effectively expressed in terms of the constituent quark, gluon and Goldstone bosons as

$$\begin{array}{c}\text{ }_{eff}=\underset{\stackrel{~}{q}}{}i\overline{\stackrel{~}{q}}\gamma ^\mu \left(_\mu +\frac{i(\stackrel{}{\zeta }\times _\mu \stackrel{}{\zeta })\stackrel{}{\tau }}{1+\stackrel{}{\zeta }^2}\frac{i\gamma _5\stackrel{}{\tau }_\mu \stackrel{}{\zeta }}{1+\stackrel{}{\zeta }^2}\right)\stackrel{~}{q}+i\overline{\stackrel{~}{q}}\gamma ^\mu \left(ig\frac{\lambda ^a}{2}A_\mu ^{a,p}\right)\stackrel{~}{q}\text{ }\hfill \\ \text{ }m\overline{\stackrel{~}{q}}\stackrel{~}{q}\frac{1}{4}G_{\mu \nu }^{ap}G^{a\mu \nu p}+_M.(16)\text{ }\hfill \end{array}$$

where $`m`$ is the constituent quark mass resulting from quark condensates. $`_M`$ is the Lagrangian of the Goldstone boson fields $`\stackrel{}{\zeta }`$ and the $`\stackrel{~}{\psi }_4`$. In this sense, the effective degrees of freedom to describe hadron structure should involve constituent quarks, gluons and Goldstone bosons. One can then recognize that the constituent quark model of hadron structure involving constituent quarks, Goldstone bosons and gluons has a quite solid QCD foundation. Furthermore, the double counting and the spurious state involved in Manohar-Georgi’s formalism disappears in the present approach.

¿From the effective action $`\mathrm{\Gamma }`$ in Eq.(2) we have

$$_M=\frac{1}{2A}_\mu \stackrel{~}{\psi }_4^\mu \stackrel{~}{\psi }_4\frac{2}{A}\stackrel{~}{\psi }_4\stackrel{~}{\psi }_4\stackrel{}{D}_\mu \stackrel{}{D}^\mu +(\stackrel{~}{\psi }_4)$$
(17)
where $`D_\mu ^i=\frac{_\mu \zeta ^i}{1+\stackrel{}{\zeta }^2}`$. $`A`$ is a dimensional constant. Considering the pseudoscalar field property of $`\zeta _i`$ and the assumption of the saturation of vacuum, we obtain that the $`_M`$ can be written as the superposition of the following two parts

$$_\sigma =\frac{1}{2}_\mu \sigma ^\mu \sigma \frac{1}{2}m_\sigma ^2\sigma ^2\frac{g_2}{3}\sigma ^3\frac{g_3}{4}\sigma ^4,$$
(18)
$$_\pi =\frac{1}{2}\frac{^\mu \stackrel{}{\pi }}{1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}}\frac{^\mu \stackrel{}{\pi }}{1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}}2f_\pi \sigma \frac{^\mu \stackrel{}{\pi }^\mu \stackrel{}{\pi }}{\left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)^2}2\sigma ^2\frac{^\mu \stackrel{}{\pi }^\mu \stackrel{}{\pi }}{\left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)^2}.$$
(19)
And the parameters $`g_2`$ and $`g_3`$ hold relation $`g_2=\frac{3m_q\sqrt{A}}{m_\sigma ^2}g_3`$ with $`m_q`$ being the current quark mass.

It has been well known that there exists $`\chi `$SB formalism at hadron level. However the relation between the formalism and that at quark level is not clear enough. To link the $`\chi `$SB formalism at hadron level with the $`\chi `$SB formalism at quark level described above, we construct a set of composite operators which have the quantum numbers of proton and neutron as

$$\begin{array}{c}\text{ }p(x)=\epsilon _{abc}(u^{aT}(x),d^{aT}(x))c\gamma ^\mu u^b(x)\gamma _5\gamma _\mu i\tau _2\left(\begin{array}{c}u^c(x)\hfill \\ d^c(x)\hfill \end{array}\right),(20a)\text{ }\hfill \\ \text{ }n(x)=\epsilon _{abc}(u^{aT}(x),d^{aT}(x))c\gamma ^\mu d^b(x)\gamma _5\gamma _\mu i\tau _2\left(\begin{array}{c}u^c(x)\hfill \\ d^c(x)\hfill \end{array}\right),(20b)\text{ }\hfill \end{array}$$

where $`a`$, $`b`$, $`c`$ are the color indices, $`\tau _2`$ acts only on the isospin indices and the row and column matrix express the doublet of the isospin. Since the nucleon in this doublet are composed of current quarks, we can refer them as current nucleon. Implementing the chiral transformation of quarks $`\stackrel{~}{q}=\mathrm{exp}(i\gamma _5\xi ^i\tau ^i)q`$, we have the transition among nucleons as

$$\stackrel{~}{N}=e^{i\gamma _5\xi ^i\tau ^i}N,$$
(21)
where $`N=\left(\begin{array}{c}p(x)\hfill \\ n(x)\hfill \end{array}\right)`$ , $`\stackrel{~}{N}=\left(\begin{array}{c}\stackrel{~}{p}(x)\hfill \\ \stackrel{~}{n}(x)\hfill \end{array}\right)`$ in which

$$\begin{array}{c}\text{ }\stackrel{~}{p}(x)=\epsilon _{abc}(\stackrel{~}{u}^{aT}(x),\stackrel{~}{d}^{aT}(x))c\gamma ^\mu \stackrel{~}{u}^b(x)\gamma _5\gamma _\mu i\tau _2\left(\begin{array}{c}\stackrel{~}{u}^c(x)\hfill \\ \stackrel{~}{d}^c(x)\hfill \end{array}\right),(22a)\text{ }\hfill \\ \text{ }\stackrel{~}{n}(x)=\epsilon _{abc}(\stackrel{~}{u}^{aT}(x),\stackrel{~}{d}^{aT}(x))c\gamma ^\mu \stackrel{~}{d}^b(x)\gamma _\mu \gamma _5i\tau _2\left(\begin{array}{c}\stackrel{~}{u}^c(x)\hfill \\ \stackrel{~}{d}^c(x)\hfill \end{array}\right),(22b)\text{ }\hfill \end{array}$$

The Eq. (21) shows that we can obtain the nucleon from the current nucleon. The nucleon observed in experiments in the low energy region can be formed only through Eq. (22).

Taking the same measure as described for the quarks and considering the Goldberg-Treiman relation and the partial conservation of the axial current(PCAC), we obtain the Lagrangian density with the nonlinear realization of the SU(2) $`\chi `$SB at hadron level as

$`_{ch}`$ $`=`$ $`\overline{\psi }\left(i\gamma _\mu ^\mu M_N{\displaystyle \frac{\stackrel{}{\tau }\left(\stackrel{}{\pi }\times \text{/}\stackrel{}{\pi }\right)}{f_\pi ^2\left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)}}+{\displaystyle \frac{2iM_Ng_A\gamma _5\stackrel{}{\tau }\stackrel{}{\pi }}{f_\pi \left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)}}{\displaystyle \frac{2g_A\gamma _5\stackrel{}{\tau }\stackrel{}{\pi }\stackrel{}{\pi }\text{/}\stackrel{}{\pi }}{f_\pi ^3\left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)^2}}\right)\psi `$ (23)
$`\sigma ^2\overline{\psi }\left({\displaystyle \frac{8ig^{}M_N\gamma _5\stackrel{}{\tau }\stackrel{}{\pi }}{f_\pi \left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)}}+{\displaystyle \frac{8g^{}\gamma _5\stackrel{}{\tau }\stackrel{}{\pi }\stackrel{}{\pi }\text{/}\stackrel{}{\pi }}{f_\pi ^3\left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)^2}}\right)\psi {\displaystyle \frac{\delta M_p+\delta M_n}{2}}\left({\displaystyle \frac{1\frac{\stackrel{}{\pi }^2}{f_\pi ^2}}{1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}}}\right)\overline{\psi }\psi `$
$`\left(\delta M_p\delta M_n\right)\left({\displaystyle \frac{\tau _3}{2}}{\displaystyle \frac{1}{f_\pi ^2}}\left({\displaystyle \frac{\pi _3}{1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}}}\right)\stackrel{}{\tau }\stackrel{}{\pi }\right)\overline{\psi }\psi g_\sigma \overline{\psi }\psi \sigma `$
$`{\displaystyle \frac{1}{2}}_\mu \sigma ^\mu \sigma U(\sigma )\left({\displaystyle \frac{1}{2}}+2f_\pi \sigma +2\sigma ^2\right){\displaystyle \frac{_\mu \stackrel{}{\pi }^\mu \stackrel{}{\pi }}{\left(1+\frac{\stackrel{}{\pi }^2}{f_\pi ^2}\right)^2}}{\displaystyle \frac{1}{2}}m_\pi ^2\stackrel{}{\pi }^2,`$
where $`M_N`$ is the mass of the nucleon caused by the $`\chi `$SB, $`\delta M_p`$ and $`\delta M_n`$ are the masses of the current proton and neutron, respectively. The $`U(\sigma )`$ stands for

$$U\left(\sigma \right)=\frac{1}{2}m_\sigma ^2\sigma ^2+\frac{1}{3}g_2\sigma ^3+\frac{1}{4}g_3\sigma ^4,$$
(24)
where the $`g_2`$ and $`g_3`$ have relation $`g_2=\frac{3\delta M_N\sqrt{A}}{m_\sigma ^2}g_3`$ with $`\delta M_N=\delta M_n=\delta M_p`$ .

To represent the repulsive interaction among nucleons and the isospin symmetry breaking in nuclear matter, we can include $`\omega `$ and $`\rho `$ mesons in the way having been widely used. Then in the mean field approximation, we have the Lagrangian

$`_{RMF}`$ $`=`$ $`\overline{\psi }\left(i\gamma _\mu ^\mu M_N\right)\psi \delta M_N\left({\displaystyle \frac{1\frac{\pi _0^2}{f_\pi ^2}}{1+\frac{\pi _0^2}{f_\pi ^2}}}\right)\overline{\psi }\psi `$ (25)
$`g_\sigma \overline{\psi }\psi \sigma _0g_\omega \overline{\psi }\gamma ^0\psi \omega _0g_\rho \overline{\psi }\gamma ^0t_3\rho _{03}\psi `$
$`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma _0^2{\displaystyle \frac{g_2}{3}}\sigma _0^3{\displaystyle \frac{g_3}{4}}\sigma _0^4{\displaystyle \frac{1}{2}}m_\pi ^2\pi _0^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _0^2+{\displaystyle \frac{1}{2}}m_\rho ^2\rho _{03}^2,`$
where the $`\sigma _0`$ is the expectation value of the isoscalar scalar field, $`\pi _0^2`$ is that of $`\stackrel{}{\pi }^2`$, $`\omega _0`$ is that of the temporal component for the isoscalar vector field since there is no spatial direction for a uniform nuclear matter at rest, $`\rho _{03}`$ is that of $`\rho `$ meson in the nuclear matter.

By fitting the saturation properties of nuclear matter, we obtain the parameters with the restriction $`M_N+\delta M_N=938`$ MeV. Two sets of the parameters are listed in Table 1. For the parameter set $`C_1`$, we get the saturation density of $`0.152\text{fm}^3`$, binding energy of $`15.297`$ MeV, a compression modulus of $`349.10`$ MeV, symmetry energy coefficient $`33.645`$ MeV and the effective mass of a nucleon of $`0.687M_N`$ for symmetric nuclear matter. The parameter set $`C_2`$ gives a saturation density $`0.151\text{fm}^3`$, a binding energy of nucleon $`15.040`$MeV, a compression modulus $`K=298.88`$ MeV, a symmetry energy coefficient $`32.593`$ MeV and an effective mass of nucleon $`0.736M_N`$ for symmetric nuclear matter. Meanwhile the curves for the equation of states (EOS) are obtained too. The numerical results show that the EOSs for the two sets of parameters are quite close to each other. we display then only the equation of states for the parameter set $`C_2`$ in Fig. 1. The figure shows evidently that the symmetric nuclear matter can exist stably. However stable pure neutron matter does not exist.

Making use of the Hellmann-Feynman theorem, one can obtain the relation between the quark condensate in nuclear matter $`Q(\rho )=0\left|\overline{q}(x)q(x)\right|0_{\rho 0}`$ and that in vacuum $`Q(0)=0\left|\overline{q}(x)q(x)\right|0_{\rho =0}`$ as

$`{\displaystyle \frac{Q(\rho )}{Q(0)}}`$ $`=`$ $`1+{\displaystyle \frac{1}{2Q(0)}}\left({\displaystyle \underset{N=n,p}{}}{\displaystyle \frac{\epsilon }{M_N}}{\displaystyle \frac{dM_N}{dm_q}}+{\displaystyle \frac{\epsilon }{M_\sigma }}{\displaystyle \frac{dM_\sigma }{dm_q}}+{\displaystyle \frac{\epsilon }{M_\pi }}{\displaystyle \frac{dM_\pi }{dm_q}}\right)`$ (26)
$`=`$ $`1{\displaystyle \frac{2\sigma _N}{(m_\pi f_\pi )^2}}\left(\rho _s(p)+\rho _s(n)\right){\displaystyle \frac{\pi _0^2}{f_\pi ^2}}{\displaystyle \frac{3g_3m_\pi ^2\sigma _0^2}{2m_\sigma ^4}}`$
where $`\sigma _N`$ is the nucleon $`\sigma `$ term, $`\rho _s=\overline{\psi }\psi `$ is the scalar density of the nucleons.

With the parameter sets $`C_2`$ determined above and $`\sigma _N=45`$ MeV, we evaluate the ratio of the in-medium quark condensate to that in vacuum. The obtained results are represented in Fig. 2. From the figure one can easily know that the quark condensate in nuclear matter decreases, in general, as the nuclear matter density increases. It manifests that the “upturn” problem in the other approaches at hadron level is removed in the present $`\chi `$SB approach. Investigating the figure more carefully, one may know that there exists a density $`\rho =0.235\text{fm}^3`$ at which the decreasing rate is enhanced. Such a behavior is consistent with the result given in Ref.. Relevant investigation shows that pion condensate in strong interaction appears as the nuclear matter density reaches $`0.235\text{fm}^3`$ and beyond. It indicates that the sudden increase of the decreasing rate or the discontinuous decrease of the quark condensate in nuclear matter may come from the appearance of pion condensate.

In summary, a nonlinear realization of SU(2) $`\chi `$SB approach is developed in the composite operator formalism of QCD. A Lagrangian including quark, gluon and Goldstone bosons of chiral quark model is obtained from the QCD Lagrangian. It indicates that, besides the constituent quarks, both the gluons and Goldstone bosons are the effective degrees of freedom to describe hadrons. We constructed also a way link the $`\chi `$SB at hadron level and that at quark level. As an application, the properties of symmetric nuclear matter and pure neutron matter are investigated in the mean-field approximation. The calculation shows that the present approach is quite powerful in describing the properties of nuclear matter. Meanwhile the quark condensate is evaluated. It shows that the quark condensate decreases with the increasing of nuclear matter density. The chiral symmetry is then restored gradually in nuclear matter as the density increases.

This work is supported by the National Natural Science Foundation of China under contact No. 10047001, 10075002, 10135030, and 19975062, the Major State Basic Research Development Program under contract No.G2000077400 and CAS Knowledge Innovation Project No. KJCX2-N11. One of the author(Liu) thanks also supports by Foundation for University Key Teacher by the Ministry of Education, China. The authors are indebted to Professors F. Wang, C. S. Chang, W. Q. Chao, P. F. Zhuang and H. S. Zong for their helpful discussions.

Table 1. The parameters in the calculation of the SU(2) chiral symmetry spontaneous breaking model with $`m_\pi =139.57MeV`$ and $`f_\pi =130MeV`$ ($`m_i(i=\sigma ,\omega ,\rho )`$, $`M_N`$ and $`\delta M_N`$ in Mev, $`g_2`$ in fm<sup>-1</sup>).

|  | $`g_\sigma (N)`$ | $`m_\sigma `$ | $`g_\omega (N)`$ | $`m_\omega `$ | $`g_\rho (N)`$ | $`m_\rho `$ | $`g_2`$ | $`g_3`$ | $`M_N`$ | $`\delta M_N`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`C_1`$ | 9.111 | 540 | 10.587 | 783 | 8.480 | 770 | $`4.0`$ | 20.0 | 888 | 50 |
| $`C_2`$ | 9.111 | 570 | 9.573 | 783 | 8.480 | 770 | $`9.0`$ | 37.5 | 888 | 50 |

Figure Captions

Obtained average energy per nucleon $`\epsilon /\rho _NM_N`$ as a function of nucleon density $`\rho _N`$ for the parameter sets $`C_2`$. The solid line denotes symmetric nuclear matter , and the dashed line is for pure neutron matter.

Obtained ratio of the quark condensate in nuclear matter to that in vacuum as a function of nucleon density for the parameter sets $`C_2`$. The solid line denotes that in symmetric nuclear matter, and the dashed line is that for pure neutron matter.