text:
# Infinite compressibility states in the Hierarchical Reference Theory of fluids. I. Analytical considerations

## 1 Introduction

Reconciling the vastly different approaches to fluid structure and thermodynamics afforded by classical integral equation (ie) formalisms and renormalization group (rg) theory, the Hierarchical Reference Theory (hrt, ) presents itself as a particularly effective instrument for studying the critical region and liquid-gas phase equilibrium in simple one-component fluids: For subcritical temperatures, $`T<T_c`$, the usual formulation of the theory  yields density intervals of rigorously flat free energy and infinite isothermal compressibility $`\kappa _T`$ the boundaries of which are readily identified with the densities $`\varrho _v`$ and $`\varrho _l`$ of the coexisting gas and liquid phases. The binodal so found terminates at some temperature $`T=T_c`$ and density $`\varrho =\varrho _c`$ in a liquid-gas critical point characterized by non-classical, partly Ising-like exponents . And far away from the coexistence region of the phase diagram, hrt reduces to one of the standard approximations of liquid state theory, viz., the popular scheme commonly known under the names of Lowest-Order $`\gamma `$ Ordered Approximation (loga, ) and Optimized Random-Phase Approximation (orpa, ). hrt’s unified treatment of thermodynamic states so diverse should be contrasted with the limitations inherent in approaches based on ies alone: Close to the critical point these generally do not have a solution or else develop major deficiencies, and the binodal is accessible only by way of a Maxwell construction . rg calculations alone, on the other hand, while invaluable for illuminating the scaling relations valid asymptotically close to the critical point, generally do not allow one to determine non-universal quantities such as, e. g., the loci of the critical point and the binodal. A theory like hrt that provides comprehensive structural and thermodynamic information both close to and away from the critical point and exhibits a non-trivial scaling limit is clearly attractive for applications calling for high-resolution data on the behavior of a fluid in the critical region.

In the present two-part series of reports we want to have a closer look at the solution of the hrt equations for thermodynamic states of diverging compressibility, i. e., at the critical point and at phase coexistence. Our motivation for this inquiry is twofold: First of all, we aim to extend our understanding of the way in which hrt achieves its remarkable description of criticality and phase separation beyond mere invocation of its conceptual ingredients, rg theory and thermodynamic consistency (v. i.) in particular. Instead, it is on the level of the partial differential equation (pde) itself and the corresponding finite difference (fd) approximations used in practical calculations that we want to understand the mechanism responsible for the suppression of van der Waals loops and the emergence of a singular limit of $`\kappa _T`$ in an extended part of the phase diagram. A second reason for our investigation lies in our earlier work on hrt and its numerical side : For all the merits of the theory, its practical application has been found to be troubled with two major difficulties that have been traced to the customary way of incorporating the core condition of vanishing pair distribution function for hard core reference systems, and to the numerical properties of the equations for high compressibility, respectively. The latter clearly is an issue of prime importance when focussing on phase separation and the immediate vicinity of the critical point where $`\kappa _T`$ diverges, and its severity can only be assessed on the basis of a thorough understanding of the numerical process in relation to the properties of the hrt pde. Evidently, such an understanding is also highly relevant to the interpretation of numerical results and the extraction of meaningful and reliable information from them.

In order to shed some light on these questions, in the present report we study the analytical properties of the pde in the presence of a singular limit of $`\kappa _T`$. Relegating some details to the appendix, after a short introduction to hrt itself and its conceptual basis in section 2 we present the pde, identify a quantity convenient for following the build-up of infinite compressibility, and infer the asymptotic scaling relations we base our work on. Employing a suitably formalized notion of smoothness, in sections 3 through 6 we find a total of three scenarios for the gradual build-up of infinite $`\kappa _T`$, clarify their respective preconditions, and infer some of their properties with a view to an eventual implementation by fd methods. According to their most prominent traits we refer to the classes of asymptotic solutions so found as “(genuinely) smooth” (section 3), “monotonous” (section 5), and “stiff” or only “effectively smooth” (section 6), respectively. Of these, only the first seems to be compatible with an earlier analysis of hrt’s scaling limit at first sight ; however, a closer investigation into the assumptions implicit in the simplifications made there leads to a reappraisal of those results that therefore cannot invalidate either of the remaining two candidate types of solution (section 7).

As the considerations outlined above do not take into account the initial and boundary conditions imposed on the pde but rather concern themselves with a summary analysis of the asymptotic behavior of the various terms in the pde and of the range of solution types compatible with these, the all-important question of which of the scenarios captures the true behavior cannot be answered in the present contribution. In part II of our investigation  we present a host of numerical evidence strongly suggestive of a pde asymptotically turning stiff for thermodynamic states of infinite compressibility while fd methods on practical discretization grids bring about a regularization that leads to an artificially smoothed solution. Assertion of stiffness also paves the way for a qualitative understanding of the relation between specific features of the Fourier transform of the potential and its numerical properties, clarifies the special standing of the hard-core Yukawa potential, and leads to a detailed and self-consistent perception of the process of solving the equations throughout the domain of the pde. Needless to say, these findings are highly relevant for the interpretation of numerical results and the methods of data analysis to be applied to them, especially when working with non-uniform high-resolution discretization grids.

## 2 Basic relations

As a starting point, let us shortly review the concepts underlying hrt when applied to simple one-component fluids, recalling some of its central notions and establishing the equations we will base our work on; for a more detailed account of the derivation, its physical justification and relation to both the ie formalism and rg theory as well as the modifications necessary for dealing with other physical systems, most notably spin models and fluid mixtures, we refer the reader to refs.  and further references therein. For consistency with our earlier work on hrt we employ a number of notational conventions summarized in the appendix: Most importantly, a tilde indicates Fourier transformation, and once a symbol has been introduced we generally omit obvious generic function arguments. The appendix also serves as a repository for some of the more cumbersome analytical expressions as well as for auxiliary definitions and relations tangential to our reasoning but necessary to make our presentation self-contained.

Working in the grand-canonical ensemble we consider a system of particles interacting via pair-wise additive forces taken to derive from a potential $`v(r)=v^{\mathrm{ref}}(r)+w(r)`$. For the sake of simplicity, the potentials are assumed $`\varrho `$ independent, and we restrict the reference fluid corresponding to $`v^{\mathrm{ref}}`$ alone to a system of hard spheres of diameter $`\sigma `$, i. e., $`v^{\mathrm{ref}}(r)`$ is infinite for $`r<\sigma `$ and vanishes otherwise. The perturbation $`w(r)`$ and temperature $`T`$ enter the calculation only in the combination $`\varphi (r)=\beta w(r)`$, where $`\beta =1/k_BT`$ and $`k_B`$ is Boltzmann’s constant.

Based on this splitting of $`v`$, a momentum space cutoff $`Q`$ is introduced by the device of a rather artificial  $`Q`$ dependent potential $`v^{(Q)}(r)`$ obtained from $`v(r)`$ by the elimination from $`w`$ of all Fourier components $`\stackrel{~}{w}(k)`$ with $`k<Q`$, i. e.,

$$\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}\end{array}$$

Clearly, the reference and target systems with potentials $`v^{\mathrm{ref}}`$ and $`v`$ are obtained in the limits $`Q\mathrm{}`$ and $`Q0`$, respectively. A rather intricate analysis of a resummed perturbation expansion for the properties of the $`(Q\mathrm{\Delta }Q)`$ system in terms of those of the $`Q`$ system at the same temperature and density in the limit $`\mathrm{\Delta }Q0`$ finally yields a non-terminating hierarchy of first-order ordinary differential equations (odes) in $`Q`$ for the free energy $`A^{(Q)}(\varrho )`$ and the $`n`$ particle direct correlation functions . Formally, these differential equations allow one to follow the evolution of structure and thermodynamics of the $`Q`$ systems when fluctuations of ever increasing wavelength $`1/Q`$ are taken into account, i. e., when $`Q`$ goes from infinity to zero and $`v^{(Q)}`$ is transformed from $`v^{\mathrm{ref}}`$ into $`v`$.

Such an infinite set of coupled odes is, of course, hardly tractable numerically, let alone analytically. As a remedy, a closure on the two-particle level resembling loga/orpa, eq. (A4) in the appendix, is customarily adopted and combined with only the first hrt equation giving the $`Q`$ dependence of $`A^{(Q)}`$. As demonstrated in ref. , if hrt’s ability to describe phase coexistence is not to be lost, it is vital to also incorporate a condition of thermodynamic consistency into the closure: In hrt’s standard formulation this takes the form of the compressibility sum rule (A5) relating $`\kappa _T^{(Q)}`$ as obtained by differentiation of the free energy to the volume integral of the direct correlation function at arbitrary $`Q`$. Due to the density derivatives so introduced the odes at fixed $`\varrho `$ give way to a single pde in $`Q`$ and $`\varrho `$ for the free energy that is to be solved on the semi-infinite strip $`𝒟`$ where $`\varrho _{\mathrm{min}}\varrho \varrho _{\mathrm{max}}\mathrm{}>Q0`$. The precise choice of the initial and boundary conditions that remain to be imposed at $`Q=\mathrm{}`$, at $`\varrho =\varrho _{\mathrm{min}}`$, and at $`\varrho =\varrho _{\mathrm{max}}`$ is of no importance for the remainder of this work, and we refer the reader to refs.  for a more detailed discussion of this point.

Discretization of the pde for the free energy so obtained is straightforward and yields a computational scheme that can, indeed, successfully be used for $`T>T_c`$; for close-to-critical and subcritical temperatures, however, attempts at a direct solution invariably fail to produce any results . In order to remedy this situation, Tau et al.  proposed an alternative formulation in terms of an auxiliary quantity $`f(Q,\varrho )`$ that is essentially the first $`Q`$ derivative of the free energy; just as in our previous work on hrt , in the present contribution we rely on a slightly different definition for $`f`$ detailed in the appendix, cf. eq. (A3). Further specializing to density-independent potentials and not explicitly including the core condition that is not expected to be relevant to the subject of our study, the pde can be written in quasilinear form,

$$\frac{f}{Q}=d_{00}[f;Q,\varrho ]+d_{02}[f;Q,\varrho ]\frac{^2f}{\varrho ^2},$$
(1)
with initial and boundary conditions that directly follow from those imposed in the original formulation.

The rather lengthy expressions for the coefficients $`d_{00}`$ and $`d_{02}`$ of eq. (1) are to be found in the appendix, as are the defining relations for a number of auxiliaries. Among these the quantity $`\epsilon (Q,\varrho )\overline{\epsilon }(Q,\varrho )+1`$ is of particular relevance to our reasoning: Essentially the exponential of $`f`$, it turns out proportional to the isothermal compressibility of the fully interacting system, cf. eq. (A6). Infinite $`\kappa _T`$ therefore directly implies attendant divergences at $`Q=0`$ in $`\epsilon `$, $`\overline{\epsilon }`$, and $`f`$, and we have to study the large-$`\overline{\epsilon }`$ behavior of the pde if we are to understand the description of the critical region afforded by hrt on the level of the pde.

On the other hand, $`f(Q,\varrho )`$ is guaranteed by the construction of the hrt hierarchy to be continuous and finite for every non-vanishing cutoff $`Q`$ and to coincide with its limit from above at $`Q0`$ wherever that limit exists. In other words, for thermodynamic states $`(T,\varrho )`$ within the coexistence part of the phase diagram $`f`$ must take on large finite values for sufficiently small non-vanishing $`Q`$, and it must diverge for $`Q0`$: If $`\kappa _T(\varrho )`$ is infinite, for every threshold $`F`$ there is a corresponding cutoff $`Q_F(\varrho )>0`$ such that $`f(Q,\varrho )>F`$ for all $`Q<Q_F(\varrho )`$, or, less formally,

$$\kappa _T(\varrho )=\mathrm{}\{\begin{array}{cc}lim_{Q0}f(Q,\varrho )=\mathrm{}\hfill & \\ f(Q,\varrho )<\mathrm{}:Q\sigma >0\hfill & \\ f(Q,\varrho )1:Q\sigma 1,\hfill & \end{array}$$
(2)
which is the very basis of the main arguments of the present report. Considering large $`\overline{\epsilon }(Q,\varrho )`$ and assuming $`f`$ to diverge more strongly than just logarithmically (v. i., section 5), inspection of the expressions in the appendix shows both of the coefficients of the quasilinear pde (1) to be of first order in $`\overline{\epsilon }`$ whereas $`f`$ is essentially the logarithm of $`\overline{\epsilon }`$:

$$\begin{array}{cc}\hfill d_{02}& =𝐎(\overline{\epsilon }),\hfill \\ \hfill d_{00}& =𝐎(\overline{\epsilon }),\hfill \\ \hfill f& =𝐎(1).\hfill \end{array}$$
(3)
$`𝐎`$ and its companion $`𝐨`$ are the usual Landau symbols, and here as in the remainder of the present series of reports $`x=𝐎(y^a)`$ is taken to actually mean that $`a`$ is the infimum of all $`b`$ for which $`x=𝐨(y^b)`$ when the implied limit is $`y\mathrm{}`$, or else the supremum when considering $`y0`$. Furthermore, qualification of some quantity $`x`$ as “essentially” independent of some other quantity $`y`$ is equivalent to the characterization of $`x`$ as of order $`𝐎(1)`$ in $`y`$ which does not rule out a weak, say, logarithmic $`y`$ dependence of $`x`$. Unless explicitly stated otherwise, all orders cited are in terms of $`\overline{\epsilon }`$.

In order to understand the properties of the solution of the pde around some point $`(Q,\varrho )𝒟`$ we still need to supplement eq. (3) with an analogous characterization of the behavior of the remaining term on the right hand side of eq. (1), viz., $`^2f/\varrho ^2`$. Lacking any a priori information to guide us, in this report we adopt the simple ansatz

$$\frac{^2f}{\varrho ^2}=𝐎(\overline{\epsilon }^r),r0,$$
(4)
a choice that is sufficiently general to admit a consistent description both of the behavior of the exact solution and of the computational process yielding a numerical approximation of it. Inserting eqs. (3) and (4) into the pde (1), comparison of both sides of the equation immediately yields

$$\begin{array}{c}\frac{f}{Q}=𝐎(\overline{\epsilon }^s),\\ s=\mathrm{max}(1,1+r)=1+r;\end{array}$$
(5)
only for $`r=0`$ is there the added possibility of a cancellation of the leading terms on the right hand side of eq. (1) that might give rise to an $`s`$ less than unity, including $`s=r=0`$. At any rate, neither $`r`$ nor $`s`$ may be negative.

Equipped with eqs. (3) through (5) we are now in a position to gain a better understanding of the pde’s workings. Although the analytical expressions presented and orders cited below hinge on the assumption (4), most of our conclusions are expected to remain qualitatively valid even in more general situations, cf. section V in part II .

## 3 Genuinely smooth solution

From the pde itself the $`Q`$ and $`\varrho `$ scales characteristic of the variations of $`f(Q,\varrho )`$ in the integration domain $`𝒟`$ are not apparent: In particular, we cannot say a priori whether they remain bounded from below in an essentially $`\overline{\epsilon }`$ and, hence, temperature independent way throughout $`𝒟`$ or else scale like some inverse power at least of $`\overline{\epsilon }`$ and so become arbitrarily small in part of $`𝒟`$ whenever the temperature falls below $`T_c`$. Restricting ourselves to analytical considerations, in this report we will study and present arguments for both of these types of solution, referring to them as smooth and non-smooth, respectively. Clearly, this distinction is highly relevant to the numerics and therefore of immediate practical interest: After all, smoothness implies that finite difference (fd) schemes are in principle well applicable to the pde at hand whereas otherwise local truncation errors will generally be unbounded and at least an estimate of the global error incurred must be obtained a posteriori in order to gauge the significance of any information extracted from fd approximations of the pde.

For this section turning to the assumption of $`f`$ being smooth in the sense stated — an assumption that we will repeatedly refer to as leading to the “(genuinely) smooth” scenario —, we immediately conclude that $`r=s=0`$: If the $`Q`$ and $`\varrho `$ scales set by $`f`$ are bounded from below in the manner indicated, for any linear differential operator $`L`$ we can choose $`\overline{\epsilon }`$ independent step sizes $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$ such that an fd approximation of $`Lf`$ becomes accurate throughout $`𝒟`$. As the estimate so obtained is just a linear combination of $`f`$ values sampled few $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$ apart, $`Lf`$ is bound to scale like $`f`$ and is thus of order $`𝐎(1)`$. Specializing to $`L^2/\varrho ^2`$ and to $`L/Q`$ we immediately conclude that $`r=0`$ and $`s=0`$, respectively. — As a corollary we note that growth of $`f`$ in proportion to an inverse power of the cutoff, $`f1/Q^a`$ at fixed $`\varrho `$ with $`a>0`$, always falls into this class of solutions.

The main advantage afforded by the presupposition of smoothness is that it allows us to understand on the level of the pde the build-up of infinite $`f`$ and $`\kappa _T`$ close to the critical point and in the coexistence region of the phase diagram: Let us assume that, at some fixed and sufficiently small cutoff $`Q`$, $`f(Q,\varrho )`$ is a function of $`\varrho `$ like that sketched in fig. 1: continuously differentiable for all densities, convex and large for densities in a range with approximate boundaries $`\varrho _1`$ and $`\varrho _2`$ but rather small elsewhere. In the region of large $`f`$ and $`\overline{\epsilon }`$ where asymptotic reasoning along the lines of eqs. (3) through (5) is applicable, the pde coefficient $`d_{00}`$ is dominated by terms related to the $`Q`$ dependence of the Fourier transforms of the potential and of the reference system direct correlation function,

$$d_{00}\frac{(\stackrel{~}{\varphi }/\stackrel{~}{𝒦})}{Q}\frac{\stackrel{~}{𝒦}^2\overline{\epsilon }^2}{\epsilon }\frac{\stackrel{~}{\varphi }_0^2}{\stackrel{~}{\varphi }^4}\text{for }\stackrel{~}{\varphi }\epsilon 1,$$
(6)
cf. eq. (A1). $`d_{00}`$ is therefore expected to be negative for most cutoffs less than the position of the first minimum of $`\stackrel{~}{\varphi }`$, and certainly for very small $`Q`$: It is by this standard that the cutoff under consideration must be “sufficiently small” as stated before. — Trivially, we also see from eq. (A1) that $`d_{02}`$ is non-positive throughout $`𝒟`$, and negative for $`\overline{\epsilon }0`$, as $`\stackrel{~}{\varphi }_0`$ must be positive for the pde to be stable at all, cf. section 2.4.1 of ref.  as well as refs. .

As for the remaining term on the right hand side of eq. (1), $`^2f/\varrho ^2`$ is negative throughout most of the density interval of large $`f`$, and positive on either side as $`f`$ falls to small values. Furthermore, the maximum of $`f`$ displayed in fig. 1 is rather flat so that $`|^2f/\varrho ^2|`$ is quite small there. Assuming this curvature so small that the $`d_{02}`$ term in eq. (1) does not reverse the sign of the right hand side, we immediately conclude that $`f/Q<0`$ where $`f`$ is large, corresponding to further growth of $`f`$ as $`Q=0`$ is approached from above. In this way we can understand how the pde implements the gradual build-up of infinite $`\kappa _T`$ and the attendant suppression of van der Waals loops in the free energy. On the other hand, it becomes clear that this mechanism depends on the preservation of the general outline of $`f`$: In particular, the stable growth of $`f`$ is likely to break down once flatness is lost in the central part of the interval $`[\varrho _1,\varrho _2]`$. A closer look at the $`d_{02}`$ term immediately reveals that it acts to stabilize this feature of the form of $`f`$: Taken by itself, $`d_{00}`$ strongly favors $`f`$ to grow most rapidly close to its maximum. This, however, is prevented by the considerable increase in the negative density curvature of the solution it would entail: The correspondingly large contribution $`d_{02}\left(^2f/\varrho ^2\right)>0`$ to $`f/Q`$ effectively counteracts $`d_{00}`$ and so ensures that $`f(Q,\varrho )`$, $`\varrho _1\varrho \varrho _2`$, remains flat even as it grows, just as postulated at the outset. Close to $`\varrho _1`$ and $`\varrho _2`$, on the other hand, the rôle of the $`d_{02}`$ term is quite different: There the curvature $`^2f/\varrho ^2`$ turns positive so that both terms on the right hand side of eq. (1) now contribute to the growth of $`f`$, thereby rendering the transition to small $`f`$ ever more sharply defined as $`Q`$ progresses towards zero. By the same token, intermediate minima in the density range of large $`f`$ are also dissolved by the $`d_{02}`$ term, which relaxes the precondition of conformance with fig. 1 somewhat.

This intuitively appearling concept of a stable mechanism of growth of $`f`$, as we will henceforth refer to it, so allows us to understand on the level of the pde the emergence of infinite $`\kappa _T`$ within a clearly demarcated density interval $`[\varrho _v,\varrho _l]`$. In addition, as the coexisting densities $`\varrho _v`$ and $`\varrho _l`$ are obtained as the limits of $`\varrho _1`$ and $`\varrho _2`$ for $`Q0`$, anything but continuous inverse compressibility is hard to accomodate within this picture, which agrees well with the known coincidence of the hrt binodal and spinodal for space dimensionality $`d=3`$ . The key rôle played by the $`d_{02}`$ term in stabilizing the solution and locating the densities of the pure phases once more highlights the importance of the thermodynamic consistency condition (A5) underlying the transition from odes at fixed $`\varrho `$ to a pde over all of $`𝒟`$: Indeed, it should come as no surprise that application of an approximation incompatible with the afore-mentioned compressibility sum rule yields pathological results including negative compressibility .

Without evidence to the contrary it is then tempting to subscribe to the seemingly natural assumption of a smooth solution, $`r=s=0`$: Not only does this provide the basis for understanding the most salient features of hrt as just discussed, it is also what has been found in an early analysis of hrt’s scaling limit , cf. section 7 below. Unfortunately, however, this view is not entirely unproblematic as $`r`$ and $`s`$ must both vanish, v. s.: As stated in connection with eq. (5), this directly implies that the right hand side of eq. (1) is affected by massive cancellation so that the sum of two terms of order $`𝐎(\overline{\epsilon })`$ each is reduced to order $`𝐎(1)`$. Indeed, the relevant signs — $`d_{02}0`$ throughout $`𝒟`$ \[eq. (A1)\], $`d_{00}<0`$ for small $`Q`$ \[eq. (6)\], and $`^2f/\varrho ^2<0`$ for most densities in $`[\varrho _1,\varrho _2]`$ (fig. 1) — do allow such a cancellation. On the other hand, the rapid growth of the coefficient functions $`d_{0i}=𝐎(\overline{\epsilon })`$ places rather stringent constraints on the shape of $`f`$ at constant cutoff: According to eqs. (1) and (3) to (5), $`f=𝐎(1)`$ must deviate by terms of order $`𝐎(1/\overline{\epsilon })`$ at most from an exact solution $`\widehat{f}(Q,\varrho )`$ of the non-linear ode

$$\frac{^2\widehat{f}}{\varrho ^2}=\frac{d_{00}[\widehat{f};Q,\varrho ]}{d_{02}[\widehat{f};Q,\varrho ]}$$
(7)
for $`\varrho _1<\varrho <\varrho _2`$. Let us now suppose that the initial and boundary conditions imposed on the hrt pde actually lead to such a near-solution of the above ode at some sufficiently small $`Q`$, assessment of the likelihood of which falls outside the scope of this report. In this case consistency with $`s=0`$ requires the residue $`f\widehat{f}`$ of eq. (7) to vanish as $`1/\overline{\epsilon }`$ when $`Q`$ further progresses towards zero, but we have not been able to support this behavior of the solution on the basis of the explicit expressions (A1) for the $`d_{0i}`$ and their properties. Whereas the cancellation necessary for the genuinely smooth scenario thus certainly cannot be ruled out, it poses far stricter preconditions than the stable mechanism of growth discussed before, and both its genesis and stability remain unclear at this point.

## 4 Refined analysis

Looking back at the preceding section, we see that the strong point of genuine smoothness, viz., the stable mechanism of growth with its consequences, does not explicitly depend on $`s=0`$, the very condition that is the source of the conceptual difficulties just discussed. It thus seems pertinent to attempt to salvage these features of the smooth scenario in a less restricted setting by eliminating the assumption of vanishing $`s`$. But rather than separately discussing the case of $`r=0`$, $`s=1`$ we now turn to the full generality afforded by eqs. (4) and (5). Of course, for $`^2f/\varrho ^2`$ to be of order $`𝐎(\overline{\epsilon }^r)`$ with non-vanishing $`r`$, $`f`$ must be a rapidly oscillating function of $`\varrho `$ and we can only hope to adapt our findings if $`f`$ resembles fig. 1 when averaged over oscillations.

Let us consider once more the evolution of $`f(Q,\varrho )`$ within the region of large $`\overline{\epsilon }`$ as $`Q`$ progresses towards zero: Combining eq. (5) with the definition (A2) of $`\overline{\epsilon }`$ we find

$$\begin{array}{c}\frac{f}{Q}=𝐎(\overline{\epsilon }^s)=𝐎(e^{fs^{}}),\\ s^{}=s\stackrel{~}{u}_0^2.\end{array}$$

Here the modified exponent $`s^{}`$ takes into account the deviation of $`\stackrel{~}{u}_0(Q)\stackrel{~}{\varphi }(Q)`$ from unity; for $`Q0`$, $`s`$ and $`s^{}`$ coincide by virtue of the normalization condition $`\stackrel{~}{u}_0(0)=1`$. We now define an auxiliary quantity

$$\begin{array}{cc}\hfill d_0(Q,\varrho )& =\frac{d_{00}[f;Q,\varrho ]+d_{02}[f;Q,\varrho ]\left(^2f/\varrho ^2\right)}{e^{fs^{}}}\hfill \\ & =𝐎(1).\hfill \end{array}$$
(8)
Existence of $`d_0`$ is guaranteed, and the pde (1) is fully equivalent to

$$\frac{f}{Q}=d_0e^{fs^{}},$$
(9)
which no longer involves a derivative with respect to $`\varrho `$. The signs of $`d_0`$ and $`f/Q`$ always coincide.

If only for the moment we assume $`s^{}`$ constant, eq. (9) can be formally integrated to yield the solution at all $`Q`$ given the initial condition that $`f`$ be $`f_1`$ at cutoff $`Q_1`$, viz.,

$$f(Q,\varrho )=\frac{1}{s^{}}\mathrm{ln}\left(e^{f_1s^{}}+s^{}_Q^{Q_1}d_0(q,\varrho )dq\right),$$
(10)
which is valid only if $`s^{}s\stackrel{~}{\varphi }^2`$ does not vanish, i. e., if $`s>0`$ and $`\stackrel{~}{\varphi }0`$. The latter condition effectively restricts our arguments to cutoffs somewhat below $`Q_{\stackrel{~}{\varphi },1}`$, the smallest of the zeros $`Q_{\stackrel{~}{\varphi },i}`$ ($`i=1,2,\mathrm{}`$, $`Q_{\stackrel{~}{\varphi },i}<Q_{\stackrel{~}{\varphi },i+1}`$) of the Fourier transform $`\stackrel{~}{\varphi }`$ of the potential; in the limit $`Q0`$ that we are interested in this certainly poses no problem. With this proviso eq. (2) can be re-cast in form of the relations

$$\begin{array}{cc}_{Q_1}^Qd_0(q)dq<\frac{e^{fs^{}}}{s^{}}& \text{for }0<Q_1<Q,\hfill \\ _0^Qd_0(q)dq=\frac{e^{fs^{}}}{s^{}}& \text{for }\kappa _T=\mathrm{};\hfill \end{array}$$
(11)
for the thermodynamic states of interest both sides of the above inequality are positive if $`QQ_1`$ is comparable to $`Q`$.

The remainder of our analytical considerations on the nature of the pde for low cutoff in the critical region will be based on eq. (11). As $`e^{fs^{}}/s^{}`$ is strictly monotonous in $`s^{}`$, interval arithmetic is trivial and allows us to tackle in a straightforward way the problem of the $`Q`$ dependence of $`s^{}`$: Given the $`\overline{\epsilon }`$ independent and rather slow variation of $`s^{}/s\stackrel{~}{\varphi }^2`$ as a function of $`Q`$, for any cutoff interval considered the range of $`s^{}`$ values is easily found and translated into an interval of $`e^{fs^{}}/s^{}`$. Arguments based upon eq. (11) like those we will present shortly are thus easily modified to take into account the non-constancy of $`s^{}`$ simply by applying the least restrictive bound for any of the $`s^{}`$ values in the $`Q`$ interval under consideration, an operation taken for granted throughout the remainder of this series of reports.

## 5 Monotonous growth and logarithmic singularity

According to its definition (8), the integrand $`d_0`$ on the left hand side of eq. (11) is always of order $`𝐎(1)`$. On the other hand, for $`f`$ to be a monotonous function of $`Q`$ in the asymptotic region, its reduced slope $`d_0`$ must always be negative there so that $`_{Q_1}^Qd_0(q)dq=𝐎(1)(QQ_1)`$. As $`Q_1`$ is an essentially free parameter to be chosen from $`(0,Q]`$, eq. (11) shows a quantity scaling like $`Q`$ to be bounded from above by another one of order $`𝐎(\overline{\epsilon }^s)`$, i. e., $`𝐎(Q)<𝐎(\overline{\epsilon }^s)`$. This prompts us to consider a power law relation between $`Q`$ and $`\overline{\epsilon }`$, say,

$$\begin{array}{c}Q=𝐎(\overline{\epsilon }^t),\\ ts,\end{array}$$

corresponding to only a logarithmic singularity of $`f`$ at $`Q=0`$. In this case, eq. (3) no longer holds due to the cutoff dependence of the prefactors in eq. (A1): To leading order we find

$$\begin{array}{cc}\hfill d_{02}& =𝐎(\overline{\epsilon }^{12t}),\hfill \\ \hfill d_{00}& =𝐎(\overline{\epsilon }^{1t}),\hfill \\ \hfill f& =𝐎(1)\hfill \end{array}$$

instead, and the balance equation is modified to

$$𝐎(\overline{\epsilon }^s)=𝐎(\overline{\epsilon }^{1t})+𝐎(\overline{\epsilon }^{12t})𝐎(\overline{\epsilon }^r).$$
(12)
If the leading terms on the right hand side of eq. (12) do not cancel, eq. (5) is thus to be replaced by the relation $`s=\mathrm{max}(1t,12t+r)`$ with $`r0`$ and $`ts0`$. For $`r>t`$, however, the assumption of monotonous growth of $`f`$ leads to inconsistencies: In this case, the right hand side of the pde (1) is asymptotically dominated by the $`d_{02}`$ term. As rapid oscillations in the $`\varrho `$ direction at fixed $`Q`$ provide the only way for $`^2f/\varrho ^2`$ to become exponentially large compared to $`f`$ itself, the second density derivative of $`f`$ is bound to be oscillatory in the density range $`\varrho _1<\varrho <\varrho _2`$ just as well. Unlike $`f`$, however, $`^2f/\varrho ^2`$ must change its sign at every swing, which immediately carries over to $`f/Q`$ and $`d_0`$ due to the asymptotic dominance of the $`d_{02}`$ term, contrary to the assumption of $`d_0(q)`$ being negative for all $`q<Q`$.

We are thus left with the possibility of $`rt`$, $`s=1t`$; combining this with the conditions $`r0`$, $`s0`$, and $`ts`$ and eliminating $`t`$ we find the admissible exponent ranges

$$\begin{array}{ccc}\hfill 0& r& 1,\hfill \\ \hfill 0& s& \mathrm{min}(\frac{1}{2},1r)\frac{1}{2}.\hfill \end{array}$$

For $`r=t`$, there is the added possibility of cancellation of the leading terms on the right hand side of eq. (12); again eliminating $`t`$ from the relations $`0r=t`$, $`0s1t`$, and $`ts`$ we obtain

$$\begin{array}{ccc}\hfill 0& r& 1,\hfill \\ \hfill 0& s& \mathrm{min}(r,1r)\frac{1}{2}.\hfill \end{array}$$

At any rate, non-vanishing exponent $`s`$ is compatible with monotonous growth of $`f`$ only in the case of a merely logarithmic singularity at $`Q=0`$, and in this case there is an upper bound of $`\frac{1}{2}`$ for $`s`$ that holds irrespective of whether the leading terms on the right hand side of eq. (12) cancel. Furthermore, in this monotonous growth scenario $`𝐎(Q)<𝐎(\overline{\epsilon }^s)`$, v. s., so that finally $`Q^2\overline{\epsilon }`$ must tend to a finite, possibly vanishing limit for $`Q0`$.

## 6 Effective smoothness from stiffness

An alternative is to give up the monotonicity assumption and allow $`d_0`$ to alternate in sign: This opens up the possibility of a partial cancellation of the positive and negative contributions to the integral in eq. (11), and the average of $`d_0`$ no longer has to be of order $`𝐎(1)`$ even though $`d_0`$ itself still is. As we will now show, this situation implies that the hrt pde turns stiff for thermodynamic states of infinite compressibility, with $`r>0`$ and $`s>1`$. In numerical applications, on the other hand, discretization of the pde on practical grids by necessity induces an artificial smoothing of the numerical solution that is then characterized by vanishing effective exponents $`r_{\mathrm{eff}}`$ and $`s_{\mathrm{eff}}`$ provided the computation is able to reach the limit $`Q0`$ at all.

Stiffness of the pde in that part of $`𝒟`$ where the divergence of $`\kappa _T`$ is built up follows from eq. (11): In the scenario outlined above, the bound $`e^{fs^{}}/Qs^{}=𝐎(\overline{\epsilon }^s)/Q`$ on the mean of $`d_0`$ over the interval $`(0,Q]`$ implies that $`d_0`$ is a rapidly oscillating function of $`Q`$ both amplitude and period of which scale like $`1/\overline{\epsilon }^sQ=𝐎(\overline{\epsilon }^s)`$ \[the possibility of a non-monotonous logarithmic divergence, $`Q\overline{\epsilon }^t`$ with $`0<t<s`$, may be accounted for by replacing $`s`$ with $`st>0`$\]. Considering eq. (9), for $`\overline{\epsilon }f1`$ the solution is then characterized by oscillations with similar scaling properties superimposed upon a large, smooth, and most likely monotonously growing regular part. These oscillations in the $`Q`$ direction are naturally accompanied by oscillations in the $`\varrho `$ direction on density scales of order $`𝐎(\overline{\epsilon }^{r/2})`$, $`r>0`$: Vanishing $`r`$ would require the evolution of $`f`$ in $`Q`$ at densities at fixed separation to remain synchronous over many, viz., on the order of $`Q\overline{\epsilon }^s`$ undulations of the function; this seems highly unlikely, and $`r=0`$ must be unstable in the setting under discussion. From eq. (5) we then obtain $`r>0`$ and $`s=r+1>1`$, where the possibility of a weak $`Q`$ and $`\varrho `$ dependence of $`r`$ and $`s`$ must also be anticipated.

The scenario just discussed differs from genuine smoothness and monotonous growth in two important ways: Firstly, the singularity of $`f`$ may be stronger than before as there is now no need for the cancellation required in section 3, nor does boundedness of the mean of $`d_0`$ imply boundedness of the mean slope of $`f`$ as in section 5. Secondly, non-vanishing exponents $`r`$ and $`s`$ mean that the scales characteristic of the variation of $`f`$ and, hence, the grid spacings $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$ appropriate for the discretization of the pde become arbitrarily small as diverging isothermal compressibility is built up, which obviously has grave repercussions for the applicability of fd methods to the pde at hand.

This last point is worth considering in some detail: Clearly, for the fd equations to be a good approximation of the pde the step sizes must go to zero as appropriate inverse powers of $`\overline{\epsilon }`$ determined by $`r`$, $`s`$, and the orders of the local truncation error in $`\mathrm{\Delta }Q`$ and $`\mathrm{\Delta }\varrho `$. Excluding a very weak divergence of $`f`$, however, the rapid growth of $`\overline{\epsilon }`$ certainly renders so fine a grid impractical , and below some cutoff $`Q_{\mathrm{\Delta }Q}(T,\varrho )>0`$ the numerics can no longer follow the evolution in $`Q`$ of the true solution of the pde, nor that in $`\varrho `$ below some $`Q_{\mathrm{\Delta }\varrho }(T,\varrho )>0`$. Except right at $`\varrho _1`$ and $`\varrho _2`$, both $`Q_{\mathrm{\Delta }x}`$ are rather well defined because of the non-zero exponents $`s>r>0`$. Below the $`Q_{\mathrm{\Delta }x}`$, the basic assumption underlying any fd method, viz., applicability of Taylor expansions of rather low order, tantamount to smoothness of the solution on the scales set by the discretization grid spacings, no longer holds so that a discretization derived on this basis is used outside its range of validity. A chosen implementation of the theory may therefore fail to proceed to $`Q<Q_{\mathrm{\Delta }x}`$ altogether due to inappropriately large step sizes: In this case no solution is ever obtained at $`Q=0`$ below the critical temperature and neither criticality nor phase separation can be described, just as reported in ref.  and appendix B.1 of ref. . Alternatively, judicious choice of the formulation of the theory and of its discretization may allow the numerical scheme to solve the fd equations all the way to vanishing cutoff; it is this very property of eq. (1) that prompted adoption of the quasilinear form of the hrt pde in the first place, cf. section 2. Taylor expansion arguments no longer being applicable, however, a priori bounds on the local truncation error are not available, and the global error may well be substantial. Furthermore, as any solution generated in a fd calculation on a given grid is well represented on that same grid by definition, we can use the same arguments as in section 3 to show that the numerical results obtained with step sizes of order $`𝐎(1)`$ necessarily reproduce vanishing effective exponents $`r_{\mathrm{eff}}=s_{\mathrm{eff}}=0`$. This reduction of the exponents from $`s>r>0`$ for the exact solution corresponds to a smoothing of $`f`$ that presumably weakens the singularity, v. s., and also suppresses any oscillations on cutoff and density scales smaller than the step sizes. Such smoothing is very well known in the numerics of pdes and in fact forms the conceptual basis for the highly efficient multi-grid methods for the solution of integral and differential equations .

In this “stiff,” or only “effectively smooth” scenario, any results obtained by fd methods then necessarily realize the genuinely smooth case of section 3, admitting, i. a., a numerical solution $`f1/Q`$. In this case the artificial smoothing attendant to the reduction of the exponents directly affects the solution for $`TT_c`$ at any cutoff lower than $`Q_{\mathrm{smooth}}(T,\varrho )\mathrm{max}(Q_{\mathrm{\Delta }Q},Q_{\mathrm{\Delta }\varrho })`$, and the final results for $`\varrho _v\varrho \varrho _l`$. Given the $`d_{02}`$ term of the pde, however, such drastic qualitative ($`r_{\mathrm{eff}}r`$) and quantitative ($`s_{\mathrm{eff}}s`$) changes inside the region of large $`\overline{\epsilon }`$ must have an effect on the location of the binodal as well as on the solution for the pure phases, too. Restriction of large $`f`$ to only a rather well-defined density interval $`[\varrho _1,\varrho _2]`$, the pde’s parabolic character, and the boundary conditions imposed at $`\varrho _{\mathrm{min}}`$ and $`\varrho _{\mathrm{max}}`$ nevertheless inspire some hope that the error incurred outside the binodal might be limited. Indeed, in the presence of stiffness and smoothing this very hope must underlie any attempt at a numerical solution of the hrt equations and certainly has to be justified at least a posteriori by combining related calculations and testing the internal consistency of the results obtained .

## 7 Stiffness and the scaling limit of HRT

So far our considerations have led us to the formulation of three different possibilities for the asymptotic behavior of the hrt pde for high compressibility states characterized by true smoothness of the solution, by monotonous growth leading to only a logarithmic singularity, and by stiffness giving rise to effective smoothness of $`f`$, respectively.

At first sight, however, anything but the genuinely smooth solution, $`r=s=0`$, seems to be at variance with the detailed analysis of hrt’s scaling limit in ref. : In section III B of that report we find the prediction that the quantity denoted $`u_Q`$ there (corresponding to $`f`$) should be quadratic in the density-like variable $`x`$ and scale like $`Q^{(d2)}`$ in $`d`$ dimensions; in our notation, $`f`$ should thus grow like $`1/Q`$ in $`d=3`$ dimensions, corresponding to $`r=s=0`$, cf. section 3.

In order to resolve this seeming contradiction let us have a closer look at the reasoning of ref.  for $`T<T_c`$: The pivotal relation is eq. (3.9) there, viz.,

$$e^{u_Q}\frac{u_Q}{Q}=2Q\frac{1}{2}\frac{^2u_Q}{x^2}Q^{d1}.$$
(13)
For this equation, the direct analogue of our eq. (1), the authors of ref.  invoke the validity of “neglecting exponentially small terms” to justify replacing its left hand side by nought; the asymptotic solution proportional to $`Q^{(d2)}`$ then follows immediately by separation of the cutoff and density dependencies. But this replacement is legal only if the slope of $`u_Q`$ actually remains small in modulus compared to its exponential or, translating to our notation, if $`(1/\overline{\epsilon })(f/Q)`$ tends to zero as $`f`$ goes to $`+\mathrm{}`$, i. e., if $`s<1`$ as is the case for the solution cited where it actually vanishes. Conversely, eq. (13) can also be read to indicate that its left hand side will be small if and only if the pde’s initial and boundary conditions have caused $`u_Q`$ at the cutoff considered to almost exactly solve a simple ode in the density-like variable $`x`$, viz.,

$$\frac{^2u_Q}{x^2}=4Q^{2d},$$

which is completely equivalent to the cancellation requirement of section 3, eq. (7) in particular. The above directly shows that $`^2u_Q/x^2`$ is proportional to $`Q^{2d}u_Q`$ and so of order $`𝐎(1)`$ in $`e^{u_Q}`$; in terms of $`f`$ this means that $`^2f/\varrho ^2`$ is of order $`𝐎(1)`$ in $`\overline{\epsilon }`$ and, hence, that $`r=0`$.

As the assumptions both of smoothness in the $`\varrho `$ direction, $`r=0`$, and of cancellation of the leading terms on the right hand side of eq. (1), $`s<1`$, are thus built right into the derivation of the scaling solution of ref. , conformance of the latter with the smooth scenario is hardly surprising and certainly does not rule out any of the other possibilities considered in the present report. As for the numerics, both genuine ($`r=s=0`$) and effective ($`r_{\mathrm{eff}}=s_{\mathrm{eff}}=0`$) smoothness suggest an $`f`$ that asymptotically grows like $`1/Q`$, although in the latter case this does not reflect the properties of the exact solution of the pde. We are thus led to a reappraisal of the analysis presented in section III B of ref.  as certainly applicable to the fdes solved numerically but not to the true pde unless $`r=s=0`$, i. e., unless the cancellation requirement of section 3 is also met; with this proviso the conclusions of ref.  remain largely unaffected. As for the intuitively appealing stable mechanism of growth responsible for suppression of van der Waals loops, it now applies to all three scenarios, if only for the numerical process unless $`r=0`$; further support for it derives from the $`x`$-dependence of $`u_Q`$ that conforms with the sketch of fig. 1 whereas eq. (7) might well be compatible with a more general outline.

In view of this resolution of the seeming contradiction between the earlier analysis of hrt’s scaling limit and the possibility of non-zero exponents $`r`$ and $`s`$ none of the three scenarios can be excluded from further consideration at this point. On the other hand, the rather summary analytical considerations presented here can only identify the solution types consistent with the asymptotic properties of the pde for infinite compressibility and explore their preconditions and consequences, which has been the subject matter of the present report. The all-important question of which of them is in fact realized in actual applications of the theory, however, crucially depends on global aspects of eq. (1) in all of $`𝒟`$ and certainly cannot be answered without considering the influence of the initial and boundary conditions that, after all, uniquely determine $`f(Q,\varrho )`$ throughout the integration domain. To arrive at a decision we therefore need to solve a discretized version of the pde, scrutinize the computational process, and interpret the numerical evidence so gained, a task we will undertake in part II of our investigation .

## Acknowledgments

The author gratefully acknowledges financial support from Fonds zur Förderung der wissenschaftlichen Forschung (Austrian Science Fund, FWF) under projects P14371-TPH, P15758-N08, and J2380-N08.

## Appendix A Notational conventions and additional relations

A most detailed account of the process of re-writing the pde for the free energy into the quasilinear eq. (1) can be found, alongside the explicit expressions for the pde coefficients $`d_{0i}`$, in appendix A of ref. . Further specializing these to take into account the elimination of the core condition and the purported density-independence of the potential we obtain

$$\begin{array}{cc}\hfill d_{00}=& +\frac{\stackrel{~}{\varphi }}{Q}\left(\frac{\stackrel{~}{\varphi }_0^2}{\stackrel{~}{𝒦}\stackrel{~}{\varphi }^2}\frac{\stackrel{~}{𝒦}\overline{\epsilon }^2\stackrel{~}{\varphi }_0^2}{\epsilon \stackrel{~}{\varphi }^4}\frac{2f}{\stackrel{~}{\varphi }}\right)\hfill \\ & +\frac{\stackrel{~}{𝒦}}{Q}\left(\frac{\overline{\epsilon }^2\stackrel{~}{\varphi }_0^2}{\epsilon \stackrel{~}{\varphi }^3}\frac{\stackrel{~}{\varphi }_0^2}{\stackrel{~}{𝒦}^2\stackrel{~}{\varphi }}\right)\hfill \\ & +\frac{^2(1/\stackrel{~}{𝒦})}{\varrho ^2}\frac{Q^2}{4\pi ^2}\frac{\overline{\epsilon }^2\stackrel{~}{\varphi }_0}{\epsilon \stackrel{~}{\varphi }},\hfill \\ \hfill d_{02}=& \frac{Q^2}{4\pi ^2}\frac{\overline{\epsilon }^2}{\epsilon \stackrel{~}{\varphi }_0}.\hfill \end{array}$$
(A1)
In the above expressions we have suppressed the obvious function arguments $`Q`$ and $`\varrho `$, and we rely on the notational convention introduced in our earlier work on hrt according to which superscripts indicate the system a quantity refers to and a tilde signals Fourier transformation. We also make use of the following auxiliary quantities:

$$\begin{array}{cc}\hfill u_0(r)& =\varphi (r)/\stackrel{~}{\varphi }_0,\hfill \\ \hfill \stackrel{~}{\varphi }_0& =\stackrel{~}{\varphi }(0),\hfill \\ \hfill \stackrel{~}{𝒦}(k,\varrho )& =\frac{1}{\varrho }+\stackrel{~}{c}_2^{\mathrm{ref}}(k,\varrho ),\hfill \\ \hfill \mathrm{ln}\epsilon (Q,\varrho )& =f(Q,\varrho )\stackrel{~}{u}_0^2(Q)\stackrel{~}{\varphi }(Q)/\stackrel{~}{𝒦}(Q,\varrho ),\hfill \\ \hfill \overline{\epsilon }(Q,\varrho )& =\epsilon (Q,\varrho )1.\hfill \end{array}$$
(A2)
As far as the definition of $`𝒦`$ is concerned, the ideal gas term involving $`1/\varrho `$ is customarily included in the definition of the hard sphere reference system direct correlation function $`c_2^{\mathrm{ref}}`$ throughout much of the literature on hrt. Furthermore, when explicitly taking into account the core condition $`𝒦`$ is typically augmented by a truncated series expansion of the deviation of the direct correlation function of the $`Q`$ system inside the core from that of the reference system, a scheme designed to allow one to side-step the need for costly numerical Fourier transformations that would otherwise arise from the inversion of the Ornstein-Zernike equation .

The connection to the thermodynamics as well as to the original formulation of the theory is afforded by the relations linking $`f`$ and the auxiliary quantities defined in eq. (A2) to the first derivative with respect to $`Q`$ of the free energy $`A^{(Q)}`$ of the system at cutoff $`Q`$: From the expressions given in ref.  one may easily show that the evolution of $`A^{(Q)}`$ is given by

$$\begin{array}{cc}\hfill \frac{}{Q}\frac{\beta A^{(Q)}}{V}& =\frac{Q^2}{4\pi ^2}\left(\mathrm{ln}\epsilon \stackrel{~}{\varphi }\varrho \right)\text{for }Q>0,\hfill \\ \hfill \frac{\beta A^{(0)}}{V}& =\underset{Q0^+}{lim}\frac{\beta A^{(Q)}}{V}\frac{\varrho ^2\stackrel{~}{\varphi }_0}{2}.\hfill \end{array}$$
(A3)
The structure of the $`Q`$ system follows from the closure relation that imposes a direct correlation function $`c_2^{(Q)}(r,\varrho )`$ of the form

$$c_2^{(Q)}=c_2^{\mathrm{ref}}+\varphi ^{(Q)}+\gamma _0^{(Q)}u_0$$
(A4)
parameterized by the single scalar $`\gamma _0^{(Q)}(\varrho )`$ that is determined throughout $`𝒟`$ from thermodynamic consistency in the form of the compressibility sum rule

$$\frac{1}{\varrho }+\stackrel{~}{c}_2^{(Q)}(0,\varrho )=\frac{^2}{\varrho ^2}\frac{\beta A^{(Q)}(\varrho )}{V}.$$
(A5)
In the spirit of loga/orpa, when implementing the core condition eq. (A4) is taken to hold only for $`r>\sigma `$, and the direct correlation function inside the core is optimized so as to approximately minimize the pair distribution function there. As shown in appendix A.3 of ref. , the isothermal compressibility $`\kappa _T^{(Q)}`$ of the $`Q`$ system is also readily evaluated and, in the limit $`Q0`$, found to reduce to

$$\kappa _T=\kappa _T^{(0)}=\frac{\beta \overline{\epsilon }}{\varrho ^2\stackrel{~}{\varphi }_0}=\frac{\overline{\epsilon }}{\varrho ^2\stackrel{~}{w}(0)}.$$
(A6)


format: text

text:
# Inelastic H+Li and H-+Li+ collisions and non-LTE Li I line formation in stellar atmospheres

## 1 Introduction

Lithium abundances in stellar atmospheres are key observational parameters in astrophysics. They tell us much about stellar evolution and nucleosynthesis, Big Bang nucleosynthesis, and the role of cosmic rays in Galactic chemical evolution, to name a few (see Lambert lambert (1993); Carlsson et al. carlsson (1994) and references therein). Stellar Li abundances are in practice determined from the 670.8 nm resonance line of Li I and in a few cases the 610.4 nm subordinate line. Thus, the ability to accurately interpret these lines is important. As Li I has a low ionization potential, non-local thermodynamic equilibrium (non-LTE) effects can be expected to be important (Steenbock & Holweger steenbock\_holweger (1984); Carlsson et al. carlsson (1994); Asplund et al. asplund (2003)).

Steenbock and Holweger (steenbock\_holweger (1984)) were the first to suggest that inelastic collisions with hydrogen might be important in thermalizing Li I. In photospheres of late-type stars neutral hydrogen atoms outnumber electrons by typically $`N_H/N_e10^4`$, even more in metal poor stars. Thus, weight of numbers may overcome the expectation that electron collisions are a more efficient thermalizing mechanism firstly due to their higher thermal velocity and therefore higher collision rate, and secondly since for speeds of interest electron cross sections are generally larger than those for neutrals since the neutral collisions are adiabatic while the electron collisions are not (e.g. Lambert lambert (1993)). In such photospheres $`kT0.2`$$`0.6`$ eV and thus for typical optical transitions the low energy collisions just above the threshold are most important in determining the collision rates. Steenbock & Holweger, in the absence of any better experimental or theoretical guidance, calculated the collisional excitation and ionization rates by modifying a formula from Drawin (drawin (1968)) for H+H which is itself a modification of the classical Thomson formula for excitation by electrons (cf. Lambert lambert (1993)). In recent years some work has been devoted to improve these highly uncertain H collision cross-sections. Low-energy experimental data have been obtained (Fleck et al. fleck (1991)) and ab initio calculations performed down to the threshold (Belyaev et al. belyaev99 (1999)) for the Na+H system with good agreement between theory and experiment. However, these results are for only the very lowest states of this one system. For non-LTE calculations, estimates are required for transitions between all states which might affect the population of the states of interest. It is for this reason that the modified Drawin formula (often with a scaling factor determined astrophysically) is still often in use despite that it is known to overestimate the Na(3s)+H$``$Na(3p)+H collision rate by several orders of magnitude.

Recently, calculations for Li+H collisions for all transitions between the states Li($`2s`$,$`2p`$,$`3s`$,$`3p`$,$`3d`$,$`4s`$,$`4p`$,$`4d`$,$`4f`$)+H($`1s`$) and Li++H- have been performed (Belyaev & Barklem paper1 (2003)). This paper should be consulted for details and justification of the approaches used. In brief summary, on the basis of the interpretation of the experimental and theoretical results for Na+H (Fleck et al. fleck (1991); Belyaev et al. belyaev99 (1999)), non-adiabatic regions associated with avoided ionic crossings in the adiabatic potential curves are expected to provide the basic mechanism for the transitions at low energy. A simple physical interpretation of this process is that at certain separations during the collision the optical electron associated with the Li atom may tunnel to the H atom, leading to a predominantly ionic charge distribution (Li+H-) of the LiH quasi-molecule. Later, at another avoided crossing, the electron may tunnel back towards the Li atom into a different covalent charge distribution leading asymptotically to a different final Li state. The electron may also remain with the H atom leading to an ion-pair production reaction (Li+H$``$Li++H-). For the transitions between the lowest four states a fully quantum mechanical approach has been used, while for the transitions involving higher states, a multichannel Landau-Zener model has been used. The model approach is justified for the higher crossings where the requirements of the Landau-Zener model are better fulfilled.

For the charge transfer reactions Li+H$``$Li++H-, referred to as ion-pair production and mutual neutralization, our calculations use the model approach. More elaborate calculations, using a diabatic molecular expansion method, have been made by Croft et al. (1999b ). This process has a remarkably large cross section near the threshold in alkali-hydrogen systems, particularly to and from the first excited $`s`$ state. This has been predicted theoretically (Bates & Boyd bates\_boyd (1956); Janev & Radulović janev (1978); Croft et al. 1999a ; Dickinson et al. dickinson (1999)) and is consistent with experiment (Peart & Hayton peart\_hayton (1994)). Despite that this process has been known to be remarkably efficient in various fields (e.g. plasma physics, early universe chemistry) for some time, to our knowledge this is the first time it has been considered in a stellar atmosphere application of Li.

## 2 Collision rate coefficients

The rate coefficients $`\sigma v`$ for excitation and de-excitation H($`1s`$)+Li($`nl`$)$``$H($`1s`$)+Li($`n^{}l^{}`$), and the charge transfer reactions Li$`(nl)`$+H$``$Li++H- are presented in Table 1. The data have been obtained by integrating a Maxwellian velocity distribution with cross-sections from Belyaev & Barklem (paper1 (2003)) in the case of excitation and de-excitation, and from Croft et al. (1999b ) in the case of mutual neutralization, with ion-pair production cross sections obtained via the principle of detailed balance. Thus the data represent, in our opinion, the best available estimates. Note that all forward and backward rates were computed separately, but should be connected by the detailed balance relation. Inspection of the data reveals the remarkably large rate coefficients for Li$`(nl)`$+H$``$Li++H-, particularly for the $`n=3`$ states.

For the temperatures here the uncertainty in the rate coefficients is determined by the uncertainty in the cross sections near the threshold. The uncertainty in the input cross sections is discussed in detail in the relevant papers. In summary, the accuracy for transitions between the lower states ($`2s`$,$`2p`$ and $`3s`$) is estimated to be better than a factor of 2. For the excitations involving higher states the data are expected to be accurate within an order of magnitude. Comparison with experimental results for the total neutralization cross-section with D- (which uses the same interaction potentials as for H-) finds the theoretical result about 20% greater at the lowest measured energy 0.68 eV (see Croft et al. 1999a ), with a general trend of the disagreement becoming larger at lower energies, recalling that the collision energies of importance here are $`E0.20.6`$ eV.

Lambert (lambert (1993)) presented unpublished estimated rate coefficients for the $`2s`$$``$$`2p`$ transition computed by Allen & Dickinson. Our data are roughly an order of magnitude smaller.

Comparisons were made with the modified Drawin formula for all optically allowed transitions at 5000 K. Except for an extreme outlier ($`3s`$$``$$`4p`$), the Drawin formula gives a rate coefficient between roughly one and six orders of magnitude greater than our result. This reinforces the view that the Drawin formula typically greatly over-estimates the collisional efficiency, as judged by the available (admittedly rather scant) evidence.

## 3 Application to non-LTE formation of Li lines

To assess the impact of the data on the interpretation of stellar spectra, we have performed non-LTE calculations for Li in cases of particular interest using the MULTI-code for 1D statistical equilibrium problems (Carlsson multi (1986)) and the MULTI3D-code (Botnen botnen (1997); Botnen & Carlsson botnen\_carlsson (1999); Asplund et al. asplund (2003)) for the 3D case. The adopted model atom is the 21-level atom compiled by Carlsson et al. (carlsson (1994)). To test the impact of the new data three model atoms were used, namely the original model atom with no H collisions (noH), the model with the new H collision data added (wH), and a model with the H collision rates halved (w0.5H); the LTE results have been obtained with the same codes using an atom with extremely large collisional cross-sections to ensure consistency with the non-LTE calculations. The 1D model atmospheres are from the MARCS code (Gustafsson et al. marcs (1975) and subsequent updates). The 3D models and calculations are as detailed in Asplund et al. (asplund99 (1999), asplund (2003)). Calculations were performed for the Sun, the metal-poor subgiant HD140283 and the metal-poor turnoff star HD84937. Equivalent widths for the 670.8 nm line are presented in Table 2. Departure coefficients for the lower states in the 1D non-LTE models for the Sun and HD140283 are presented in Fig. 1.

Although quantitatively different, the overall effects of the H collisions are similar in the 1D and 3D cases. Almost all of the difference resulting from the inclusion of the H collision data is caused by the reaction Li$`(3s)`$+H$``$Li++H-. Removing all H collisions but this reaction results in equivalent widths within 1% of the results for the wH model. This reaction causes a stronger collisional coupling between the ionization stages, nearly thermalizing the $`3s`$ state. Photon suction (Bruls et al. bruls (1992)) results in the other levels being affected, in particular the ground state is pushed towards LTE, or even into overpopulation as seen in the case of HD140283 in 1D, at around $`\mathrm{log}\tau _{500}=2`$. As clear from Table 2, the results are not particularly sensitive to the details of the H collisional cross-sections with only small differences resulting when decreasing all H cross-sections by 50% (w0.5H). The remaining uncertainties in the charge transfer reactions are thus expected to have limited impact on the Li I line formation.

In terms of abundances, the inclusion of charge transfer reactions lowers the derived values by about 0.05–0.10 dex for all three stellar models, whether in 1D or 3D. While the 1D non-LTE abundance corrections (1D non-LTE – 1D LTE) are small ($`<`$0.05 dex) for the relevant stars (Table 2, see also Carlsson et al. carlsson (1994)), the 3D non-LTE effects (3D non-LTE – 3D LTE) are substantial for metal-poor stars (Asplund et al. asplund (2003)). However, since the differences between 3D LTE and 1D LTE are also significant and of opposite sign (Asplund et al. asplund99 (1999)), the net difference between 3D non-LTE and 1D non-LTE is less dramatic. For the two metal-poor stars this difference amounts to about 0.05–0.10 dex lower abundance than in 1D111Note that it is not appropriate to directly compare the equivalent widths in Table 2 to estimate the net abundance difference between 1D and 3D as the 3D results are here only shown for one snapshot. The proper procedure is to add the 3D non-LTE corrections obtained here to the 3D LTE corrections published by Asplund et al. (asplund99 (1999)) based on a long time-sequence with better numerical resolution.. This has a direct bearing on the question of the primordial Li abundance. The observational evidence suggests a plateau at $`\mathrm{log}ϵ_{\mathrm{Li}}=2.12.3`$ with 1D non-LTE abundance corrections taken into account (Ryan et al. ryan (2000); Bonifacio bonifacio (2002)), with the abundance range partly reflecting different effective temperature scales for halo stars. By slightly decreasing these abundances, the inclusion of 3D models and charge transfer reactions slightly aggravates the discrepancy between the observed Li abundances and the predicted primordial Li abundance from Big Bang nucleosynthesis. High-redshift deuterium measurements and analyses of the cosmic microwave background yield a relatively high baryon density of the Universe (e.g. Burles et al. burles (2001)). In particular, the recent WMAP analysis (Bennett et al. wmap (2003)) corresponds to a primordial Li abundance of $`\mathrm{log}ϵ_{\mathrm{Li},\mathrm{p}}=2.69\pm 0.13`$ based on standard Big Bang nucleosynthesis (Burles et al. burles (2001)). The small amount of scatter in observed halo Li abundances is very difficult to reconcile with uniform stellar Li depletions of 0.4 dex or more. Likewise, it would require improbably high effective temperatures to resolve the disagreement. Clearly, halo stars and their spectrum formation are not as well understood as perhaps thought.

## 4 Discussion

We have found that while excitation and de-excitation by inelastic H collisions are not important in thermalizing Li, the charge exchange reactions involving the Li($`3s`$) state, are quite important. This mechanism is almost certainly important also in other alkalis (Na, K, Rb etc), where large cross-sections are also found (Dickinson et al. dickinson (1999); Janev & Radulović janev (1978)). For other elements, further detailed calculations will be required to determine whether inelastic H collisions, including excitation and charge exchange, are important. The main obstacle for such calculations at present is the quantum chemical potentials and couplings required for a large number of molecular states. Though systems of avoided ionic crossings exist in all metal hydride quasi-molecules, we note other mechanisms may be important. As suggested by Lambert (lambert (1993)), the importance of details in the molecular dynamics, such as avoided ionic crossings, in this process makes a simple general recipe impossible, in contrast to collisional broadening and depolarization by H where such effects play a lesser role (Anstee & O’Mara anstee (1995); Barklem & O’Mara barklem (2001); Derouich et al. derouich (2003)). Ab initio potentials for alkaline-earth-hydrides are possible with current techniques (Chambaud & Lévy chambaud\_levy (1989)). However, some cases of particular astrophysical interest such as O+H and Fe+H present a considerable challenge.

###### Acknowledgements.

We are indebted to F.X. Gadéa and A.S. Dickinson for providing us with the quantum chemical data used in computing the cross sections, and their results for mutual neutralization. We thank D. Kiselman for encouraging us to investigate inelastic H collisions. We gratefully acknowledge support from the Swedish Research Council, the Swedish Royal Academy of Sciences, the Göran Gustafsson Foundation and the Australian Research Council.


text:
# Photon generation from vacuum in nondegenerate cavities with regular and random periodic displacements of boundaries

## 1 Introduction

Among many interesting classical and quantum phenomena in cavities with moving boundaries, the effect of photon creation from the vacuum seems to be the most impressive. The first rough (although not always quite correct) estimations were made by several authors in Refs. . More detailed and consistent studies performed for the past decade in the frameworks of different approaches  (for other references see review ) predicted, indeed, an exponential growth of the field energy in the case of periodic motion and under the resonance conditions, when the walls vibrate at the frequency which is a multiple of the unperturbed field eigenfrequency.

Now, an experimental verification of theoretical predictions becomes a difficult but realizable task. An idea which seems to be the most realistic at the moment is to realize an effective moving mirror by means of periodic creation of a conducting layer from a semiconductor film posed on the surface of the cavity wall and irradiated by powerful laser impulses . Different prototypes of this idea were discussed in , but not in a periodic regime (for other schemes see, e.g., ; however, their realizations seem to be more difficult). The effective mirror scheme has several advantages over the primary idea of an acoustic excitation of harmonical surface vibrations of the wall .

The main of them is the removal of severe limitations on the maximal dimensionless amplitude of the displacement of the surface, which resulted from tremendous internal stresses inside the wall whose surface vibrates at the necessary frequency having an order of several GHz (or higher). In the effective mirror scheme, the maximal displacement of the boundary depends of the thickness of the semiconductor film and the laser power, and it can be made several orders of magnitude larger than the maximal displacements (not exceeding $`10^8`$cm ) which could be achived in the scheme based on the acoustic excitation of the surface vibrations. As a consequence, one can diminish (in the same proportion) the number of oscillations of the boundary necessary to produce photons in the cavity, thus relaxing the requirements for the $`Q`$-factor of the cavity and the admissible detuning from the exact resonance . Another advantage is a possibility to use periodic pulses of arbitrary shape, whose period can be much longer than the period of oscillations of the chosen mode of electromagnetic field.

However, a practical realization of the scheme depends on its robustness against inevitable experimental imperfections, such as, e.g., nonexact periodicity of laser pulses (creating free carriers in the semiconductor layer) or fluctuations in their power. Therefore, the aim of our article is to give some estimations of possible deviations from ideal results caused by such imperfections.

## 2 Basic relations

The most part of theoretical studies on the nonstationary (dynamical) Casimir effect was performed within the framework of the model of a one-dimensional (Fabry–Perot) cavity with equidistant spectrum of unperturbed eigenfrequencies of the field, when all modes are coupled in the resonance case. However, in realistic three-dimensional cavities the spectrum is, as a rule, nonequidistant. This results in great simplifications in the resonant case, when periodic motion of the boundary excites only one resonant mode, whereas the response of all other modes can be neglected in the long-time limit . Therefore, we confine ourselves here to this simplest case (interesting phenomena which could occur in cavities with few accidentally resonant modes, when the ratio of their unperturbed eigenfrequencies is integral number due to some additional symmetry, were considered recently in ).

In this case, as was shown in , the problem of photon generation in the selected mode is reduced to the problem of excitation of a quantum oscillator with a time dependent frequency $`\omega (t)`$, which is determined by the instantaneous geometry of the cavity. For example, in a rectangular cavity with fixed dimensions $`L_x,L_y`$ and a variable dimension $`L_z(t)`$ we have

$$\omega (t)=\pi c\left[\frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2}+\frac{n_z^2}{L_z^2(t)}\right]^{1/2}.$$
(1)
For fifty years passed after the seminal paper by Husimi , this last problem was studied in detail in numerous publications (for reviews see, e.g., ). It appears that all properties of the quantum oscillator are determined, as a matter of fact, by the fundamental set of solutions of the classical equation of motion

$$\ddot{x}+\omega ^2(t)x=0.$$
(2)
Here we actually need only one consequence of the general theory, namely, that the mean number of quanta generated from the initial oscillator ground state due to the time dependence of frequency is given by the energy reflection coefficient from an effective “potential barrier” given by the function $`\omega ^2(t)`$. More precisely, if $`\omega (t)=\omega _i`$ for $`t\mathrm{}`$ and $`\omega (t)=\omega _f`$ for $`t\mathrm{}`$ (this means that the cavity wall is supposed to assume some fixed positions before and after the experiment), then one has to calculate the coefficients of the asymptotical form of the solutions to Eq. (2),

$$x(t\mathrm{})=\omega _f^{1/2}\left[\xi e^{i\omega _ft}+\eta e^{i\omega _ft}\right],$$
(3)
satisfying the initial condition $`x(t)=\omega _i^{1/2}e^{i\omega _it}`$ for $`t\mathrm{}`$. Due to the unitarity of evolution, these coefficients obey the identity

$$|\xi |^2|\eta |^2=1.$$
(4)
The mean number of quanta for $`t\mathrm{}`$ can be expressed in terms of $`\xi `$ and $`\eta `$ as follows ,

$$𝒩=|\eta |^2=R/T,$$
(5)
where $`R|\eta /\xi |^2`$ and $`T1R|\xi |^2`$ can be interpreted as energy reflection and transmission coefficients from the effective potential barrier.

Eq. (5) shows that no photons can be created during a single pass of the wall from one position to another, because for monotonous functions $`\omega (t)`$ the reflection coefficient is limited by the Fresnel formula

$$R_1=\left(\frac{\omega _i\omega _f}{\omega _i+\omega _f}\right)^2$$
(6)
corresponding to very rapid change of position (for the time much less than the period of field oscillations) . Since real boundaries in laboratory experiments can only move with velocities much less than the speed of light, the process is almost adiabatic, and $`R_11`$. However, it is well known that the reflection coefficient can be made very close to unity in the case of periodic variations of parameters due to interference effects.

For a periodic motion of the wall, the function $`\omega ^2(t)`$ is also periodic. We assume that it has a form like that shown in Fig. 1. The most important assumption is that each “effective potential barrier” is well separated from the next one by some interval of time where $`\omega =const`$ (i.e., that the wall moves for some period of time, returns to its initial position, stays at this position for some time, and then repeats the cycle). Thus we exclude monochromatic oscillations of the boundary, which have been already studied in detail in .

In such a case, each “barrier” can be completely characterized by two complex amplitude reflection coefficients and two complex amplitude transmission coefficients, which connect the “plane waves” coming from the “left” and from the “right”. Namely, if some “barrier” begins at $`t=0`$ at terminates at $`t=t_{}`$, then one can write two independent solutions of Eq. (2) outside the “barrier” as (hereafter we consider for simplicity the case of equal initial and final constant values of the frequency $`\omega `$)

$`x^{()}(t)=\{\begin{array}{cc}e^{i\omega t}+r^{()}e^{i\omega t},\hfill & t<0\hfill \\ t^{()}e^{i\omega t},\hfill & t>t_{}\hfill \end{array}`$ (9)
$`x^{(+)}(t)=\{\begin{array}{cc}t^{(+)}e^{i\omega t},\hfill & t<0\hfill \\ e^{i\omega t}+r^{(+)}e^{i\omega t},\hfill & t>t_{}\hfill \end{array}`$ (12)
(We use the letter $`t`$ without supescripts to denote the time variable, whereas the same letter supplied with supescripts means the amplitude transmission coefficient; we hope this will not lead to a confusion.) The coefficients $`r^{(\pm )}`$ and $`t^{(\pm )}`$ are not independent, because equation (2) is invariant with respect to complex conjugation (the frequency $`\omega (t)`$ is real). The following relations hold (see, e.g., ; the simplest way to obtain these relations is to calculate Wronskians for suitable pairs of independent solutions):

$`t^{()}=t^{(+)}t^{(\pm )},`$ (13)
$`r^{()}t^{(\pm )}+r^{(+)}t^{(\pm )}=0,`$ (14)
$`|r^{()}|^2+|t^{()}|^2=|r^{(+)}|^2+|t^{(+)}|^2=1.`$ (15)
Suppose now that we have several well separated “barriers”, the $`n`$th “barrier” (possessing coefficients $`r_n^{(\pm )}`$ and $`t_n^{(\pm )}`$ if it begins at $`t=0`$) being shifted in time by the value $`t_{n1}`$ (with $`t_0=0`$) with respect to the beginning of the first one. Then one can combine the first $`n`$ barriers in a single effective barrier, characterized by the amplitude reflection and transmission coefficients $`\rho _n^{(\pm )}`$ and $`\tau _n^{(\pm )}`$. There are many different methods enabling to express these coefficients in terms of $`r_k^{(\pm )}`$ and $`t_k^{(\pm )}`$, $`k=1,\mathrm{},n`$: see, e.g.,  and references therein. We use the method of transfer matrix  (see also  for historical review and  for the most recent applications).

Let us write the solutions to Eq. (2) after the $`k`$th barrier (in the region where $`\omega (t)=\omega =const`$) as

$$x^{(k)}(t)=a_ke^{i\omega t}+b_ke^{i\omega t},$$
(16)
with suffix $`0`$ related to the solution before the first barrier. Evidently, every two sets of the nearest constant coefficients, $`(a_{k1},b_{k1})`$ and $`(a_k,b_k)`$, are related by means of a linear transformation

$$\left(\begin{array}{c}a_k\\ b_k\end{array}\right)=M_k\left(\begin{array}{c}a_{k1}\\ b_{k1}\end{array}\right).$$
(17)
Comparing Eqs. (16) and (17) with (9) and (12), and taking into account the identities (13)-(15), one can express the elements of matrix $`M_k`$ as follows:

$$M_k=\begin{array}{cc}f_k^{}& g_k^{}\\ g_k& f_k\end{array},$$
(18)
where

$$f_k[t_k^{(\pm )}]^1,g_kr_k^{()}/t_k^{(\pm )}.$$
(19)
Consequently, each matrix $`M_k`$ is unimodular:

$$detM_k=|f_k|^2|g_k|^21.$$
(20)
We shall use the notation

$$_n=\begin{array}{cc}F_n^{}& G_n^{}\\ G_n& F_n\end{array}$$
(21)
for the total transfer matrix of $`n`$ barriers, shifted in time by $`t_k`$, $`k=1,\mathrm{},n1`$, with respect to the initial instant $`t=0`$. Applying consecutively the transformations (17) with account of the phase shifts $`\theta _k\omega t_k`$ (with $`\theta _0=0`$), one can obtain the following matrix formula:

$$_n=\mathrm{\Phi }_{n1}M_n\mathrm{\Phi }_{n1}^{}\mathrm{\Phi }_{n2}M_{n1}\mathrm{}\mathrm{\Phi }_1M_2\mathrm{\Phi }_1^{}M_1,$$
(22)
where

$$\mathrm{\Phi }_k\begin{array}{cc}\mathrm{exp}(i\theta _k)& 0\\ 0& \mathrm{exp}(i\theta _k)\end{array}.$$
(23)
Eq. (22) is equivalent to the recursive relations

$$_k=\mathrm{\Phi }_{k1}M_k\mathrm{\Phi }_{k1}^{}_{k1},$$
(24)
$$F_k=g_kG_{k1}^{}e^{2i\theta _{k1}}+f_kF_{k1},$$
(25)
$$G_k=g_kF_{k1}^{}e^{2i\theta _{k1}}+f_kG_{k1},$$
(26)
where $`F_0=1`$ and $`G_0=0`$. One can verify that Eqs. (25) and (26) preserve the identities

$$|F_k|^2|G_k|^21.$$
(27)
The total mean number of quanta (5) created after $`n`$ impulses is nothing but

$$𝒩_n|G_n|^2.$$
(28)
The coefficients $`\rho _n`$ and $`\tau _n`$ themselves obey the nonlinear recurrence relations

$$\rho _n=\frac{\rho _{n1}+r_ns_{n1}\mathrm{exp}(2i\theta _{n1})}{1+\rho _{n1}^{}r_ns_{n1}\mathrm{exp}(2i\theta _{n1})},$$
(29)
$$\tau _n=\frac{\tau _{n1}t_n^{(\pm )}}{1+\rho _{n1}^{}r_ns_{n1}\mathrm{exp}(2i\theta _{n1})},$$
(30)
where $`s_n\tau _n/\tau _n^{}`$, $`\rho _1=r_1^{()}`$, $`\tau _1=t_1^{(\pm )}`$.

## 3 Strictly periodic motion of boundaries

Formula (22) becomes especially useful in the case of strictly periodic motion of boundaries, when all one-barrier matrices $`M_k`$ coincide with $`M_1`$, and $`\mathrm{\Phi }_k\mathrm{\Phi }_1^k`$:

$$_n=\mathrm{\Phi }_1^n(\mathrm{\Phi }_1^{}M_1)^n.$$
(31)
Obviously, the matrix $`\mathrm{\Phi }_1^{}M_1`$ is unimodular, as well as matrix $`M_1`$. Therefore, we can use the well-known formula for the powers of any two-dimensional unimodular matrix $`S`$ (see, e.g., ):

$$S^n=U_{n1}(z)SU_{n2}(z)E,z\frac{1}{2}\text{Tr}S,$$
(32)
where $`E`$ means the unit matrix and $`U_n(z)`$ is the Tchebyshev polynomial of the second kind.

In the case discussed, introducing the new parameter $`\nu `$ according to the relations (hereafter $`\theta \theta _1`$, $`gg_1`$, $`ff_1`$)

$$z=\text{Re}\left(fe^{i\theta }\right)=\pm \mathrm{cosh}\nu ,$$
(33)
and using the hyperbolic representation of the Tchebyshev polynomial

$$U_n(\mathrm{cosh}\nu )\frac{\mathrm{sinh}[(n+1)\nu ]}{\mathrm{sinh}(\nu )},$$
(34)
we arrive at the following explicit expressions for the elements of matrix $`_n`$ (for simplicity, we suppose that the sign of variable $`z`$ in Eq. (33) is positive; this sign does not affect the final result (37) for the number of created photons):

$$G_n=g\frac{\mathrm{sinh}(n\nu )}{\mathrm{sinh}(\nu )}e^{i\theta (n1)},$$
(35)
$$F_n=f\frac{\mathrm{sinh}(n\nu )}{\mathrm{sinh}(\nu )}e^{i\theta (n1)}\frac{\mathrm{sinh}[(n1)\nu ]}{\mathrm{sinh}(\nu )}e^{i\theta n}.$$
(36)
Consequently, the number of created photons after $`n`$ impulses equals

$$𝒩_n=|g|^2\frac{\mathrm{sinh}^2(n\nu )}{\mathrm{sinh}^2(\nu )}.$$
(37)
The generation of photons is possible provided parameter $`\nu `$ is real, i.e.,

$$|f||\mathrm{cos}(\phi \theta )|>1,f|f|e^\phi .$$
(38)
Since $`|f|>1`$, one has to ajust the phase shift $`\theta `$ to the phase of the inverse transmission coefficient $`\phi `$. The maximal effect is achieved for

$$\theta =\theta _{res}\phi +\pi m$$
(39)
(with $`m`$ an integer). This is equivalent to the following relation between the periodicity of impulses $`T=\theta /\omega `$, the period of the electromagnetic field oscillations in the mode concerned $`T_0=2\pi /\omega `$ and the phase $`\phi `$:

$$T=\frac{T_0}{2}\left(m+\phi /\pi \right).$$
(40)
In particular, for nonmonochromatic oscillations of the boundary, the field mode can be excited not only under the condition of the parametric resonance $`T=T_0/2`$, but the period of motion of the boundary can be greater than the period of field oscillations. This fact may be important for possible experimental arrangements.

Under the condition (40), $`\mathrm{cosh}(\nu )=|f|`$ and $`\mathrm{sinh}(\nu )=|g|`$, so that

$`𝒩_n^{(res)}`$ $`=`$ $`\mathrm{sinh}^2\left(n\nu \right)`$
$``$ $`{\displaystyle \frac{1}{4}}\left[\left({\displaystyle \frac{1+|r|}{1|r|}}\right)^{n/2}\left({\displaystyle \frac{1|r|}{1+|r|}}\right)^{n/2}\right]^2,`$
where $`rr_1^{()}`$ is the amplitude reflection coefficient from each barrier (remember the relations (19) between different coefficients). If $`|r|1`$, then one can replace $`\nu `$ by $`|r|`$. Then for $`n|r|>1`$ we have

$$𝒩_n^{(res)}\frac{1}{4}\mathrm{exp}(2n|r|).$$
(42)
In the case of a harmonically oscillating boundary near parametric resonance,

$$\omega (t)=\omega _0\left[1+2\epsilon \mathrm{sin}(\mathrm{\Omega }t)\right],\mathrm{\Omega }=2(\omega _0+\delta )$$
(43)
($`\omega _0`$ is the unperturbed field eigenfrequency, and $`\mathrm{\Omega }`$ is the frequency of the wall vibrations), the number of quanta created from vacuum is given by the expression  (under the conditions $`|\delta |\omega _0`$ and $`|\epsilon |1`$)

$$𝒩=\frac{\mathrm{sinh}^2\left(\omega _0\epsilon \gamma t\right)}{\gamma ^2},\gamma ^2=1\frac{\delta ^2}{(\omega _0\epsilon )^2}.$$
(44)
For $`\delta =0`$ (exact resonance ), (44) coincides with (3), if one identifies $`\nu _{ef}=\pi \epsilon `$ (the number of full cycles equals $`n=\mathrm{\Omega }t/(2\pi )=\omega _0t/\pi `$).

We see that the concrete form of the law of motion of the boundary $`L(t)`$ turns out to be not very important for the exponential dependence of the number of created quanta on the number of wall’s oscillations $`n`$ for $`n\nu 1`$ (although it influences the concrete value of the coefficient $`\nu `$). What is important, it is the fulfillment of the condition (40).

The generation of quanta in the harmonic case becomes impossible if $`\delta >\omega _0\epsilon `$. This means that the phase deviation from the resonance value accumulated for one cycle, $`\mathrm{\Delta }\theta =\delta T=\delta \pi /\omega _0`$, should not exceed the critical value $`\mathrm{\Delta }\theta _c=\pi \epsilon =\nu _{ef}`$. The same result holds, in fact, in the general nonmonochromatic case. Indeed, for small values of $`|g|1`$ we have $`|f|=\sqrt{1+|g|^2}1+|g|^2/2`$. Then Eq. (33) results in $`\nu ^2|g|^2|\mathrm{\Delta }\theta |^2`$, where $`\mathrm{\Delta }\theta \theta \theta _{res}`$ (see Eq. (39) for the definition of the resonance phase $`\theta _{res}`$), so that

$$𝒩_n\frac{|g|^2}{|g|^2|\mathrm{\Delta }\theta |^2}\mathrm{sinh}^2\left(n\sqrt{|g|^2|\mathrm{\Delta }\theta |^2}\right).$$
(45)
The critical regular phase shift per cycle equals $`|\mathrm{\Delta }\theta _c|=|g|=\mathrm{sinh}(\nu _{res})\nu _{res}`$. By the order of magnitude, taking into account the Fresnel limit (6), we can estimate the maximal phase deviation per cycle (from the resonance value) as $`|\mathrm{\Delta }\theta _c|<\mathrm{\Delta }L/(2L)`$, where $`\mathrm{\Delta }L`$ is the change of the cavity length.

## 4 Influence of irregularities of periodic motion

In real experiments, it can be difficult to maintain the resonance conditions $`g_kg_1`$, $`f_kf_1`$, $`\theta _k=k\theta _{res}`$ exactly. Therefore, it is important to evaluate the influence of systematic and random deviations (fluctuations) from these conditions.

### 4.1 Amplitude variations

It seems that fluctuations of absolute values of the effective reflection and transmission coefficients, or coefficients $`|g_k|`$ and $`|f_k|`$, are less important than fluctuations of phases $`\theta _k`$ (i.e., unequal intervals between pulses). To show this, let us notice that for identical barriers, the resonance condition (39) is equivalent to the requirement that the second term in the denominators of fractions in Eqs. (29) and (30) is real. In such a case, all coefficients $`\rho _n`$ have identical phases. Moreover, one can adjust the moments of time $`t_k`$ in such a way that this property holds even for barriers of different heights. Under such “generalized resonance conditions”,

$$\zeta _{n+1}\zeta _n2\phi _n+2(\theta _{n+1}\theta _n)=0\text{mod}(2\pi ),$$
(46)
where $`r_n=|r_n|\mathrm{exp}(i\zeta _n)`$ and $`f_n=|f_n|\mathrm{exp}(i\phi _n)`$, Eq. (29) assumes the form

$$|\rho _n|=\frac{|\rho _{n1}|+|r_n|}{1+|\rho _{n1}||r_n|}.$$
(47)
Writing $`|r_n|=\mathrm{tanh}(\nu _n)`$, one can verify that

$$|\rho _n|=\mathrm{tanh}\left(\underset{k=1}{\overset{n}{}}\nu _k\right),$$
(48)
and the number of quanta in the selected mode monotonously increases with time as

$$𝒩=\mathrm{sinh}^2\left(\underset{k=1}{\overset{n}{}}\nu _k\right).$$
(49)
For strictly periodic processes ($`\nu _n=\nu =const`$) formula (49) coincides with (3).

### 4.2 Phase fluctuations

Now we suppose that reflection and transmission coefficients are exactly the same for all barriers, but the strict periodicity is broken:

$$\theta _k=k\theta +\chi _k,f_kf,g_kg,$$
(50)
where $`\chi _k`$ is a stochastic variable. To have a qualitative understanding, what can happen in this situation, let us consider first the case, where $`\chi _k`$ is in fact regular function, namely,

$$\chi _{2k+1}=+\chi ,\chi _{2k}=\chi ,$$

which nonetheless destroys a strict periodicity of pulses. This problem still can be solved analytically, because each two consecutive barriers can be combined in new effective barriers, forming strictly periodic sequence. The total transfer matrix after $`2n`$ pulses can be represented in the following form:

$$_{2n}=\mathrm{\Phi }_1^{2n1}X\left(M^{}X^2M^{}X^2\right)^{n1}M^{}X^1M^{}\mathrm{\Phi }_1$$

where

$$M^{}=M_1\mathrm{\Phi }_1^1=\begin{array}{cc}\stackrel{~}{f}^{}& \stackrel{~}{g}^{}\\ \stackrel{~}{g}& \stackrel{~}{f}\end{array},X=\begin{array}{cc}e^{i\chi }& 0\\ 0& e^{i\chi }\end{array}$$

and $`\stackrel{~}{f}=fe^{i\chi }`$, $`\stackrel{~}{g}=ge^{i\chi }`$. Using again formula (32), we can express elements of the total transfer matrix in terms of the Tchebyshev polynomials of the argument $`z=\text{Re}(\stackrel{~}{f}^2)+|g|^2\mathrm{cos}(4\chi )`$. We omit intermediate calculations and bring here the final simplified result, which is valid for $`|g|1`$, $`n1`$ and $`|\mathrm{\Delta }\theta |1`$, Under these restrictions (which correspond to expected experimental conditions), we can neglect the difference between $`U_n(z)`$ and $`U_{n1}(z)`$ for $`z`$ close to $`1`$ and replace $`\mathrm{cos}(2\mathrm{\Delta }\theta )`$ by $`12(\mathrm{\Delta }\theta )^2`$. Moreover, we can neglect the diference between the values $`𝒩_{2n}`$ and $`𝒩_{2n1}`$. Finally, we arrive at the following simple formula:

$$𝒩_n\frac{\mathrm{cos}^2(2\chi )\mathrm{sinh}^2\left(|g|n\sqrt{\mathrm{cos}^2(2\chi )(\mathrm{\Delta }\theta /|g|)^2}\right)}{\mathrm{cos}^2(2\chi )(\mathrm{\Delta }\theta /|g|)^2}.$$
(51)
We see that the presence of phase “jumps” diminishes the rate of increase of the number of photons, although the growth remains exponential for sufficiently large number of pulses. The admissible level of jumps is determined by the inequality $`|\mathrm{cos}(2\chi )|>|\mathrm{\Delta }\theta |/|g|`$, which shows that “random” deviations from the resonance are much less dangerous than systematic ones. Indeed, the maximal admissible value of the systematic phase shift per period equals $`\mathrm{\Delta }\theta _{max}=|g|`$, in complete acordance with evaluations made at the end of section 3. On the other hand, for zero systematic shift from the exact resonance ($`\mathrm{\Delta }\theta =0`$), exponential growth of the number of photons can be observed for any amplitude of “phase jumps” $`\chi `$, although the growth becomes slower with increase of $`\chi `$ in the interval $`(0,\pi /4)`$:

$$𝒩_n\mathrm{sinh}^2\left(|g|n|\mathrm{cos}(2\chi )|\right).$$
(52)
We have performed numerical experiments on the basis of recursive relations (25) and (26), choosing parameters $`g`$ and $`f=\sqrt{1+g^2}`$ in Eq. (50) real (so that the resonance condition is $`\theta _{res}=\pi m`$; obviously, this choice does not affect the qualitative picture), namely, we chose $`g=0.01`$ (in order to reduce the time of calculations and make the results more vizualizable). The coefficients $`\chi _k`$ were chosen with the aid of a program generating random numbers in the interval $`(\chi ,\chi )`$, for different maximal values $`\chi `$ and systematic displacements $`\theta =\mathrm{\Delta }\theta `$ from $`\theta _{res}=0`$ (which is equivalent to $`2\pi m`$, of course). The results are given in the figures 2-5.

Fig. 2 shows the mean number of photons in the cavity after a fixed number of $`n=500`$ pulses (which generate approximately $`\mathrm{sinh}^2(5)5500`$ photons for $`g=0.01`$ in the case of strict resonance) versus the maximal amplitude of random phase fluctuations $`\chi `$ for fixed values of systematic phase deviations $`\theta `$. The averaging was performed in all cases over $`700`$ realizations consisting of $`500`$ random numbers $`\chi _k`$. The average values $`\chi _k=n^1_{k=1}^n\chi _k`$ varied from $`10^7`$ to $`10^5`$ for each realization (so that no systematic shift was added due to random choice of $`\chi _k`$). Fig. 3 gives the mean number of photons versus $`\theta `$ for different fixed values of $`\chi `$. Fig. 4 shows the relative (and in the insertion – absolute) mean-squared deviations from the average numbers, as functions of $`\chi `$ for $`\theta =0`$. In Fig. 5, the upper curve shows the dependence of the critical amplitude of random phase fluctuations (defined by the condition that the average number of created photons after $`n=500`$ pulses is between $`9`$ and $`10`$) on the regular phase shift $`\theta `$; the lower curve gives the maximal amplitude of phase “jumps”, when the photon generation stops according to Eq. (51). It is clearly seen, that even very high level of random phase fluctuations does not stop the generation of photons, contrary to systematic phase shifts. Moreover, the critical amplitude of truely random fluctuations remains nonzero even for systematic phase shifts exceeding the value $`|g|`$.

### 4.3 Analogy with resistance of disordered chains

It is worth noting that the problem of transmission through a sequence of one-dimensional barriers with randomly varying parameters (in particular, interbarrier distances) has been extensively studied in the context of solid state physics . There, one of objectives is the resistance $``$ of some chain of scatteres, expressed (in dimensionless units) by exactly the same formula (5) which gives the mean photon number in our case. In particular, Landauer  used a consequence of Eq. (29),

$$\frac{|\rho _n|^2}{1|\rho _n|^2}=\frac{|\rho _{n1}|^2+|r_n|^2+2|r_n\rho _{n1}|\mathrm{cos}(\mathrm{\Psi }_n)}{\left(1|\rho _{n1}|^2\right)\left(1|r_n|^2\right)},$$
(53)
$$\mathrm{\Psi }_n=\zeta _n+2\theta _{n1}+2\mathrm{arg}(s_{n1})\mathrm{arg}(\rho _{n1}),$$

and, supposing that the phase shift $`\mathrm{\Psi }_n`$ is a random variable, uniformly distributed over a large interval, performed averaging over $`\mathrm{\Psi }_n`$ directly in Eq. (53), having replaced the term with $`\mathrm{cos}(\mathrm{\Psi }_n)`$ by zero. One can easily verify that the change of variables, $`|\rho _n|^2=\mathrm{tanh}(y_n)`$, $`|r_n|^2=\mathrm{tanh}(z_n)`$, transformes the set of recurrence relations (53) with $`\mathrm{cos}(\mathrm{\Psi }_n)`$ replaced by zero, to the equations $`y_n=y_{n1}+z_n`$, whose obvious solution is $`y_n=_{k=1}^nz_k`$. Thus we arrive at the following formula for the resistance of a chain with completely random distances between scatterers, uniformly distributed over sufficiently large intervals:

$$_n𝒩_n=\frac{1}{2}\underset{k=1}{\overset{n}{}}\frac{1+|r_k|^2}{1|r_k|^2}\frac{1}{2}.$$
(54)
Here $`|r_k|^2`$ is the energy reflection coefficient from the $`k`$th obstacle (barrier) and $`n`$ is the number of obstacles. If $`|r_n|^2=|r|^2=const`$, then (54) goes to the Landauer formula

$$_n𝒩_n=\frac{1}{2}\left(\frac{1+|r|^2}{1|r|^2}\right)^n\frac{1}{2}.$$
(55)
It is worth comparing Eq. (55) with Eq. (3), paying attention to different powers of the reflection coefficient $`|r|`$. For $`n|r|^2>1`$ formula (55) yields

$$𝒩_n^{(rand)}\frac{1}{2}\mathrm{exp}(2n|r|^2).$$
(56)
Consequently, totally random phase fluctuations without systematic phase shift must also lead to an exponential growth of the photon number (chain resistance), although with much more slow rate than in the resonance case, if $`|r|1`$. For example, for $`|r|=0.01`$ one needs about 500 pulses to create 5000 photons under the resonance conditions, according to Eq. (42), whereas Eq. (56) suggests that the same number of photons can appear after about $`510^4`$ pulses in the case of completely random fluctuations (with zero systematic phase shift).

However, removing the term proportional to $`\mathrm{cos}(\mathrm{\Psi }_n)`$ from the right-hand side of Eq. (53) is not well justified operation, because averaging over fluctuations of phase $`\mathrm{\Psi }_n`$ must be performed in the final expression for $`𝒩_n`$, and not at some intermediate steps (as soon as $`f(x)f(x)`$). Moreover, if the random variable $`\mathrm{\Psi }_n`$ is uniformly distributed over some interval $`(\chi ,\chi )`$, then the average value $`\mathrm{cos}(\mathrm{\Psi }_n)=\mathrm{sin}(\chi )/\chi `$ is equal (or close) to zero either if the ratio $`\chi /\pi `$ is an integer or if $`\chi 1`$. In Fig. 6 we show the results of numerical calculations of the dependence of the natural logarithm $`\mathrm{ln}𝒩_n`$ of the average number of photons created after $`n`$ pulses with completely random phase $`\chi _k`$ (defined by Eq. (50)), uniformly distributed over some large interval $`(\chi ,\chi )`$, for $`\chi >1`$ and $`\theta =0`$ (i.e., in the absence of systematic shifts). For each choice of $`\chi `$, averaging was performed over $`20`$ realizations.

Qualitatively, we see an agreement with Landauer’s predictions: the growth of the number of photons (or equivalent resistance) is close to exponential, although the necessary number of pulses is much greater than in the case of $`\chi <1`$ illustrated in the Figs. 2-5. However, the results of numerical experiments are close to the asymptotical Landauer’s dependence (56), $`\mathrm{ln}𝒩_n=2n|r|^2\mathrm{ln}2`$ (given by the dashed straight line labeled as “Landauer”), only for integral ratios $`m=\chi /\pi `$ (take into account the logarithmic scale), whereas the inclinations of curves corresponding to half-integral values of the ratio $`\chi /\pi =m+1/2`$ (when $`|\mathrm{cos}(\chi _n)|=[\pi (m+1/2)]^1`$) are much bigger for small and moderate values of $`m`$, approximating Landauer’s limit only when $`m>100`$.

The results of this section are in qualitative agreement with other studies, e.g.,  (where quasiperiodic motions of the boundaries were considered),  (parametric resonance in periodic paraxial optical systems), or  (devoted to quantum oscillators with fluctuating frequencies). A possibility of generating coherent phonons in solids with time-dependent lattice strains was studied in . Mathematically, this problem is reduced to quantum oscillator with time-dependent (fluctuating) frequency. Another mechanism of producing coherent phonons by femtosecond laser pulses, which is mathematically equivalent to the problem of quantum oscillator with time-dependent force, was considered recently in .

## 5 Conclusion

We have considered a simple model of photon generation from vacuum due to nonstationary (dynamical) Casimir effect in a cavity with nonequidistant spectrum of eigenfrequencies, which is reduced to the problem of quantum harmonic oscillator with time-dependent frequency. We have shown that periodic displacements of the (effective) boundary must result in exponential growth of the mean number of quanta (photons) in the selected mode under certain resonance conditions. The oscillations of the boundary need not to be monochromatic, but they must be close to periodic. However, although the period of motion of the boundary, $`T_b`$, should be adjusted to the period of field oscillations, $`T_f`$, and to some phase, depending on the concrete form of pulse, $`T_b`$ may be much greater than $`T_f`$. This result seems to be very important from the point of view of facilitating performing experiments. The concrete form of the trajectory of the wall enters the final result only through effective reflection and transmission coefficients from some “barriers” in the time dependence of the eigenfrequency.

We have studied the influence of deviations from the strict resonance, caused by systematic and random perturbations of the optimal resonance trajectory of the boundary. It appears that amplitude perturbations are less important than the phase ones. In turn, among the phase perturbations, the most dangerous for preserving the regime of photon generation are systematic deviations, which should not exceed rather low level of the order of the frequency modulation depth. On the contrary, random phase fluctuations can have much greater amplitude (even of the order of unity for zero systematic deviations) without qualitative changes in the behaviour of the system. This circumstance also seems to be important for planning future experiments, showing that some requirements might be not so rigid as one could suppose, thus diminishing the number of technical problems to be solved.

The results obtained can be easily generalized to the case when the initial state of the field mode was not vacuum, but a thermal state. Actually, one should simply multiply the right-hand side of Eq. (5) by the factor $`1+2n_{th}`$, where $`n_{th}`$ is the mean number of thermal photons in the initial state .

However, the presented study should be considered only as a qualitative model, because many important things have not been taken into account. For example, we considered the excitation of a single selected mode, neglecting its possible interactions with other modes. Such an approach seems to be justified for cavities with nonequidistant spectra, when the boundary performs harmonic oscillations at the resonance frequency . But the spectrum of anharmonic oscillations contains many frequencies. On the one hand, this fact, perhaps, explains why photons can be generated even if the period of displacements of the boundary is greater than the period of the field oscillations. On the other hand, in the presence of many frequencies some other modes can be excited, too. This question needs a detailed investigation. Also, the effects of polarization of true electromagnetic field (not its scalar model) could be important, as was shown in other examples . And, of course, the problem of consistent account of losses in cavities with moving nonideal mirrors must be solved (perhaps, following the lines generalizing the approaches of Refs. ).

## Acknowledgement

The authors are grateful to the Brazilian agency CNPq for the support.