text:
# The Hierarchical Reference Theory as applied to square well fluids of variable range

## I Introduction

In a large part of the density-temperature plane, integral equation theories are a reliable tool for studying thermodynamic and structural properties of, among others, simple one-component fluids ; unfortunately, in the vicinity of a liquid-vapor critical point, integral equations are haunted by a host of difficulties, leading to a variety of shortcomings such as incorrect and non-matching branches of the binodal, classical values at best for the critical exponents, or other deviations from the correct behavior at the critical singularity . Asymptotically close to the critical point, on the other hand, renormalization group (rg) theory is the instrument of choice for describing the fluid; in general, however, rg approaches do not allow one to derive non-universal quantities from microscopic information only, i. e. from knowledge of the forces acting between the fluid’s particles alone. One of the theories devised to bridge the conceptual gap between these complementary approaches is the Hierarchical Reference Theory (hrt) first put forward by Parola and Reatto : In this theory the introduction of a cut-off wavenumber $`Q`$ inspired by momentum space rg theory and, for every value of $`Q`$, of a renormalized potential $`v^{(Q)}(r)`$ means that only non-critical systems have to be considered at any stage of the calculation; consequently, integral equations may successfully be applied to every system with $`Q>0`$, and critical behavior characterized by non-classical critical exponents is recovered only in the limit $`Q0`$.

While applicability of hrt to a number of interesting systems, ranging from a lattice gas or Ising model  to various one-component fluids  even including three-body interactions , internal degrees of freedom , or non-hard-core reference systems , was demonstrated early on, the main focus of research on hrt has since shifted to the richer phase behavior of binary systems . Nevertheless, in the light of hrt’s high promise and low penetration into the liquid physics community, further study and critical assessment of this theory seem worthwhile, even and foremost in the case of simple one-component fluids: indeed, it is in this comparatively simple setting that we may gain important insights into the numerical side of the theory, and barring special mechanisms relevant to some specific model system only, any problems uncovered here must be expected to haunt more advanced applications of hrt, too. In our work we have found it convenient to restrict ourselves even further, implementing hrt in its usual formulation  for purely pairwise additive interactions via a potential $`v(r)`$ obtained from the superposition of an infinitely repulsive hard sphere serving as reference system, $`v^{\mathrm{ref}}(r)=v^{\mathrm{hs}}(r)`$, and a predominantly attractive tail $`w(r)`$, $`\stackrel{~}{w}(0)<0`$. Here we have made use of the notation introduced previously : superscripts always denote the system a quantitiy refers to (here, “ref” and “hs” for the reference system and hard spheres, respectively; similarly, “$`(Q)`$” for the system with cut-off $`Q`$), and a tilde indicates Fourier transformation. In the present contribution we apply our recent re-implementation of the theory  to one of the simplest potentials exhibiting phase separation, viz. the square well potential $`v^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}`$ (cf. sub-section III B of ref. ):

$$\begin{array}{c}v^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}(r)=v^{\mathrm{hs}[\sigma ]}(r)+w^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}(r)\\ v^{\mathrm{hs}[\sigma ]}(r)=\{\begin{array}{ccc}+\mathrm{}& :& r<\sigma \\ 0& :& r>\sigma \end{array}\\ w^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}(r)=\{\begin{array}{ccc}ϵ& :& r<\lambda \sigma \\ 0& :& r>\lambda \sigma .\end{array}\end{array}$$
(1)
Considering density-independent potentials only and chosing the hard core diameter $`\sigma `$ and the well’s depth $`ϵ`$ as units of length and energy, respectively, the attractive well’s range $`\lambda `$ is the sole remaining parameter; in this report we will study values of $`\lambda `$ from slightly above unity up to 3.6.

With just one parameter, viz. $`\lambda `$, to vary, square wells obviously make for a convenient test case of hrt and, indeed, of liquid state theories in general; consequently, a great many simulational and theoretical efforts have been directed at this system, and studies of its phase behavior abound . But square wells are also of interest in their own right, serving as — albeit somewhat crude — models of a wide variety of physical systems including, e. g., $`{}_{}{}^{3}\mathrm{He}`$, $`\mathrm{Ne}`$, $`\mathrm{Ar}`$, $`\mathrm{H}_2`$, $`\mathrm{CO}_2`$, $`\mathrm{CH}_4`$, $`\mathrm{C}_2\mathrm{H}_6`$, $`n`$-pentane and $`n`$-butane  while current interest in this potential derives mainly from the finding that square wells capture the essential features of the interactions found in colloidal systems . Yet another motivation for this first application of hrt to square wells comes from a recent, very accurate simulation study  confirming and quantifying the presence in the system with $`\lambda =1.5`$ of the Yang-Yang (YY) anomaly expected and experimentally found for asymmetric fluids .

Due to the extensive amount of data available in the literature the more recent of which will shortly be presented later on, and in view of some of the limitations of hrt in its current formulation we cannot expect to gain new insight into the system at hand with a level of precision comparable to that of the more sophisticated simulation schemes. Instead, in the present contribution our focus of interest lies on some aspects of hrt’s numerical side, specifically on those that are sensitive to the potential’s range: indeed, as stated already in ref. , for a potential as pronouncedly short-ranged as square wells some of the numerical problems should show up much more prominently than in other systems like, e. g., the hard-core Yukawa fluid previously considered  where they are, of course, in principle still present but do not manifest themselves as clearly.

In accordance with the preceding remarks, another reason why application of hrt to square wells might be worthwhile lies in the closure underlying seemingly all applications so far of hrt to simple one-component fluids with hard sphere reference part: As the usual formulation of hrt in these cases relies on an ansatz for the two-particle direct correlation function $`c_2(r)`$ very much in the spirit of Stell’s Lowest-Order $`\gamma `$-ordered Approximation  (loga) or the equivalent Optimized Random-Phase Approximation  (orpa) by Andersen and Chandler, the direct correlation function can never extend to larger $`r`$ values than the potential $`v(r)`$ itself. In particular, for the square well potential $`v^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}(r)`$ we necessarily have $`c_2(r)=0`$ for $`r>\lambda \sigma `$ so that all moments of $`c_2(r)`$, i. e. $`_^3\mathrm{d}^3rc_2(r)r^n`$, $`n0`$, exist, which is obviously at variance with the correct behavior near the critical point ; furthermore, at intermediate $`Q`$ the direct correlation function can hardly be considered satisfactory, especially  close to $`r=\lambda \sigma `$. While some earlier publications  already blamed unsatisfactory aspects of hrt results on this inadequacy of the closure, square wells should bring out related problems of hrt with the usual loga/orpa-style closure much more clearly, and the numerical procedure’s response should provide us with a signature to be looked out for in other systems, too; also, even within the loga/orpa-style approximation the implementation of the core condition via approximate ordinary differential equations (odes) for the relevant expansion coefficients easily shown to be inadequate for very short-ranged potentials  casts some doubt on the range of $`\lambda `$ values amenable to an hrt treatment in the current formulation of the theory. Determination of the admissible $`\lambda `$-range, on the other hand, is particularly interesting in the light of refs.  and  as well as in view of the global renormalization scheme  originally developed by White and co-workers as an extension of Wilson’s phase-space cell method  to the liquid state; it is only by combining tests internal to the theory and comparison with data available by other means that we are able to answer this question.

In this contribution, after a sketchy presentation of the underlying theory itself (section II) and the implementation used (section III), in section IV we turn to the results of applying hrt to square well systems of variable range. After a short summary of the critical point’s location as obtained from simulation-based and other purely theoretical studies of square wells for various values of $`\lambda `$ (sub-section IV A) we first look into the core condition’s implementation, which provides us with a first hint regarding the range of $`\lambda `$ values accessible to hrt in its current formulation and once more highlights the decoupling assumption’s rôle (sub-section IV B). The latter is also implicated in the correct choice of the boundary condition imposed at high density as discussed, alongside the boundary condition’s location’s effect, in sub-section IV C. — A particularly grave aspect of hrt’s numerical side is the stiffness of the partial differential equation (pde) for close-to-critical and sub-critical temperatures (sub-section IV D), the vestiges of which are evident in the results obtained for quasi-continually varying $`\lambda `$ as presented in sub-section IV E. A short summary of our findings and conclusions ends our contribution (section V).

## II The Theory

The definite resource on hrt is the review  by the theory’s original authors, summarizing its formalism as developed in a series of earlier publications  so that we here present only an overview of the theory used and of its formulation, recapitulating some of our earlier findings ; the notation employed here of course co-incides with that of our preceding contribution .

As mentioned already in section I, hrt’s mainstay is the implementation of the suppression of long-wavelength fluctuations characteristic of rg methods by means of a renormalized potential $`v^{(Q)}`$. Thus, rather than directly going from a reference fluid the properties of which are assumed known to the fully interacting system, i. e. from pair potential $`v^{\mathrm{ref}}(r)`$ to $`v(r)=v^{\mathrm{ref}}(r)+w(r)`$, hrt proceeds via a succession of rather artificial  intermediate potentials $`v^{(Q)}(r)`$: For every value of the cut-off wave number $`Q`$, $`v^{(Q)}`$ is given by

$$\begin{array}{c}v^{(Q)}(r)=v^{\mathrm{ref}}(r)+w^{(Q)}(r)\\ \stackrel{~}{w}^{(Q)}(k)=\{\begin{array}{ccc}\stackrel{~}{w}(k)& :& k>Q\\ 0& :& k<Q\end{array}\\ \stackrel{~}{w}^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}(k)=4\pi ϵ\frac{\mathrm{sin}\lambda \sigma k\lambda \sigma k\mathrm{cos}\lambda \sigma k}{k^3},\end{array}$$

where the last line specializes to the square well potential of eq. (1); obviously, $`v^{\mathrm{ref}}`$ and $`v`$ are recovered in the limits $`Q\mathrm{}`$ and $`Q0`$,

$$\begin{array}{c}v^{(\mathrm{})}(r)=v^{\mathrm{ref}}(r)=v^{\mathrm{hs}}(r)\\ v^{(0)}(r)=v(r)=v^{\mathrm{sw}}(r),\end{array}$$

allowing hrt to gradually turn on fluctuations of ever increasing wavelength by lowering $`Q`$ from $`\mathrm{}`$ to zero (numerically , from $`Q_{\mathrm{}}`$ to $`Q_0`$); as mentioned before, criticality (together with non-classical critical exponents) and phase separation (with isotherms rigorously flat in the two phase region) are obtained in the limit $`Q0`$. In this procedure it is essential to maintain the differential picture implied by rg theory and to make sure that the transition from $`Q`$ to infinitesimally smaller cut-off $`Q\mathrm{d}Q`$ is continuous even in the limit $`Q0`$. The latter requirement necessitates replacing the usual free energy $`A^{(Q)}`$ and two-particle direct correlation function $`c_2^{(Q)}(r)`$ of the hypothetical system with cut-off $`Q`$ and potential $`v^{(Q)}(r)`$, the “$`Q`$-system”, by suitably modified quantities, viz.

$$\begin{array}{c}\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)\\ 𝒞^{(Q)}(r)=c_2^{(Q)}(r)+\varphi (r)\varphi ^{(Q)}(r)\\ \varphi =\beta w\beta =1/k_BT,\end{array}$$

where $`\varrho `$ is the number density of the system at hand; the higher order correlation functions $`c_n^{(Q)}(r_1,\mathrm{},r_n)`$, $`n3`$, are free from such problems. (Note that all the direct correlation functions including $`𝒞^{(Q)}(r)`$ are taken to include the ideal gas terms .)

With this set of quantities continuous even in the limit $`Q0`$, viz. $`𝒜^{(Q)}`$, $`𝒞^{(Q)}`$, and the $`c_n^{(Q)}`$, $`n3`$, hrt is derived as a non-terminating hierarchy of coupled odes at fixed density $`\varrho `$, calculating the properties of the $`Q`$-system by treating the system at infinitesimally higher cut-off $`Q+\mathrm{d}Q`$ as a reference system; of these equations, usually only the evolution equation for $`𝒜^{(Q)}`$, viz.

$$\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),$$
(2)
as well as the important compressibility sum-rule

$$\stackrel{~}{𝒞}^{(Q)}(0)=\frac{^2}{\varrho ^2}\left(\frac{\beta 𝒜^{(Q)}}{V}\right)$$
(3)
valid for any cut-off $`Q`$ directly enter practical calculations.

When combined with a closure on the two-particle level, eqs. (2) and (3) define a pde in the $`(Q,\varrho )`$-plane; it is this pde that we will concern ourselves with in the remainder of this text. Said closure, reminiscent of loga/orpa but adding one free parameter to allow imposing thermodynamic consistency as embodied in eq. (3), is given, just as in our earlier contribution , by

$$\begin{array}{cc}\hfill 𝒞^{(Q)}(r,\varrho )& =\varphi (r,\varrho )+\gamma _0^{(Q)}(\varrho )u_0(r,\varrho )+𝒦^{(Q)}(r,\varrho )\hfill \\ \hfill 𝒦^{(Q)}(r,\varrho )& =𝒢^{(Q)}(r,\varrho )+c_2^{\mathrm{ref}}(r,\varrho )\hfill \\ \hfill 𝒢^{(Q)}(r,\varrho )& =_{n=1}^{\mathrm{}}\gamma _n^{(Q)}(\varrho )u_n(r,\varrho ),\hfill \end{array}$$
(4)
where we have generalized to density-dependent potentials. Basis function $`u_0(r,\varrho )`$ is chosen to coincide with  $`w(r,\varrho )/\stackrel{~}{w}(0,\varrho )`$, and the higher basis functions $`u_n(r,\varrho )`$, $`n1`$, vanish outside the core; for our specific choice of basis functions see appendix B of ref. . In order to ensure that both the core condition, i. e. $`g(r,\varrho )=0`$ for $`r<\sigma (\varrho )`$ where $`g(r,\varrho )`$ is the pair distribution function, and sum-rule (3) are met it is necessary to choose the correct set of expansion coefficients $`\gamma _n^{(Q)}(\varrho )`$, $`n0`$, at every cut-off $`Q`$ and for every density $`\varrho `$; assuming their validity for $`Q=\mathrm{}`$ and adopting the short hand notations

$$\begin{array}{c}\alpha ^{(Q)}(\varrho )=\frac{^3}{Q\varrho ^2}\left(\frac{\beta 𝒜^{(Q)}}{V}\right)\\ \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}\end{array}$$

(here, $`\psi `$ is an arbitrary function of $`k`$ and $`\varrho `$), both relations can be combined to

$$\begin{array}{c}\underset{n=1}{\overset{\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}\hfill \\ =\alpha ^{(Q)}(\varrho )\widehat{}^{(Q)}[\stackrel{~}{u}_j(k,\varrho )\stackrel{~}{u}_0(k,\varrho ),\varrho ]\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}$$
(5)
this set of equations must, of course, be truncated to a finite number $`1+N_{\mathrm{cc}}`$ of basis functions, and it is also necessary to neglect non-local contributions to $`\widehat{}^{(Q)}[\psi (k,\varrho ),\varrho ]/Q`$ to allow convenient evaluation at arbitrary $`Q`$. Both of these approximations have been discussed at length in our previous contribution , and while the value of $`N_{\mathrm{cc}}`$ was found to strongly influence the quality of the results obtained, determination of the $`\gamma _n^{(Q)}(\varrho )`$ from eq. (5) and said approximation for the slowly converging $`\widehat{}`$-integrals’ $`Q`$-dependence leads to systematic deficiencies at small $`r`$ in $`g(r)`$ as determined from the Ornstein-Zernike relation. — Unfortunately, for numerical reasons  it is necessary to also adopt the so-called decoupling assumption , viz. $`\alpha ^{(Q)}(\varrho )=0`$; as can easily be seen, this is not only mathematically incompatible with thermodynamic consistency but even suffices to decouple the pde implied by eqs. (2) and (3) to a set of unrelated odes at fixed density only lacking thermodynamic consistency and thus unable to predict clear phase boundaries ; furthermore, we cannot rule out that decoupling may have a significant influence on the solution generated , which is particularly troublesome as the much longer range of $`u_0(r)\varphi (r)`$ when compared with the other basis functions was originally invoked as justification for setting $`\alpha ^{(Q)}(\varrho )=0`$: for square wells, this assumption is certainly even less justified than for the rather long-ranged hard-core Yukawa system ($`z=1.8/\sigma `$) considered in ref. .

Returning to the pde, for the numerical implementation’s benefit we, too, adopted a re-formulation in terms of an auxiliary function $`f(Q,\varrho )`$ simply related to the modified free energy’s derivative with respect to $`Q`$. The details of the procedure leading to a pde of the form

$$\begin{array}{c}\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 ),\\ f(Q,\varrho )\stackrel{~}{u}_0^2(Q,\varrho )=\mathrm{ln}\left(1\frac{\stackrel{~}{\varphi }(Q,\varrho )}{\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )}\right)+\frac{\stackrel{~}{\varphi }(Q,\varrho )}{\stackrel{~}{𝒦}^{(Q)}(Q,\varrho )},\end{array}$$
(6)
and the coefficient functions $`d_{0i}`$ themselves can be found in appendix A of ref. , q. v. ref. .

The above formulation (6) of the problem, of course still coupled to the odes implementing the core condition, obviously has to be amended by initial and boundary conditions; while the former easily follow from $`\gamma _n^{(Q_{\mathrm{}})}=0`$, $`n0`$, (which is sufficient to also determine $`f(Q_{\mathrm{}},\varrho )`$), choice of appropriate boundary conditions is slightly more complicated: if, as is the case in most of the calculations reported here (exceptions see below), the low density boundary is located at $`\varrho _{\mathrm{min}}=0`$, we can make use of the divergence of the ideal gas term $`1/\varrho `$ in $`\stackrel{~}{c}_2^{\mathrm{ref}}`$ to derive not only $`f(Q,0)=0`$ but also $`f(Q,0)/\varrho =0`$ which alone is, in principle, sufficient to uniquely determine the solution up to arbitrarily high density; for computational reasons, however, it is preferable to instead only impose vanishing $`f`$ at $`\varrho _{\mathrm{min}}`$ and to supply an approximate condition for calculating $`f`$ at $`\varrho _{\mathrm{max}}`$. Among the candidates for the constraint to be imposed upon the solution at $`\varrho _{\mathrm{max}}`$ in addition to the core condition there are two we should mention here: Starting with ref. , the so-called orpa-condition, viz. $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})=0`$, has been used extensively; it should, however, be noted that this condition is incompatible with both thermodynamic consistency and the decoupling assumption . An alternative first considered in our previous report  is the decoupling assumption $`\alpha ^{(Q)}(\varrho _{\mathrm{max}})=0`$ itself; of course, this condition is still incompatible with the compressibility sum-rule (3) but this is less of a problem at a boundary where partial derivatives with respect to $`\varrho `$ cannot be evaluated anyway. Another option (not pursued in this contribution) is to give up the core condition altogether, retaining only the lowest basis function $`u_0`$ in the closure (4) and thus effectively replacing eq. (5) by eq. (3); this has the added advantage of mathematical consistency while still retaining the structure of a pde so important for obtaining clear phase boundaries , v. s.

It is one of hrt’s main achievements to allow calculating a fluid’s binodal (coinciding with the spinodal in three dimensions ) without resorting to Maxwell constructions , for subcritical temperatures yielding isotherms rigorously flat in density intervals the boundaries of which are readily identified with the coexisting densities $`\varrho _v`$ and $`\varrho _l`$. Thus, as the isothermal compressibility $`\kappa _T`$ of the fully interacting system, readily found to be proportional to $`\mathrm{exp}(f(\stackrel{~}{\varphi }(0)/\stackrel{~}{𝒦}^{(0)}))1`$ (cf. appendix A of ref. ), diverges in the two-phase region, so must the auxiliary function $`f(Q,\varrho )`$ in the limit $`Q0`$. As a direct consequence of this, the transition from the modified free energy $`𝒜^{(Q)}(\varrho )`$ to $`f(Q,\varrho )`$ is not only computationally convenient but also allows us to follow the isothermal compressibility’s build-up much more easily; even more importantly, a simple analysis  of the behavior of the pde’s coefficients for large $`f(Q,\varrho )`$ readily characterizes the pde as stiff: for any density $`\varrho [\varrho _v,\varrho _l]`$ and close to $`Q=0`$, the true solution $`f(Q,\varrho )`$ oscillates rapidly on a $`Q`$-scale of the order of $`\mathrm{exp}(f)`$, with both an upper bound on the oscillations’ amplitudes and $`f`$’s average slope growing like $`1/Q`$ — needless to say that this behavior cannot be reproduced numerically (v. i. sub-section IV D; q. v. ref. ). Note, however, that it is not an artifact of the re-writing of the pde in the form (6) but rather a problem inherent to hrt itself in a formulation based upon eq. (4) .

## III Numerical procedure

The numerical study of hrt for square well systems of varying range parameter $`\lambda `$ in section IV has only become feasible due to our recent re-implementation of this theory, discussed at length in refs.  and ; we will make use of results obtained with this program exclusively. From a practical point of view, our software provides a means of solving a finite-difference approximation to the pde (6) in an iterated full-approximation scheme, imposing boundary conditions at densities $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$ as well as initial conditions at $`Q=Q_{\mathrm{}}`$, generating a solution for $`Q`$ as low as $`Q_0`$ while ensuring numerical soundness of every step by employing a number of criteria. The pivotal parameter governing all of the numerics is a small quantity denoted $`ϵ_\mathrm{\#}`$ characteristic of the maximum admissible relative error introduced in a single step in the $`Q`$ direction; due to the paramount importance of derivatives with respect to $`\varrho `$, $`ϵ_\mathrm{\#}`$ is strictly related to the coarseness of the density grid.

The only exception to the general strategy of ensuring a numerical quality of $`ϵ_\mathrm{\#}`$ at every step in the calculation is the choice of step sizes $`\mathrm{\Delta }Q`$ in the $`Q`$ direction, at least for sub-critical and close-to-critical temperatures: indeed, in that part of the $`(Q,\varrho )`$-plane where the divergence of the isothermal compressibility builds up, the pde’s stiffness (v. s.) renders fixed-precision arithmetic and relative errors bounded by $`ϵ_\mathrm{\#}`$ incompatible. Consequently, for the calculations reported below we resort to step sizes $`\mathrm{\Delta }Q`$ pre-determined in a way analogous to that employed in earlier applications ; still, monitoring and assessing suitable components of the solution vector in terms of $`ϵ_\mathrm{\#}`$ as described in sub-section III E of ref.  may yield a wealth of information on the numerical process and its evolution.

Most of the calculations reported here have been performed on an equispaced density grid of $`N_\varrho =100`$ density intervals spanning the range from $`\varrho _{\mathrm{min}}=0`$ to $`\varrho _{\mathrm{max}}=1/\sigma ^3`$, corresponding to a value of $`ϵ_\mathrm{\#}=10^2`$; $`N_{\mathrm{cc}}`$ was usually set to 7; and the pre-determined step sizes started from $`\mathrm{\Delta }Q=10^2/\sigma `$ at $`Q_{\mathrm{}}=80/\sigma `$, plunging to a mere $`510^6/\sigma `$ when approaching $`Q_0=10^4/\sigma `$. — When locating the binodal via the divergence of the isothermal compressibility $`\kappa _T`$ we did not require an actual overflow to occur but instead looked for a $`\kappa _T`$-ratio at neighbouring densities exceeding $`10^4`$, which is a rather reliable indicator for the binodal’s location as $`\kappa _T`$ typically jumps by far less than two or by at least some twenty orders of magnitude within one $`\mathrm{\Delta }\varrho `$; the reported values for $`\varrho _v`$ and $`\varrho _l`$ are the mid-points of the density intervals so found. In principle this allows us to locate the coexisting densities and the critical temperature and density to arbitrary precision, even though the computational cost rises sharply with falling $`ϵ_\mathrm{\#}`$.

## IV Application to square wells

As mentioned before, much of the motivation for applying hrt in the formulation outlined in section II to the simple square well model potential is based upon various observations indicating possible limitations of this approximate formulation of hrt for short-ranged potentials. A case in point is the recent work of Caccamo et al.  entirely devoted to several thermodynamically consistent theories’ ability to deal with narrow hard-core Yukawa systems; sure enough, in the case of hrt the shortcomings of the loga/orpa-style closure (4) and, presumably, of the accompanying decoupling assumption underlying the core condition’s implementation via eq. (5) were manifest already in refs.  and  and recently confirmed by us , q. v. ref. .

Of course, any of the problems discussed below only relate to hrt when implemented along the lines of sections II and III and not to hrt proper; however, for reasons discussed in ref. , alternative formulations almost certainly render the numerics far more demanding and open up a whole new suite of problems regarding the numerical implementation’s soundness, especially when performing Fourier transformations of cut-off affected functions .

In the following sub-sections we will complement the discussion of ref.  by further investigation into the numerical nature of hrt; before that, however, it seems pertinent to re-iterate some of the points raised in that publication as far as they concern the reasoning to be put forward in the following. In particular, according to section IV of ref. , for the numerical results to be meaningful the coexisting densities $`\varrho _v`$ and $`\varrho _l`$ must maintain a separation of at least several density grid spacings $`\mathrm{\Delta }\varrho `$ from the boundaries at $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$; consequently, $`\beta `$ should never exceed some maximum value, $`\beta <\beta _{\mathrm{max}}`$, and for the systems considered here and in ref.  and for the typical choices for $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$ the binodal’s proximity to the low density boundary renders $`\beta _{\mathrm{max}}`$ largely density grid- and $`ϵ_\mathrm{\#}`$-independent. — Not to be confused with $`\beta _{\mathrm{max}}`$ is the lowest temperature $`k_B/\beta _{\mathrm{max},\mathrm{\#}}`$ numerically accessible to the program with pre-determined step sizes: this is the temperature below which the program of section III never reaches $`QQ_0`$ or produces abnormal results; note that $`\beta _{\mathrm{max},\mathrm{\#}}`$ may be larger or smaller than $`\beta _{\mathrm{max}}`$, depending on the chosen combination of physical potential, approximations in the formulation used (the boundary conditions in particular), and the choice of parameters affecting the numerical work.

Regarding the implementation of the core condition as sketched in section II, the main conclusion of ref.  was that a minimum of $`N_{\mathrm{cc}}=7`$ basis functions in addition to $`u_0`$ were necessary for acceptable results despite residual defects of $`g(r)`$ close to the origin; a short discussion of the core condition’s slightly different rôle for square wells will be presented below (sub-section IV B).

The critical density $`\varrho _c`$ predicted by hrt, it should be noted, is virtually always in reasonable agreement with literature data as shortly presented in sub-section IV A; indeed, hrt is even able to reproduce the marked rise in $`\varrho _c`$ predicted by refs. , , and  for $`\lambda 1+`$ as opposed to the rigorously constant value in ref. . Due to the satisfactory $`\varrho _c`$-values obtained numerically we will henceforth exclude $`\varrho _c`$ from the discussion; for a demonstration of both $`\varrho _c`$’s insensitivity to variation of parameters of the numerical procedure and the quantitative agreement with the data of sub-section IV A see figures 1 and 2.

In this context it may be of interest that the hrt estimate for the critical density presents no difficulties for the hard-core Yukawa fluid considered in ref. , either, nor is there any mention of such difficulties in any of the other publications on this topic that we are aware of; indeed, the theory’s numerical problems primarily lie in the solution’s small-$`Q`$ behavior for close-to-critical and sub-critical temperatures on the one hand and the use of mutually incompatible assumptions prompted by the need to employ decoupling without giving up thermodynamic consistency on the other hand. Both of these aspects of hrt pertain to different parts of the $`(Q,\varrho )`$-plane, located close to the high density boundary for the rôle mathematical inconsistencies play and at not too large $`Q`$ and $`\varrho \varrho _c`$ for the pathological behavior related to coexistence; they will be discussed in sub-sections IV C and IV D, respectively, and their vestiges will also be seen in the results of applying hrt in the formulation of section II to square wells of quasi-continually varying range $`\lambda ]1,3.6]`$ in sub-section IV E.

### A Non-hrt results for square wells

For comparison purposes we have collected in tables I and II the critical temperatures of various square well systems as obtained from simulations (table I) or by purely theoretical means (table II); the data included have been published within the last decade.

Of the simulation based results included in table I, only those of ref.  for $`\lambda \{1.25,1.375,1.5,1.75,2\}`$ have been obtained by molecular dynamics (MD); most of the other simulation studies rely on one or the other variant of the Monte Carlo (MC) method: Among these, the Gibbs ensemble MC (GEMC) calculations of ref.  set out to determine critical exponents, $`\beta `$ in particular; that work’s finding of $`\beta 1/2`$ for $`\lambda =2`$ as opposed to the expected $`\beta 1/3`$ found for $`\lambda `$ up to 1.75 prompted re-examination of the square well fluid with $`\lambda =2`$ by GEMC augmented by finite-size scaling (FSS) techniques , refuting the mean field value for the effective exponent. — Especially in the critical regime, grand canonical MC (GCMC) simulations incorporating histogram reweighting and FSS offer some advantage over GEMC due to the latter’s restriction to fixed temperature; such an approach has been applied to square wells with $`\lambda =1.5`$ and 3 in ref. ; a more elaborate GCMC scheme not biased towards the Ising universality class and taking into account the YY anomaly has recently been applied to  $`\lambda =1.5`$, v. s. Yet another method goes under the name of thermodynamic- or temperature-and-density-scaling MC (TDSMC); it was applied to the case of $`\lambda =1.5`$ and analyzed in terms of an effective Hamiltonian in refs.  and . — Also included in table I are the results of ref. , employing an MC scheme modified to take advantage of a speed-up possible by combining simulation data with an analytical ansatz for the chemical potential; the efficiency of this approach originally devised to study phase separation allows a large number of systems to be considered. (The error bounds given for these “modified MC” results in table I have been obtained from the different results displayed in ref.  for different parameter settings.)

The theoretical predictions for the critical temperature listed in table II comprise a second-order analytic perturbation theory  (APT2) applicable to $`1<\lambda 2`$ and claimed accurate for $`\lambda 1.4`$ as well as the hard-sphere van der Waals (HSvdW) equation of state . In addition, though not listed in table II, we have utilized the non-square-well-specific Okumura-Yonezawa (OY) estimate  for $`\beta _c`$, primarily as a starting value when looking for the critical temperature in our hrt calculations; for $`v^{\mathrm{sw}[ϵ,\lambda ,\sigma ]}`$, the OY prediction is $`k_BT_c/ϵ=0.203(2\pi /3)\lambda ^30.273`$.

### B The core condition

Ever since application of hrt to continuous fluids started, the implementation of the core condition has been a major issue, probably motivating adoption of the closure (4) and variants thereof for non-hard-sphere reference systems  despite its known deficiencies in the first place; indeed, it is no coincidence that several studies  primarily concerned with the rg aspect of the theory chose to completely eliminate the core condition. When applying hrt as a regular liquid state theory, on the other hand, this is not an option: too great is the effect this may have on both correlation functions and phase behavior . From table III where we compile the critical temperature $`T_c=1/k_B\beta _c`$ for various square well potentials as functions of the number $`N_{\mathrm{cc}}+1`$ of basis functions in the closure (4), just as in ref.  we find virtually constant critical temperatures for $`1N_{\mathrm{cc}}4`$; on the other hand, the amount of variation seen upon further increasing $`N_{\mathrm{cc}}`$ strongly depends on $`\lambda `$, which immediately carries over to the pair distribution function $`g^{(Q_0)}(r,\varrho )`$ and its compatibility with the core condition: For $`\lambda =3`$, the longest ranged potential considered in table III, $`g^{(Q_0)}(r,\varrho )=0`$, $`r<\sigma `$, holds reasonably well except very close to $`r=0`$ even for $`N_{\mathrm{cc}}=1`$; when increasing the number of basis functions all the way to $`N_{\mathrm{cc}}=10`$, the pair distribution function has to be corrected for very small $`r`$ only, yielding a $`|g^{(Q_0)}(r,\varrho )|`$ that remains bounded by some $`10^2`$ of the contact value $`g^{(Q_0)}(\sigma +,\varrho )`$ for all $`r<\sigma `$; the corresponding small change in $`g^{(Q_0)}(r,\varrho )`$ and $`𝒞^{(Q_0)}(r,\varrho )`$ is reflected in the near-constant predictions for $`\beta _c`$ evident from table III. Similarly, for $`\lambda =1.5`$ and $`\lambda =2`$ and within the $`N_{\mathrm{cc}}`$-range considered, the implementation of the core condition does not convincingly improve except for supercritical temperatures and intermediate densities; this time, however, the pair distribution functions remain far from compatible with the core condition even for $`N_{\mathrm{cc}}=10`$, and neither $`\beta _c`$ nor $`g^{(Q_0)}(r,\varrho )`$ itself nor, for that matter, the final values of the loga/orpa expansion coefficients $`\gamma _n^{(Q_0)}(\varrho )`$ indicate that the expansion (4) for $`\stackrel{~}{𝒞}^{(Q)}(k,\varrho )`$ might be close to convergence. But if the quality of $`g^{(Q_0)}(r,\varrho )`$ improves only little if at all, the remaining deficiencies are probably to be blamed on the approximation for the poorly convergent integrals’ derivative with respect to $`Q`$ mentioned earlier (cf. eq. (12) of ref. ) rather than on an insufficient number of basis functions; on the other hand, even though the decoupling assumption cannot directly affect the pair distribution function’s compliance with the core condition, the approximation of neglecting the non-local term in $`\widehat{}^{(Q)}[\psi (k,\varrho ),\varrho ]/Q`$ is on the same level as that of setting $`\alpha ^{(Q)}(\varrho )=0`$, as was stressed by the authors of ref.  upon jointly introducing these two assumptions. Thus, combining the above findings regarding the core condition with the analogous analysis of sub-section IV of ref.  and with that contribution’s investigation into the decoupling assumption’s possible effect (cf. fig. 2 of ref. ) we are led to the conclusion that decoupling poses certainly no less a problem here than for the hard-core Yukawa potential studied there.

### C High-density boundary condition

Numerically, there are two ways for the implementation of section III to fail to reach $`Q=Q_0`$, both, of course, easily detected by the “monitoring” variant of our code (cf. sub-section III E of ref. ): due to the solution’s pathological behavior wherever $`f(Q,\varrho )`$ is large (cf. sub-section IV D), or because of inappropriate boundary conditions at high density. As for the latter — an issue intimately linked to the decoupling assumption —, the immediate reason for the program’s failure is a near-discontinuity in the numerical solution close to the boundary: For the moment setting aside the decoupling assumption and other approximations, in the application of hrt with the closure (4) at any point $`(Q,\varrho )`$ in the interior of the pde’s domain the core condition uniquely determines the $`\gamma _n^{(Q)}(\varrho )`$, $`n1`$, for given $`\gamma _0^{(Q)}(\varrho )`$; this expansion coefficient is then determined by imposing thermodynamic consistency as embodied in the compressibility sum-rule (3). At a boundary, i. e. for $`\varrho \{\varrho _{\mathrm{min}},\varrho _{\mathrm{max}}\}`$, however, the second density derivative cannot be evaluated reliably so that some other condition must be imposed; in the calculations reported here (with the obvious exception of those for fig. 4) we always choose $`\varrho _{\mathrm{min}}=0`$ so that the divergence of the ideal gas term in $`\stackrel{~}{c}_2^{(Q)}`$ provides $`f(Q,\varrho )=0`$ as a convenient and unproblematic boundary condition. For $`\varrho =\varrho _{\mathrm{max}}`$, on the other hand, we are in principle free to use any suitable approximation for the structural and thermodynamic properties of the $`Q`$-system and to calculate $`f(Q,\varrho _{\mathrm{max}})`$ from said approximation, thereby providing the necessary boundary condition for the pde (6); but for practical reasons it is desirable to use the same loga/orpa-form for the $`Q`$-system’s direct correlation function at $`\varrho _{\mathrm{max}}`$ as in the rest of the problem’s domain so that, in particular, the loga/orpa prescription $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})=0`$ is a natural choice of boundary condition. In general, however, due to the pde’s diffusion-like character any condition imposed at $`\varrho _{\mathrm{max}}`$ that is incompatible with the solution for $`\varrho <\varrho _{\mathrm{max}}`$ by necessity induces a corresponding near-discontinuity in $`f(Q,\varrho )`$ close to the boundary; within the framework of a finite difference scheme this is reflected in a mismatch of $`f(Q,\varrho _{\mathrm{max}})`$ and the solution at densities close by, i. e. $`f(Q,\varrho _{\mathrm{max}}i\mathrm{\Delta }\varrho )`$ for small $`i1`$, and the mismatch’s severity may serve as a direct measure for the inappropriateness of the boundary condition at $`\varrho _{\mathrm{max}}`$ in relation to the approximations applied at densities in $`]\varrho _{\mathrm{min}},\varrho _{\mathrm{max}}[`$.

On the other hand, the numerics become intractable unless we adopt the decoupling assumption, and the only way to consistently use $`\alpha ^{(Q)}(\varrho )=0`$ without abandoning the core condition is to decouple the hrt-pde to a set of odes at fixed density only ; this, unfortunately, removes all traces of thermodynamic consistency from the equations and thereby precludes obtaining clear phase boundaries . It is therefore necessary to restrict decoupling to the implementation of the core condition only while retaining the structure of a pde together with the compressibility sum rule (3) despite the latter’s incompatibility with decoupling. Thus, for $`\varrho _{\mathrm{min}}<\varrho <\varrho _{\mathrm{max}}`$, both $`\stackrel{~}{𝒞}^{(Q)}(0,\varrho )=^2\left(\beta 𝒜^{(Q)}/V\right)/\varrho ^2`$ and $`\alpha ^{(Q)}(\varrho )=0`$ are used for different parts of the problem; at $`\varrho _{\mathrm{max}}`$, however, again any approximation allowing calculation of $`f(Q,\varrho _{\mathrm{max}})`$ may be used, so that it is tempting to once again resort to the loga/orpa-condition of vanishing $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})`$ or variants thereof.

But due to the decoupling assumption’s possibly large effect, any boundary condition that does not incorporate $`\alpha ^{(Q)}(\varrho _{\mathrm{max}})=0`$ — and bear in mind that $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})`$ and $`\alpha ^{(Q)}(\varrho _{\mathrm{max}})`$ cannot both vanish at the same time for generic cut-off $`Q`$ — will once again incur a fatally large mismatch; if, however, we must resort to decoupling anyway, it seems preferable to consistently apply it for the boundary condition rather than to inconsistently combine it with a condition alien to the theory; also, though the mismatches’ magnitudes from imposing $`\alpha ^{(Q)}(\varrho _{\mathrm{max}})=0`$ alone or from mixing it with the loga/orpa condition $`\gamma _0^{(Q)}(\varrho _{\mathrm{max}})=0`$ generally do not differ much as long as the pde’s stiffness does not play a rôle (e. g., for $`\lambda =1.5`$, as long as we restrict ourselves to $`Q8/\sigma `$ or higher, or to $`\beta \beta _c`$), the former approach fares better than the other one more often than not. It is only in this sense, i. e. presupposing a loga/orpa-like ansatz even at $`\varrho _{\mathrm{max}}`$ and application of decoupling in the implementation of the core condition according to eq. (5) at all $`\varrho `$, that the results are largely independent of the choice of boundary condition as claimed, for $`\beta <\beta _c`$, in ref. .

In the numerical work we find that such a mismatch is present whenever the calculation proceeds via mathematically inconsistent or conflicting approximations; in the case of square wells with their comparatively short potential range, however, the problems are much more severe than in other systems so that $`\beta _{\mathrm{max},\mathrm{\#}}`$ is rather small and even drops below $`\beta _c`$ for most of the $`\lambda `$ interval from 1 to 2 (cf. sub-section IV E). Restricting ourselves to $`\beta <\beta _{\mathrm{max},\mathrm{\#}}`$ and $`Q=Q_0`$, the mismatch is typically reflected in an increase by one order of magnitude in the three-point finite-difference estimate of, e. g., $`|^2f(Q_0,\varrho )/\varrho ^2|`$ right at the boundary over the near-constant values at slightly lower densities; apart from a positive correlation with $`ϵ_\mathrm{\#}`$, the mismatch’s severity is qualitatively unaffected by a change in parameters of the numerical procedure or the choice and location of the boundary condition (with the above provisions).

Another effect worth mentioning in connection with the boundaries is the influence their locations, viz. $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$, may have. The basic mechanism and its implications for the coexisting densities were already mentioned in the opening remarks of this section; here we only want to point out that the non-criticality enforced by the boundary conditions not only may unduely distort the binodal predicted by hrt as demonstrated in fig. 1, very small $`\varrho _{\mathrm{max}}`$ may also allow one to reach $`Q=Q_0`$ at higher $`\beta `$, thus effectively raising $`\beta _{\mathrm{max},\mathrm{\#}}`$ while lowering $`\beta _{\mathrm{max}}`$. — Sometimes, however, the expectation of the binodal keeping a separation from the boundary of several $`\mathrm{\Delta }\varrho `$ at least does not hold, and a preposterous two-phase region appears very close to $`\varrho _{\mathrm{max}}`$ or, very rarely, close to $`\varrho _{\mathrm{min}}`$; e. g., for $`\lambda =1.88`$ and $`\beta =0.392/ϵ`$ the equations can be solved all the way down to $`Q=Q_0=10^4/\sigma `$, predicting an unrealistic two-phase region extending from $`0.845(5)/\sigma ^3`$ to $`0.995(5)/\sigma ^3`$. This behavior turns out to come from the interplay of the mismatch at $`\varrho _{\mathrm{max}}`$ and the numerical treatment of the stiffness of the pde.

### D The region of large $`f(Q,\varrho )`$

For subcritical temperatures the hrt-pde’s true solution’s erratic behavior in that part of the $`(Q,\varrho )`$-plane where $`f(Q,\varrho )`$ is large and the isothermal compressibility’s divergence is built up (cf. section II, q. v. ref. ) obviously eludes reliable numerical realization; in particular, while $`ϵ_\mathrm{\#}`$ still characterizes the level of accuracy in auxiliary calculations, the same can no longer be true for the accuracy of the pde’s discretization as this would require step sizes $`\mathrm{\Delta }Q`$ so small as to cause floating point underflow upon evaluating, e. g., $`Q(Q\mathrm{\Delta }Q)`$, thus rendering finite differences numerically insignificant.

Consequently, in this respect we have to give up our strategy of controlling the numerical procedure so as to locally ensure a quality of $`ϵ_\mathrm{\#}`$ at least, turning to pre-determined step sizes  $`\mathrm{\Delta }Q`$ in addition to fixed $`\mathrm{\Delta }\varrho `$, to which similar concerns apply ; on such a coarse mesh of $`(Q,\varrho )`$-points underlying the finite difference scheme, however, the true solution cannot even be represented adequately, and the numerical approximation for $`f(Q,\varrho )`$ obtained from the pde’s discretization with these far too large step sizes cannot be trusted to faithfully represent even the average behavior of $`f(Q,\varrho )`$.

This inadequacy of the step sizes is reflected in various peculiarities of the solution vector obtained in the numerical procedure; indeed, when monitoring the evolution of $`f(Q,\varrho )`$ and the core condition coefficients $`\gamma _n^{(Q)}(\varrho )`$, our code readily detects the plummeting step sizes necessary and signals the incompatibility of the behavior seen with the assumption of smoothness underlying finite difference schemes. Another telltale sign is iterated corrector steps’ failure to converge when $`f(Q,\varrho )`$ is large: even though implicit schemes like the one we employ  are the standard treatment for stiff systems, the rapid growth of the oscillations’ amplitudes renders the finite difference equations themselves unstable under iteration; only when resigning on any control of the numerical error and refraining from iterations of the corrector step do the step sizes $`\mathrm{\Delta }Q`$ chosen allow one to force advancing $`Q`$ all the way to $`Q_0`$ in remarkably many cases. Also, comparison of $`f(Q_0,\varrho )`$ as obtained with different sets of step sizes $`\mathrm{\Delta }Q`$ reveals that, for $`\varrho _v<\varrho <\varrho _l`$, the evolution of $`f(Q,\varrho )`$ seen numerically is driven by the number and size of $`Q`$-steps only and certainly does not correspond to an average over oscillations ; the same mechanism is also responsible for a small $`\mathrm{\Delta }Q`$-dependence of the critical temperature $`\beta _c`$. — By the same token, due to the $`d_{01}`$\- and $`d_{02}`$-terms in eq. (6), the pde’s stiffness and the related problems have a direct bearing on the solution outside the coexistence region even if the numerical predictions there turn out rather insensitive to variation of parameters of the numerical procedure; in particular, we expect a gradual but non-negligible distortion (in addition to the effects of numerical differentiation close to the near discontinuity) of the binodal, increasing with falling temperature.

### E hrt results for square wells

In the light of the preceding exposition as well as of the discussion in ref.  it may at first seem surprising that hrt in the formulation of section II has a record of being well applicable to a variety of systems (cf. section I); also, as we shall see in a moment, even for square wells, a system expected to be particularly vulnerable to the problems just outlined, we find reasonable estimates of the critical points’ locations for a wide range of $`\lambda `$-values. Still, the mechanisms sketched in subsections IV C and IV D as well as the difficulties presented in ref.  remain and manifest themselves numerically in a number ways.

For a first orientation, let us look at the results summarized in figs. 2 and 3 where the critical temperature $`T_c`$ and density $`\varrho _c`$ are shown as functions of $`\lambda `$; the underlying calculations have been obtained with $`ϵ_\mathrm{\#}=10^2`$, imposing decoupling in a consistent way at $`\varrho _{\mathrm{max}}=1/\sigma ^3`$ and with $`N_{\mathrm{cc}}+1=7+1`$ basis functions in the expansion (4) of the loga/orpa-function $`𝒢^{(Q)}`$. With the exception of some spurious results at $`\lambda 1.1`$, whereever $`\beta _c<\beta _{\mathrm{max},\mathrm{\#}}`$ the critical temperature in general compares quite favorably with the data of tables I and II; from the calculations we have performed for a large number of systems in the range $`1<\lambda 3.6`$ and ignoring some isolated results, a critical point is found for $`1.06\lambda 1.24`$, for $`1.45\lambda 1.53`$, and for $`\lambda 1.939`$; calculations with $`N_{\mathrm{cc}}=5`$ yield analogous results , with $`\beta _c<\beta _{\mathrm{max},\mathrm{\#}}`$ in a somewhat larger part of the parameter range, viz. for $`1.09\lambda 1.58`$ and for $`\lambda 1.896`$, but will not be considered in the following in view of the considerations of sub-section IV B and of other defects that turn out to be larger than for $`N_{\mathrm{cc}}=7`$.

For the moment setting aside the data for $`\lambda <1.939`$, hrt’s predictions for the critical temperature are generally found to be in satisfactory agreement with the $`\beta _c(\lambda )`$-curve expected from the simulation-based and theoretical results presented in sub-section IV A. Embedded into this regular overall behavior of $`\beta _c`$ as a function of $`\lambda `$, however, we find a number of depressions and elevations of $`\beta _c`$, some of which cannot be seen on the scale of the plot 2 but from the numeric results only ; others, however, are so strong as to render the critical temperature a non-monotonic function of $`\lambda `$, which is certainly not expected on the grounds of the literature presented in sub-section IV A, the data of refs.  in particular.

In the light of sub-section IV D it is of course tempting to simply attribute this behavior to the difficulties previously discussed, especially since the critical point is located in the region of large $`f(Q_0,\varrho )`$ by definition; the peculiar distribution of $`\lambda `$-values affected, however, suggests that these problems of the numerical procedure are triggered by a special mechanism. Indeed, a closer look at the core condition function $`\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )`$ for fixed density $`\varrho `$ reveals, for every single one of the $`\lambda `$ values implicated that we checked, that the combination of terms pertaining to $`w(r)`$ or $`v^{\mathrm{hs}}(r)`$ alone (of ranges $`\lambda \sigma `$ and $`\sigma `$, respectively) regularly and quite frequently reduces the amplitude of this function’s swings about the ideal gas value of $`1/\varrho `$; the same happens only occasionally for $`\lambda `$-values removed from these irregularities so that it is, in fact, possible to quite reliably determine whether or not a given $`\lambda `$ is affected from a plot of $`\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )`$ for $`\varrho \varrho _c`$ alone as illustrated in fig. 4. It will come as no surprise that most of the irregularities occur when $`\lambda `$, the ratio of the two characteristic lengths present in the model, is close to a simple fraction: among the shifts in $`T_c`$ most obvious are those where $`\lambda `$ is close to $`2`$ (cf. fig. 3), $`2\frac{1}{4}`$, $`2\frac{1}{7}`$, $`2\frac{1}{9}`$ and $`2\frac{1}{12}`$, and in retrospect it seems justified to also include the small parameter range around $`\lambda =1\frac{1}{2}`$ in this list, v. i.; the effect is less obvious from fig. 2 but still discernible at $`2\frac{1}{2}`$, $`2\frac{1}{3}`$ and $`2\frac{2}{3}`$, whereas for $`2\frac{1}{4}`$ and $`2\frac{3}{4}`$ it is so small as to make the plot of $`\beta _c(\lambda )`$ appear smooth while the irregularities are still evident from the numerical values; also note that, once again, $`\varrho _c`$ is hardly affected.

All these observations seem to indicate that indeed it is the interplay of the two different length scales and the resulting partial oppression of a significant portion of the oscillations of $`\stackrel{~}{𝒞}^{(Q)}(Q,\varrho )`$ that cause the discrepancy of hrt and literature results for the critical temperature around certain $`\lambda `$ values even though a smooth interpolation of hrt’s predictions from $`\lambda `$ values nearby is well compatible with the data presented in sub-section IV A. Even though we currently cannot pinpoint the precise mechanism by which this unphysical behavior of hrt arises and, in particular, cannot distinguish between the closure’s inadequacy and the pde’s stiffness as the main culprit — though the latter is certainly implicated to some degree —, two conclusions may be drawn quite safely: for one, as long as we stay clear of values of $`\lambda \text{ }>2`$ that are close to simple fractions, or restrict ourselves to $`\lambda \text{ }>2.7`$ where the effects are rather small, we can probably trust the numerical results — with the caveats of ref.  and sub-sections IV C and IV D — to the same degree of confidence as those obtained for the hard-core Yukawa system in ref. . And secondly, it is only in the presence of discontinuities in the potential that certain lengths feature prominently in the relevant functions’ Fourier transforms and can so give rise to problems of the kind outlined above; consequently, as long as we confine ourselves to continuous $`w(r,\varrho )`$, which still includes most of the potentials popular in liquid state physics, the unphysical shifts in $`\beta _c`$ seen for certain parameter combinations are likely not an issue, whereas the same problems are expected to resurface, e. g., for the multi-step potential also defined in sub-section III B of ref. .

Another lesson to be drawn from the findings presented here as well as in ref.  is that, as a general rule, conclusions should never be drawn from isolated results alone; it is only through the combination and meticulous scrutiny of a set of related calculations that meaningful information can be extracted from hrt calculations: due to the problems related to the implementation of the core condition, to the nature and location of the boundary conditions, and to the pde’s stiffness, any single calculation must be considered as of uncertain standing. As an example , the analogue of fig. 1 for $`\lambda =1.5`$ shows a considerably larger variation in the binodal and the critical point’s location, which is consistent with the above conclusions regarding the reason for $`\beta _{\mathrm{max},\mathrm{\#}}`$ rising above $`\beta _c`$ in a narrow region around this $`\lambda `$ value, whereas any one of the phase boundaries found in itself is a perfectly plausible candidate for the “true” hrt binodal.

## V Conclusion

In conjunction with the findings of ref. , the discussion of section IV provides quite coherent a picture of hrt’s numerical side as well as of some peculiarities encountered for square wells. Most prominently, we see a marked dependence of the quality of the results on the potential’s range, confirming the trend of decreasing accuracy for narrower potentials reported  for the hard-core Yukawa fluid; it has long been accepted  that the simplistic but computationally convenient  closure (4) has a part in this, and an improved closure has recently been proposed . Still, as far as numerical application of hrt is concerned the closure cannot be discussed without reference to the decoupling assumption and to the approximate implementation of the core condition via odes coupled to the hrt-pde; and while the former has been found problematic both for square wells (present contribution) and for the hard-core Yukawa fluid considered in ref.  and should probably not be trusted easily for any system, the severity of the difficulties brought about by the simplified treatment of the core condition sensitively depends on the potential type and parameters chosen: for the continuous and rather long ranged Yukawa potential with $`z=1.8/\sigma `$, $`g^{(Q_0)}(r,\varrho )`$ can be made sufficiently small within the core, and the square well fluid with $`\lambda =3`$ fares equally well at least; from the discussion of sub-section IV B and the data of table III, however, it becomes apparent that smaller $`\lambda `$ — we have looked at $`\lambda =1.5`$ and $`\lambda =2`$ in particular — incurs substantial problems, with residual defects in the pair distribution functions attributed to the ill-justified approximation of neglecting some slowly converging integrals  in the core condition.

But table III demonstrates not only the $`\lambda `$-dependence of the results’ sensitivity to the number $`N_{\mathrm{cc}}+1`$ of basis functions retained in the truncated eq. (5) when varying $`N_{\mathrm{cc}}`$ in the range $`0N_{\mathrm{cc}}10`$: while the virtually constant critical temperature predicted for $`\lambda =3`$ seems trustworthy and is, indeed, well compatible with simulation results (cf. table I), the amount of variation in $`\beta _c`$ for $`\lambda =1.5`$ and, to a much lesser degree, for $`\lambda =2`$ precludes accurate determination of the critical temperature; this is a first indication that the theory might be able to handle square wells with $`\lambda =3`$ quite reliably whereas problems cannot be denied for $`\lambda =2`$, and $`\lambda =1.5`$ seems largely out of reach for hrt in the present formulation. This is confirmed by the results obtained by quasi-continuous variation of $`\lambda `$ in the range $`1<\lambda 3.6`$ as shown, for $`N_{\mathrm{cc}}=7`$, in fig. 2: the critical point is accessible only in part of this parameter range, and not only the critical temperature at fixed $`\lambda `$ but also the boundaries of the $`\lambda `$-intervals where hrt is able to reach temperatures as low as $`T_c`$ strongly depend on $`N_{\mathrm{cc}}`$ (cf. sub-section IV E).

Applicability of hrt to only a restricted $`\lambda `$-range is, of course, again related to the pronounced short-rangedness of the square well potential; to explain it, however, we have to invoke not only the loga/orpa-style closure and the approximate implementation of the core condition but also the other difficulties encountered in the numerical procedure as highlighted in this and our preceding contribution, viz. the decoupling assumption (ref. ), inappropriate boundary conditions (sub-section IV C) and, most importantly, the pde’s stiffness for thermodynamic states of high compressibility (sub-section IV D). All of these are, in principle, always present to some degree in numerical applications of hrt; it may prove valuable that sub-sections IV C and IV D provide distinct signatures readily detected by the implementation of section III that might be looked out for in more advanced applications of the theory, too. — Related to these difficulties is a peculiar effect specific to square wells (sub-section IV E): close to certain $`\lambda `$ values, simple fractions in particular, we see shifts in the critical temperature that render hrt’s predictions much less compatible with simulations and other theoretical descriptions of the square well fluid than would be expected from the results obtained for $`\lambda `$ values close by; the mechanism for triggering these local distortions of the $`\beta _c(\lambda )`$-curve is illuminated at least to the point of linking it to the presence of a discontinuity in the potential’s perturbational part. All in all, the numerical evidence as well as comparison with literature data suggest that the formulation of hrt sketched in section II is well able to deal with square wells and to locate their critical points to reasonable accuracy for $`\lambda \text{ }>2`$ as long as certain values are avoided, or else for $`\lambda \text{ }>2.7`$.

## VI Acknowledgments

The authors gratefully acknowledge support from Österreichischer Forschungsfonds under project number P13062-TPH.


format: text

text:
# Ultralocal solutions for quantum integrable nonultralocal models

## Abstract

A challenge in the theory of integrable systems is to show for every nonultralocal quantum integrable model, a possible connection to an ultralocal model. Some of such gauge connections were discovered earlier. We complete the task by identifying the same for the remaining ones along with two new models. We also unveil the underlying algebraic structure for these nonultralocal models.

PACS numbers: 02.30.Ik, 03.65.Fd, 02.20.Uw, 11.10.Lm

Key Words: Nonultralocal quantum integrable models, Braided Yang-Baxter equation, Discrete current algebra, Gauge connection with ultralocal models

1. Introduction

Impressive success has been achieved in the theory of quantum integrable systems during the last twenty years. However the major progress made and the important results obtained in the subject are limited mostly to a class of models, known as the ultralocal (UL) class. The Liouville model, sine-Gordon (SG) model in laboratory coordinates, relativistic Toda chain, nonlinear Schrödinger equation (NLS) etc. are few well known examples of this class. The UL models must satisfy the ultralocality condition, i.e. their representative Lax operators at different lattice sites must commute: $`L_{2k}^{ul}(\mu )L_{1j}^{ul}(\lambda )=L_{1j}^{ul}(\lambda )L_{2k}^{ul}(\mu )\text{for }jk,`$ which is essential for extending the Yang-Baxter equation (YBE): $`R_{12}(\lambda \mu )L_{1j}^{ul}(\lambda )L_{2j}^{ul}(\mu )=L_{2j}^{ul}(\mu )L_{1j}^{ul}(\lambda )R_{12}(\lambda \mu )`$ to its global form: $`R_{12}(\lambda \mu )T_1^{ul}(\lambda )T_2^{ul}(\mu )=T_2^{ul}(\mu )T_1^{ul}(\lambda )R_{12}(\lambda \mu )`$, for the monodromy matrix $`T_a^{ul}(\lambda )=_{j=1}^NL_{aj}^{ul}(\lambda )`$. Note that $`a=1,2`$ labels here the auxiliary or matrix spaces, while $`j=1,2,\mathrm{},N`$ denotes the associated quantum spaces. For the models with periodic boundary condition the global YBE leads in turn to the required trace identity $`[\tau (\lambda ),\tau (\mu )]=0`$, for the transfer matrix $`\tau (\lambda )=tr_a(T_a^{ul}(\lambda ))`$, establishing the quantum integrability of the system by showing mutual commutativity: $`[c_n,c_m]=0`$, for the set of conserved operators obtained from $`\tau (\lambda )=_nc_n\lambda ^n`$ or from any other function of it. The eigenvalue problem for all these conserved operators including the Hamiltonian of such models can be solved exactly by the powerful algebraic Bethe ansatz, which effectively exploits the global form of YBE . The ultralocality condition reflects generally the fact that the basic fields involved in the Lax operator of an UL model are of canonical nature, which commute trivially at space-like separated points $`xy.`$

However, there exists on the other hand another rich class of integrable models which violate the ultralocality condition and make their quantum description through standard YBE formulation difficult. In such nonultralocal (NUL) models the basic fields are either of noncanonical type having nontrivial commutators at $`xy`$ or their associated Lax operators contain space-derivatives of canonical fields. Nevertheless there is now a considerable number of NUL systems, namely, nonabelian Toda chain , the quantum mapping , model related to the Coulomb gas picture of CFT , current in WZWN , integrable model on moduli space  and the quantum mKdV , for which the quantum integrability is established and the braided extensions of YBE have been formulated . We shall call such NUL models genuine quantum integrable models. However, there are some NUL models like the SG in light-cone coordinates and the nonlinear $`\sigma `$-model, for which the quantum integrability through braided YBE still remains unsolved.

Nevertheless, for the complete understanding of NUL models and for fully exploiting the powerful machinery of quantum integrable systems developed already for UL models, it is highly desirable to find a possible connection for these NUL systems to UL models. Though in general such a connection is not guaranteed, an UL solution was first discovered for the nonabelian Toda chain using a operator dependent local gauge transformation . Subsequently, similar relation was unveiled also for other quantum NUL systems like current algebra in WZWN model , Chern-Simon theory related integrable model on moduli space  and the mKdV model .

It is therefore tempting to conjecture that such a connection to UL model should exist for every NUL system, at least for the genuine quantum integrable models that respect the braided YBE. However the problem is that, the well known NUL models like the Coulomb gas picture of CFT (CGCFT)  and the quantum mapping , though fit well into the braided YBE formulation  are not known to have any established UL connection. Moreover for the quantum light-cone SG (LCSG) model, not only its explicit UL relation is unknown, but also, as mentioned already, its quantum integrability through braided YBE could not be formulated. Our aim here is to find the required NUL-UL connection for all these models and thus to resolve the eminent problem for the validity of the above conjecture, at least for the genuine quantum integrable models discovered till today. Along with the above result we find also a new NUL sine-Gordon (NSG) like system and able to identify its connection with a Liouville-like UL model on a lattice. A preliminary announcement on two new models investigated here, e.g. LCSG and NSG, has been made in .
2. Nonultralocal models and their ultralocal connections

The NUL models of our concern satisfy the braiding relation

$$L_{2j+1}^{nul}(\mu )Z_{12}^1L_{1j}^{nul}(\lambda )=L_{1j}^{nul}(\lambda )L_{2j+1}^{nul}(\mu )$$
(1)
involving nearest-neighboring (NN) sites, together with the braided YBE (BYBE)

$$R_{12}(\lambda \mu )L_{1j}^{nul}(\lambda )L_{2j}^{nul}(\mu )Z_{12}^1=L_{2j}^{nul}(\mu )L_{1j}^{nul}(\lambda )Z_{21}^1R_{12}(\lambda \mu ).$$
(2)
Note that along with the $`R`$-matrix as in UL models an additional braiding matrix $`Z`$ enters into these equations. Remarkably the nonultralocality (1) can be incorporated in the BYBE to obtain its global form for the corresponding monodromy matrix $`T_a^{nul}(\lambda )=L_{aN}^{nul}(\lambda )\mathrm{}L_{a1}^{nul}(\lambda )`$ with $`a=1,2`$ as

$$R_{12}(\lambda \mu )T_1^{nul}(\lambda )Z_{21}^1T_2^{nul}(\mu )Z_{12}^1=T_2^{nul}(\mu )Z_{12}^1T_1^{nul}(\lambda )Z_{21}^1R_{12}(\lambda \mu ),$$
(3)
which leads again to the trace identity establishing the quantum integrability for the NUL models. It may be noticed that unlike the YBE the global form of BYBE differs a bit from its local form. This happens because for the periodic models with $`N+11`$, the Lax operators $`L_{aN}^{nul}(\lambda )`$ and $`L_{b1}^{nul}(\lambda )`$, $`a,b=1,2`$ become in fact NN operators exhibiting nontrivial commutation relations due to (1). It can be seen easily that by putting $`Z=1`$ in the above NUL equations, i.e. braiding relation and the BYBEs, one obtains the corresponding equations, i.e. the ultralocality condition and the standard YBEs for the UL models.

The BYBE (2) provides similar Hopf algebra structure to the integrable NUL models, though braiding relation (1) induces now braided algebra property by modifying the multiplication rule . At the same time the trace identity for the NUL models as obtained from (3) is also modified to $`tr_{12}(T_1^{nul}(\lambda )Z_{21}^1T_2^{nul}(\mu )Z_{12}^1)=tr_{12}(T_2^{nul}(\mu )Z_{12}^1T_1^{nul}(\lambda )Z_{21}^1)`$ and consequently its crucial factorizability depends heavily on the the structure of the braiding matrix $`Z`$, adding another difficulty in tackling such models . Therefore though some progress has been made , as mentioned above, in direct analysis of integrable NUL models, identifying explicit gauge relation of such systems to UL models would be significantly important for getting insight into their underlying algebraic structures as well as their exact solutions. Following the idea of  we aim to find such UL source models for genuine quantum integrable NUL models through a gauge transformation like $`L_j^{nul}=D_{j+1}L_j^{ul}D_j^1`$, such that $`L_j^{ul}`$ satisfies the YBE along with the ultralocality. This in turn is guaranteed, as can be shown by some simple algebraic manipulation, by the condition

$$D_{aj}L_{bj}^{nul}D_{aj}^1=L_{bj}^{nul}Z_{ab}^1,a,b=1,2,ab,$$
(4)
provided the gauge operator $`D_j`$ satisfies the YBE together with a related equation as

$$R_{12}(\lambda \mu )D_{1j}D_{2j}=D_{2j}D_{1j}R_{12}(\lambda \mu ),R_{12}(\lambda \mu )D_{2j}D_{1j}=D_{1j}D_{2j}R_{12}(\lambda \mu ).$$
(5)
However in the particular case when $`[D_{1j},D_{2j}]=0,`$ both the equations (5) reduce to a single condition $`[R_{12},D_{1j}D_{2j}]=0`$, which actually holds true for gauge operators in all the examples we have considered here. Moreover when $`D_{bj}`$ and $`Z_{ab}`$ mutually commute (which is true for all our examples except that with the quantum mapping), condition (4) holds also for the gauge related $`L_j^{ul}`$. Interestingly the monodromy matrices for these periodic models are linked now as

$$T^{nul}=\underset{j}{\overset{N}{}}L_j^{nul}=D_{N+1}(\underset{j}{\overset{N}{}}L_j^{ul})D_1^1=D_1T^{ul}D_1^1,$$
(6)
and satisfy the relation

$$D_{a1}T_b^{nul}D_{a1}^1=T_b^{nul}Z_{ab}^1,a,b=1,2,ab,$$
(7)
as can be shown starting from (4) with the assumption $`[D_{aj+1},L_{bj}^{nul}]=0,ab,`$, which is fulfilled in our cases. We should mention here that though the global BYBE (3) can be derived from its local form and the braiding relation independent of its conjectured relation with the UL model, by using such relations (6) and (7) one can actually derive (3) starting from the standard YBE for the corresponding $`T^{ul}`$. Inspired by the fact that the UL model must have the same $`R`$-matrix as its NUL counterpart related by a gauge transformation, we intend to seek our $`L^{ul}`$ solutions using the generating scheme of  based on the $`R`$-matrix classification.

Recall that the well known trigonometric $`4\times 4`$ $`R(\lambda )`$-matrix is defined by its nontrivial elements

$$R_{11}^{11}=R_{22}^{22}=\mathrm{sin}(\lambda +\alpha ),R_{12}^{12}=R_{21}^{21}=\mathrm{sin}\lambda ,R_{21}^{12}=R_{12}^{21}=\mathrm{sin}\alpha ,$$
(8)
while the ancestor Lax operator of quantum UL models associated with (8) and satisfying the YBE may be given  by

$$L_k^{tanc(ul)}(\lambda )=\left(\begin{array}{c}c_1^+\xi q^{S_k^3}+c_1^{}\frac{1}{\xi }q^{S_k^3}S_k^{}\\ S_k^+c_2^+\xi q^{S_k^3}+c_2^{}\frac{1}{\xi }q^{S_k^3}\end{array}\right),q=e^{i\alpha },\xi =e^{i\lambda }$$
(9)
or multiplying it by $`\sigma ^i,i=1,2,3`$ due to a symmetry of (8) as $`[R,\sigma ^i\sigma ^i]=0`$. The basic operators in (9) satisfy a quadratic quantum algebra

$$[S_k^+,S_l^{}]=4\delta _{kl}\mathrm{sin}\alpha (M^+\mathrm{sin}(\alpha 2S^3)+M^{}\mathrm{cos}(\alpha 2S_k^3)),[S_k^3,S_l^\pm ]=\pm \delta _{kl}S_k^\pm ,[M^\pm ,]=0,$$
(10)
which is a Hopf algebra with well defined coproduct and with $`c_a^\pm ,a=1,2`$ and $`M^\pm =\pm \frac{1}{2}\sqrt{\pm 1}(c_1^+c_2^{}\pm c_1^{}c_2^+)`$ being sets of central elements. Reductions of the general $`L`$-operator (9) (modulo its product with $`\sigma ^i`$), as shown in , generate in a systematic way the Lax operators of integrable UL models having the same $`R`$-matrix (8). At the undeformed $`q1`$ limit the whole procedure is repeated for the rational $`R`$-matrix.

We focus now on the NUL models under consideration for explicit solution of their UL sources. Only the simplest $`SU(2)`$ case is considered here; more general group admissible to some models can be treated analogously.
2.1. Coulomb gas picture of CFT

The NUL nature of this CGCFT model  follows from the Drinfeld-Sokolov problem $`_xQ=L(x)Q`$, since the linear operator

$$L(x)=v(x)\sigma ^3+\sigma ^+$$
(11)
contains a current-like field $`\{v(x),v(y)\}=i\alpha \delta ^{^{}}(xy).`$ Quantum and discrete version of this linear operator: $`L_k`$ was shown to satisfy the spectral parameter-less analogs of braided equations (1) and (2) with the braiding matrix $`Z`$ as presented below and the limit of (8): $`R(\lambda \mathrm{})R_q^+`$ acting as the $`R`$-matrix.

Note that by fixing the central elements in (9) as $`c_1^+=\mathrm{\Delta }`$, with $`\mathrm{\Delta }`$ being the lattice constant and all the rest $`c`$’s as zeroes, the algebra (10) reduces simply to

$$[S_k^+,S_l^{}]=0[S_k^3,S_l^\pm ]=\pm \delta _{kl}S_k^\pm $$
(12)
and we may find a realization for the generators as

$$S_k^{}=e^{\pm ip_k},S_k^3=u_k\frac{1}{\alpha }p_k$$
(13)
reducing the Lax operator (9) (after multiplying it from right by $`\sigma ^1`$) to

$$L_k^{ul}(\xi )=e^{ip_k\sigma ^3}+\mathrm{\Delta }\xi e^{i(\alpha u_kp_k)}\sigma ^+,$$
(14)
in canonical variables $`[u_k,p_j]=i\delta _{kj}`$. We may conclude that (14) represents a relativistic Toda chain (RTC) like integrable UL model  and satisfies the YBE with the $`R`$-matrix (8), as a consequence of the general argument of . We observe that a local gauge matrix $`D_j^{(1)}=e^{\frac{i}{2}\alpha \sigma ^3u_j}`$ along with a change of basic operators:

$$v_k^+=\frac{\alpha }{2}(u_{k+1}u_k)p_k,v_k^{}=\frac{\alpha }{2}(u_{k+1}u_k)p_k$$
(15)
can transform (14) to

$$L_k^{nul}(\xi )=e^{iv_k^{}\sigma ^3}+\mathrm{\Delta }\xi e^{iv_k^+}\sigma ^+,$$
(16)
NUL properties of which are induced by the discrete version of the current algebra for $`v_k^\pm `$:

$$[v_k^\pm ,v_l^\pm ]=\pm i\frac{\alpha }{2}(\delta _{k,l+1}\delta _{k+1,l}),[v_k^+,v_l^{}]=i\frac{\alpha }{2}(\delta _{k+1,l}2\delta _{k,l}+\delta _{k,l+1}).$$
(17)
One may check that (14) as well as (16) satisfy (4) with the braiding matrix $`Z^{(1)}=e^{\frac{i}{2}\alpha \sigma ^3\sigma ^3}`$, as found in . This guarantees in turn that the NUL integrability relations (1), (2) and (3) must hold for the Lax operator (16). We observe further that at the continuum limit $`\mathrm{\Delta }0,`$ $`v_k^\pm \sqrt{\pm }1\mathrm{\Delta }v(x)`$ yield the current-like field and consequently (16) reproduces through $`L_k^{nul}(\xi )=I+\mathrm{\Delta }L(x,\xi )+O(\mathrm{\Delta }^2)`$, $`L(x,\xi )=v(x)\sigma ^3+\xi \sigma ^+`$, i.e. the field Lax operator (11) of , though generalized here to include spectral parameter $`\xi `$ and satisfy the standard YBE, which shows also its novelty. Thus (14) is proved to be the required UL source solution for the NUL integrable CGCFT model.
2.2. NUL sine-Gordon like model

For introducing a new discrete NUL sine-Gordon (NSG) type model, we start from an UL Lax operator

$$L_k^{ul}(\lambda )=e^{i(\lambda \sigma ^3p_k)}+e^{i(\alpha u_kp_k)}\sigma ^++e^{i(\alpha u_k+p_k)}\sigma ^{},\xi =e^{i\lambda }.$$
(18)
Note that like the previous model it can again be derived from the same general Lax operator (9), though now with a different parameter choice: $`c_1^{}=c_2^+=0,c_1^+=c_2^{}=1,`$ which reduces the underlying general algebra (10) to an interesting exponentially deformed algebra

$$[S_k^+,S_l^{}]=2i\delta _{kl}e^{2i\alpha S_k^3}\mathrm{sin}\alpha ,[S_k^3,S_l^\pm ]=\pm \delta _{kl}S_k^\pm .$$
(19)
Note that such reductions, which must represent integrable UL models associated with the trigonometric $`R`$-matrix (8) fall in the class of discrete Liouville model  and with a realization of (19):

$$S_k^{}=e^{\pm i(\alpha u_kp_k)},S_k^3=\frac{1}{\alpha }p_k$$
(20)
we can derive (18) directly from (9).

Remarkably, in spite of the different underlying algebra and its realization from the previous model, (18) permits us to choose the gauge operator again in the same form as the previous one: $`D_j^{(2)}D_j^{(1)}=e^{\frac{i}{2}\alpha \sigma ^3u_j}.`$ Performing therefore this gauge transformation on (18) and changing to the current-like operators $`v_k^\pm `$ using the map (15) we derive finally the Lax operator of our new NSG model as

$$L_k^{nul}(\lambda )=\left(\begin{array}{c}\xi W_k^+W_k^+\\ W_k^{}\frac{1}{\xi }W_k^{}\end{array}\right),W_k^\pm =e^{iv_k^\pm },\xi =e^{i\lambda }.$$
(21)
The validity of (4) by both (18) and (21) gives the explicit form of the braiding matrix $`Z_{12}^{(2)}=e^{i\frac{\alpha }{2}\sigma ^3}I`$, which differs clearly from the previous case. Therefore, as established above, this discrete NSG (21) must be a genuine quantum integrable NUL model, which satisfies (1) and the braided YBE (2). For finding the full set of its mutually commuting conserved operators: $`c_{\pm j},j=1,2.\mathrm{},N`$ we have to factorize the trace identity, which however becomes easy for the present $`Z`$ matrix yielding $`\stackrel{~}{\tau }(\xi )=tr(e^{i\frac{\alpha }{2}\sigma ^3}T^{nul}(\xi ))`$ as the generating function: $`\stackrel{~}{\tau }(\xi )=_{\pm j}^Nc_{\pm j}\xi ^{\pm j}`$. We may define the Hamiltonian of this new integrable model as $`H=\frac{1}{2}((c_N)^1c_{N2}+(c_N)^1c_{(N2)})=_j^N\mathrm{cos}(v_j^+v_j^{}+\alpha )`$ and name it as nonultralocal SG model, due to its resemblance with the standard SG. Note that noncanonical commutation relations (17) involving different sites, which must be used for deriving the dynamical equations, induce nearest-neighbor interaction in the model. Notice also that by transforming $`\xi \frac{1}{\xi },W_k^+W_k^{}`$ in (21) we would get a dual though similar integrable NSG model.
2.3. Quantum light-cone sine-Gordon model

Since quantum formulation through the braided YBE is not available for this well known nonultralocal classical model, our first goal is to discover such a formulation, for which we follow our above strategy and look for its UL source model as a particular realization of the general UL model (9). We observe that the choice of central elements as $`c_a^{}=0,c_a^+=\mathrm{\Delta }`$ or its complementary set $`c_a^{}=\mathrm{\Delta },c_a^+=0,`$ with $`a=1,2`$ leads again to the simple algebra (12). Therefore we can use realizations similar to (13):

$$S_k^{}=e^{\pm ip_k},S_k^3=u_k+ϵ\frac{1}{\alpha }p_k$$
(22)
with $`ϵ=`$ for the two complementary (left , right) cases and generate from (9) (after multiplying with $`\sigma ^1`$) the left ( right ) Lax operator $`L_k^{()}(\lambda )`$ as

$$L_k^{()ul}(\lambda )=e^{ip_k\sigma ^3}+\mathrm{\Delta }\xi ^{\pm 1}\left(e^{i(\pm \alpha u_kp_k)}\sigma ^++e^{i(\alpha u_k+p_k)}\sigma ^{}\right),\xi =e^{i\lambda },$$
(23)
representing again RTC type integrable UL model. Acting now, for example, on the left Lax operator $`L_k^{()}(\lambda )`$ of (23) by a gauge matrix $`D_j^{(3)}(\alpha )=e^{i\alpha u_j\sigma ^3}`$, which may be noticed to be simply the square of that found in the previous cases: $`D_j^{(3)}(\alpha )=(D_j^{(2)})^2`$, we construct a NUL Lax operator

$$L_j^{()lcsg}(\lambda )=e^{i(p_j\alpha u_j)\sigma ^3}+\mathrm{\Delta }\xi \left(e^{i(p_j+\alpha u_{j+1})}\sigma ^++e^{i(p_j+\alpha u_{j+1})}\sigma ^{}\right),u_ju_{j+1}u_j.$$
(24)
We may check directly that $`L_k^{()ul}(\lambda )`$ in (23) as well as $`L_j^{()lcsg}(\lambda )`$ in (24) respect condition (4) with the braiding matrix $`Z_{12}^{()(3)}=e^{i\alpha \sigma ^3\sigma ^3}`$ and therefore conclude that the NUL Lax-operator (24) together with $`R`$-matrix (8) and this $`Z^{()(3)}`$-matrix satisfies braiding relation (1) and the BYBEs (2), (3). This thus proves the genuine quantum integrability of the associated model which, as we see below, may be considered as the exact lattice version of the quantum light-cone SG (LCSG) model. Therefore, (23) indeed is an UL source model for the LCSG, which is new as a quantum integrable NUL model. At the field limit $`\mathrm{\Delta }0,`$ when we have $`p_j\mathrm{\Delta }_tu(x),\alpha u_ju(x),\alpha u_{j+1}u(x)+\mathrm{\Delta }_xu(x)`$, defining $`_tu\pm _xu=\frac{1}{2}_\pm u`$, it is not difficult to see that (24) reduces to $`L_j^{()lcsg}(\lambda )I+\mathrm{\Delta }\mathrm{U}_{}(x)`$, where $`\mathrm{U}_{}(x)=\frac{i}{2}_{}u(x)\sigma ^3+\xi (e^{iu(x)}\sigma ^++e^{iu(x)}\sigma ^{})`$ yields one of the Lax pair : $`_{}\mathrm{\Phi }=\mathrm{U}_{}\mathrm{\Phi }`$ of the well known LCSG field model.

It is important to note that the right operator $`L_k^{(+)}(\lambda )`$ in (23) can be obtained formally from the left $`L_k^{()}(\lambda )`$ through simple mapping $`\xi \frac{1}{\xi },\alpha \alpha `$ and therefore the corresponding results for the complementary right LCSG can be recovered easily from that of the left one through the same mapping. Thus we derive the NUL model $`L_j^{(+)lcsg}(\lambda )`$ which at the field limit gives similarly the other component of the LCSG Lax pair: $`\mathrm{U}_+(x)`$. Zero curvature condition: $`_{}\mathrm{U}_+_+\mathrm{U}_++[\mathrm{U}_+,\mathrm{U}_{}]=0`$ involving both these Lax operators yields finally the well known form of the sine-Gordon field equation in light-cone coordinates: $`_+^2u=2\mathrm{sin}2u`$.

It is intriguing to remark that, the UL source models found recently for the quantum left/right mKdV  can be given exactly by the same Lax operators (23) discovered here for the light-cone SG model, though the gauge operator related to the mKdV model coincides rather with $`\sqrt{D}_j^{(3)}(\alpha )`$.
2.4. Quantum mapping

A rational quantum mapping related to the lattice KdV type equation may be described by a Lax operator

$$U_n=L_{2n}^{nul}L_{2n1}^{nul},\text{where}L_{aj}^{nul}(\lambda _a)=\frac{1}{2}(1+\sigma _a^3)v_j+\sigma _a^++\lambda _a\sigma _a^{},a=1,2$$
(25)
which exhibits NUL nature due to the presence of discrete current-like operator $`v_jv_j^+`$ (17). $`L_j^{nul}`$ and hence $`U_n`$ in (25) were shown  to satisfy both braiding relation and BYBE with the rational $`R`$-matrix $`R_{12}(\lambda _1\lambda _2)=\mathrm{𝟏}+i\frac{\alpha }{2}\frac{P_{12}}{\lambda _1\lambda _2},`$ (which may be obtained from (8) at $`\alpha 0,\lambda 0`$) and a spectral parameter dependent braiding matrix $`Z_{21}^{(4)}(\lambda _2)=\mathrm{𝟏}\frac{i\alpha }{2\lambda _2}\sigma ^{}\sigma ^+`$ (similarly $`Z_{12}(\lambda _1)`$).

For finding the UL connection as asserted by the conjecture, we look again for the source model as a realization of the general construction and try to identify a suitable gauge matrix, as done in the above cases. However since the present model is associated with the rational $`R`$-matrix we have to start now from the ancestor Lax operator of quantum integrable rational UL models , which is obtained as a $`\alpha 0`$ (or $`q1`$) limit of (9) in the form

$$L^{ranc(ul)}(\lambda )=\left(\begin{array}{c}c_1^0(\lambda +s^3)+c_1^1s^{}\\ s^+c_2^0(\lambda s^3)c_2^1\end{array}\right).$$
(26)
The basic operators satisfy the undeformed quadratic algebra

$$[s^+,s^{}]=2m^+s^3+m^{},[s^3,s^\pm ]=\pm s^\pm $$
(27)
with $`m^+=c_1^0c_2^0,m^{}=c_1^1c_2^0+c_1^0c_2^1`$ and $`c_a^n,a=1,2,n=0,1`$ as central elements. For our construction we choose first these elements as $`c_1^0=i\frac{\alpha }{2},c_2^0=0,c_1^1=c_2^1=1`$ to derive the underlying algebra from (27) in a simple form

$$[s^+,s^{}]=i\frac{\alpha }{2},[s^3,s^\pm ]=\pm s^\pm .$$
(28)
We can therefore realize this algebra simply as

$$s_k^+=i\psi _k,s_k^{}=\varphi _k,s_k^3=i\frac{2}{\alpha }\varphi _k\psi _k$$
(29)
through operators $`[\psi _j,\varphi _k]=\frac{\alpha }{2}\delta _{jk}`$ to generate from (26) (after multiplying from left by $`\sigma ^1`$) an UL Lax operator

$$L_j^{ul}(\lambda _1)=\left(\begin{array}{c}i\psi _j1\\ \lambda _1i\varphi _j\psi _j\varphi _j\end{array}\right),\lambda _1=1+i\frac{\alpha }{2}\lambda .$$
(30)
Remarkably this Lax operator is associated with a quantum integrable simple lattice NLS model  and satisfies the YBE with rational $`R`$-matrix. We find next the gauge operator as $`D_j^{(4)}=I\varphi _j\sigma ^{}`$, use a simple realization of the current-like operator satisfying (17) as $`v_j^+=i\psi _j+\varphi _{j1}`$ and perform the gauge transformation on our UL model (30) to derive exactly the Lax operator $`L_j^{nul}`$ of the quantum mapping (25). Note however that since we intend to recover the known form of $`L_j^{nul}`$ as given in , we had to use here a slightly different form of gauge transformation: $`D_jL_j^{ul}(\lambda _1)D_{j1}^1=L_j^{nul}(\lambda _1)`$, compared to the rest of our cases. As a result $`L_j^{nul}`$ would satisfy NUL integrable relations, which were used in  and are complementary to those presented here and also in place of (4) a similar condition $`D_{aj}L_{bj}^{nul}D_{aj}^1=Z_{ab}^1L_{bj}^{nul}`$. This establishes therefore (30) to be a perfect gauge related UL source solution for the well known NUL quantum mapping.
3. Conclusion

We conjecture an ultralocal (UL) connection for every nonultralocal (NUL) quantum integrable model and establish the yet unestablished such connections including some new ones. We find explicitly for the well known quantum NUL systems, e.g. CGCFT  and quantum mapping , the corresponding gauge related UL source models as the relativistic Toda chain (RTC) and a lattice NLS, which are well studied UL models. At the same time we discover two new quantum integrable NUL models, namely a discrete NUL type SG (NSG) and an exact lattice light-cone SG (LCSG) model and identify their corresponding UL connections as given through a Liouville and a RTC model. Such an explicit UL gauge relation of a NUL model would help in extracting all necessary knowledge of a more complicated NUL system including its Bethe ansatz result through a more tractable and in most cases already known UL integrable model. The NUL-UL connection should become particularly useful, when we are unable to tackle a NUL integrable model directly, due to nonfactorizability of its trace identity or other reasons. In such cases, thank to the present conjecture, we could switch over to its gauge related UL model and use relation like (6): $`\stackrel{~}{T}^{nul}=T^{ul}`$ for its UL description. The $`\stackrel{~}{T}^{nul}=L_{aN}^{nul}\mathrm{}L_{a2}^{nul}\stackrel{~}{L}_{a1}^{nul}`$ however would be an approximation, where one of its boundary terms is modified as $`\stackrel{~}{L}_{a1}^{nul}=D_1^1L_{a1}^{nul}D_1`$. It is remarkable that, though we have considered NUL models of completely different nature with no apparent relations between them, their gauge related UL sources could be derived from the same ancestor model (or its $`q1`$ limit). Such universality of this ancestor is due to its general form, which is capable of generating all possible quantum integrable UL models within the specified class . Therefore we can conjecture further that all other quantum integrable NUL models, not studied here, also should have their gauge related UL models derivable from the same ancestor. Various reductions of this ancestor model, which generate different descendant UL models and serve as the source models for the NUL systems, are due to different choices of the central elements leading to varied underlying algebras and hence various allowed realizations. For generating the corresponding NUL models we need to find further the gauge operators, which depend in general on the structure of individual models as do also the associated braiding matrices. Among the four NUL models we considered, CGCFT and LCSG have the same underlying algebra, while CGCFT and NSG have the coinciding gauge operator. The braiding matrices of all these models however turn out to be different.

Since for our construction we exploit the general algebraic scheme of , we also unravel easily the important underlying algebraic structure for each of the NUL models. This algebraic information along with their gauge connections to UL models should make quantum integrable NUL models more interesting and easier objects to study, boosting their development which indeed deserves and needs more attention.