text:
# Study of Scalar Top Quarks in the Neutralino and Chargino Decay Channel

## Introduction

The study of the scalar top quarks is of particular interest, since the lighter stop mass eigenstate is likely to be the lightest scalar quark in a supersymmetric theory. The mass eigenstates are $`m_{\stackrel{~}{\mathrm{t}}_1}`$ and $`m_{\stackrel{~}{\mathrm{t}}_2}`$ with $`m_{\stackrel{~}{\mathrm{t}}_1}<m_{\stackrel{~}{\mathrm{t}}_2}`$, where $`\stackrel{~}{\mathrm{t}}_1=\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}\stackrel{~}{\mathrm{t}}_\mathrm{L}+\mathrm{sin}\theta _{\stackrel{~}{\mathrm{t}}}\stackrel{~}{\mathrm{t}}_\mathrm{R}`$ and $`\stackrel{~}{\mathrm{t}}_2=\mathrm{sin}\theta _{\stackrel{~}{\mathrm{t}}}\stackrel{~}{\mathrm{t}}_\mathrm{L}+\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}\stackrel{~}{\mathrm{t}}_\mathrm{R}`$ with the mixing angle $`\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}`$. We study the experimental possibilities to determine $`m_{\stackrel{~}{\mathrm{t}}_1}`$ and $`\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}`$ at a high-luminosity $`\mathrm{e}^+\mathrm{e}^{}`$ linear collider like the TESLA project $`^\mathrm{?}`$ with polarized $`\mathrm{e}^+`$ and $`\mathrm{e}^{}`$ beams.

The simulated production process is $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mathrm{t}}_1\overline{\stackrel{~}{\mathrm{t}}_1}`$ with two decay modes $`\stackrel{~}{\mathrm{t}}_1\stackrel{~}{\chi }^0`$c and $`\stackrel{~}{\mathrm{t}}_1\stackrel{~}{\chi }^+`$b. A 100% branching fraction in each decay mode is simulated. The first scalar top decay into a c-quark and the lightest neutralino results in a signature of two jets and large missing energy. The second investigated stop decay mode leads also to large missing energy and further jets from the chargino decay. The neutralino channel is dominant unless the decay into a chargino is kinematically allowed. Details of the event simulation with SGV $`^\mathrm{?}`$ tuned for a TESLA detector $`^\mathrm{?}`$ are given in Ref. $`^\mathrm{?}`$. The signals and a total of 16 million Standard Model background events are simulated (Table 1) for $`=500`$ fb-1.

## Neutralino Channel

The reaction $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mathrm{t}}_1\overline{\stackrel{~}{\mathrm{t}}_1}\stackrel{~}{\chi }^0\mathrm{c}\stackrel{~}{\chi }^0\overline{\mathrm{c}}`$ has been studied for a 180 GeV scalar top and a 100 GeV neutralino. After a preselection, 278377 background events remain $`^{\mathrm{?},\mathrm{?}}`$. In order to separate the signal from the background, the following selection variables are defined: visible energy, number of jets, thrust value and direction, number of clusters, transverse and parallel imbalance, acoplanarity and invariant mass of two jets $`^{\mathrm{?},\mathrm{?}}`$. An Iterative Discriminant Analysis (IDA) $`^\mathrm{?}`$ optimized the selection. For unpolarized beams and 12% efficiency, 400 background events are expected.

The polarization of the $`\mathrm{e}^+`$ and $`\mathrm{e}^{}`$ beams at a future linear collider offers the opportunity to enhance or suppress the left- or right-handed couplings of the scalar top signal and to determine mass and mixing angle independently. The production cross section of each background process depends differently on the polarization. It is therefore important for a high-statistics analysis to study the expected background channels individually. The expected cross sections $`^{\mathrm{?},\mathrm{?}}`$ are given in Table 2 for different beam polarization states. The IDA analysis was repeated for $`0.9`$ and $`0.9`$ polarization $`^\mathrm{?}`$bbbFor a polarization of $`0.9`$, 95% of the $`\mathrm{e}^{}`$ are left-polarized. In the previous analyses $`^{\mathrm{?},\mathrm{?},\mathrm{?},\mathrm{?}}`$ it was assumed that only 90% of the $`\mathrm{e}^{}`$ were polarized.. Here, we recalculate all rates for the new machine polarization. Figure 1 shows the number of background events as a function of the signal efficiency for $`0.8/0.6`$ (left) and $`0.8/0.6`$ (right) polarization. For 12% detection efficiency, 1194 background events are expected leading to $`\sigma _{\mathrm{left}}=81.8\pm 1.3`$ fb, and 208 background events giving $`\sigma _{\mathrm{right}}=76.4\pm 1.2`$ fb, where $`\mathrm{\Delta }\sigma /\sigma =\sqrt{N_{\mathrm{signal}}+N_{\mathrm{background}}}/N_{\mathrm{signal}}`$.

## Chargino Channel

The reaction $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mathrm{t}}_1\overline{\stackrel{~}{\mathrm{t}}_1}\stackrel{~}{\chi }^+\mathrm{b}\stackrel{~}{\chi }^{}\overline{\mathrm{b}}\stackrel{~}{\chi }^0\mathrm{W}^+\mathrm{b}\stackrel{~}{\chi }^0\mathrm{W}^{}\overline{\mathrm{b}}`$ has been studied with focus on the hadronic W decays for a 180 GeV scalar top, a 150 GeV chargino and a 60 GeV neutralino, where the chargino decays 100% into a W boson and a neutralino. A preselection similar to that for the study of the neutralino channel is applied and 209051 background events remain $`^\mathrm{?}`$. In order to separate the signal from the background, the following selection variables are defined: visible energy, number of jets, thrust value, number of clusters, transverse and parallel imbalance, and the isolation angle of identified leptons. For an unpolarized beam, Fig. 2 shows the final IDA output variable and the resulting number of background events as a function of the signal efficiency. For 12% efficiency, only 20 background events are expected. Allowing a background rate of 400 events, as in the neutralino channel, the efficiency is 44%, from which we derive the relative error on the cross section to be 0.75%. In this case the number of expected signal events is much larger than the expected background, thus no separate tuning of the IDA for left- and right-polarization is required. The background is neglected for the determination of mass and mixing angle. Note that the dependence of the background rates on the polarization has to be taken into account for stop masses closer to the kinematic threshold.

## Results

For the neutralino and chargino decay channel of scalar top quarks, we have determined the expected Standard Model background rate as a function of the signal efficiency. The total simulated background of about 16 million events is largely reduced, which allows a precision measurement of the scalar top production cross section with a relative error of better than 2% in the neutralino channel and about 1% in the chargino channel. Based on experiences gained at LEP, we expect that detection efficiencies for other mass combinations are similar as long as the mass difference between the scalar top and the neutralino is larger than about 20 GeV. Figure 3 shows the corresponding error bands and the error ellipse in the $`m_{\stackrel{~}{\mathrm{t}}_1}`$ – $`\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}`$ plane for both decay channels for $`0.8/0.6`$ left- and right-polarization of the $`\mathrm{e}^{}/\mathrm{e}^+`$ beams and a luminosity of 500 fb-1 each. The statistical errors are a factor 7 better in the neutralino channel and about a factor 14 better in the chargino channel than reported previously $`^{\mathrm{?},\mathrm{?},\mathrm{?},\mathrm{?}}`$, and improve further when in addition e+ polarization is included. Detailed results are given in Table 3. The highest statistical precision is obtained in the chargino channel with an error $`\mathrm{\Delta }m_{\stackrel{~}{\mathrm{t}}_1}=0.4`$ GeV for $`m_{\stackrel{~}{\mathrm{t}}_1}=180`$ GeV and $`\mathrm{\Delta }\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}=0.003`$ $`\mathrm{cos}\theta _{\stackrel{~}{\mathrm{t}}}=0.570`$. Based on the experience from direct searches at LEP, the systematic errors on the event selection are less than 1%; precise investigations require the detailed detector layout and a full simulation. The stop generator $`^\mathrm{?}`$ has been interfaced with the SIMDET $`^\mathrm{?}`$ simulation to allow an independent test of the detector simulation and related systematic errors. Another uncertainty could arise from the luminosity measurement, the measurement of the polarization, and the theoretical uncertainty of the production cross section, which will be determined in the future. A high-luminosity linear collider with the capability of beam polarization has a great potential for precision measurements in the scalar quark sector predicted by Supersymmetric theories.

## References


format: text

text:
# Thermohydrodynamics for a van der Waals fluid

## I Introduction

Boiling of water in a pot, the formation of a cloud, or the rise of bubbles in a pint of beer are fascinating everyday experiences that involve liquid-vapor coexistence in non-trivial flow situations. Knowledge of the dynamics of this phenomena is crucial in less pleasant situations like the prediction and control of nuclear accidents in refrigerated nuclear reactors. The complexity of the problem requires the aid of computer simulations in order to extract practical information. Engineering conventional computational fluid mechanics approaches have relied on effective hydrodynamic equations in which the presence of bubbles or drops is taken into account through the introduction of void or vapor fraction which require empirical constitutive equations not always available . In these much coarse grained approaches the detailed interface dynamics of bubbles and droplets is not available.

There is a great recent interest in the computer simulations of interfacial dynamics for liquid-vapor coexistence. Prove of that is the large variety of techniques used to tackle the problem: Lattice-Boltzmann method , Direct Montecarlo method , Smoothed Particle Hydrodynamics (SPH), Finite Difference discretization of Navier-Stokes equations ,, and the Volume Of Fluid method . Some of these techniques suffer from ad-hoc assumptions, more prominently isothermal behavior. Others require quite specific model technicalities in order to get sensible results or the specification of boundary conditions in complex topologies. It seems, therefore, appropriate to review here the theoretical foundation of the hydrodynamic equations for the van der Waals fluid.

Although hydrodynamics equations for the van der Waals fluid were posed as early as 1901 by Korteveg ,, a clarification of the structure and thermodynamic consistency of the equations is convenient in order to avoid potential problems. In the late 60’s, for example, a variant of the theory was introduced by Kawasaki , following the ideas of van Kampen . It can be shown that Kawasaki theory is thermodynamically inconsistent. In Kawasaki theory, the long range attractive mean field potential between molecules acts as an external force in the momentum equation. A gradient approximation of this mean field potential leads to purely local equations involving third order spatial derivatives of the density field. Such a theory was proposed by Felderhof who constructed a thermodynamically consistent set of hydrodynamic equations for a van der Waals fluid . However, Felderhof theory deals with the inviscid fluid and no dissipation nor fluctuations where included in his theory so a generalization of his theory seems appropriate. For an excellent review of diffuse interface models described by third order derivative terms see Ref. . On the other hand, the discretization of third order derivatives is subject to high numerical errors. From a computational point of view, a treatment of the diffuse interfaces though a non-local integral term as in Kawasaki theory might represent an advantage worth exploring.

The essential theoretical tool for describing equilibrium properties of liquid-vapor interfaces is the Density Functional Theory (DFT) . A well-known local realization of the DFT known as the Density Gradient Theory successfully describes the equilibrium properties of liquid-vapor interfaces like surface tension and density profiles, and the correct relation between different scaling exponents near the critical point . This is due to the relatively smooth behavior of these interfaces as compared with fluid-solid interfaces for which a non-local expression for the free energy functional is needed . The work presented in this paper can be understood as a generalization of the local DFT to non-equilibrium situations in which an entropy functional instead of a free energy functional is used to describe both flow and thermal transport .

The strategy we follow is to derive the entropy functional from microscopic principles following ideas of van Kampen, who derived a coarse grained free energy functional for the van der Waals gas . This entropy functional is one of the building blocks of a recently developed framework for non-equilibrium thermodynamics known as GENERIC . This two generator formalism ensures thermodynamic consistency of its dynamic equations. Under a set of well-defined assumptions, the GENERIC formalism can be derived from first principles by means of a projection operator technique . It is therefore not surprising that all well-established dynamic equations for non-equilibrium systems fit into the GENERIC formalism: Relativistic and non-relativistic hydrodynamics, kinetic theory of gases (the Boltzmann equation) and of polymer systems in non-isothermal situations, chemical reactions, just to mention a few, have the GENERIC structure . Significant new dynamical models for a number of complex systems, such as polymer/surface interactions, liquid crystals, or polymer blends have already been obtained, and even new cosmological models have been developed. In particular, the recognition of the role of the reversible dynamics in the system evolution equations has led to interesting new conclusions regarding closure approximations and polymer reptation models. It seems, therefore, sensible to formulate a thermohydrodynamic theory for liquid-vapor coexistence according to the GENERIC framework. The basic bonus for this procedure is that strict respect for the first and second law of thermodynamics is guarded.

The paper is distributed as follows. In Section II we present a brief summary of the GENERIC formalism and of the entropy functional for a van der Waals fluid. The entropy functional is computed explicitly in the Appendix VI. This appendix can be regarded as providing solid microscopic ground to the plausible assumptions made by Felderhof when postulating the entropy functional. In Section III we propose the hydrodynamic equations of a van der Waals fluid according to the GENERIC formalism and compare it with the theory of Kawasaki. In Section IV the local gradient theory is presented and compared with Felderhof’s theory. Some concluding remarks are given in Section V, while some further information is given in Appendix VII, and VIII.

## II The GENERIC framework and the entropy

A brief summary of the GENERIC structure is presented and we refer to the original references for further details . The first essential step in the description of a system is the selection of the proper relevant variables $`x`$ used to describe the system at a given level of description. The GENERIC dynamical equation for $`x`$ is

$$\dot{x}=L(x)E(x)+M(x)S(x).$$
(1)
The first term in the right hand side produces the reversible part of the dynamics whereas the second term is responsible for the irreversible dissipative dynamics. Here, $`E(x),S(x)`$ are the energy and entropy of the system expressed in terms of the variables $`x`$, $``$ is the gradient operator in $`x`$-space, and $`L(x),M(x)`$ are matrices that satisfy the following degeneracy requirements,

$$LS=0,ME=0.$$
(2)
In addition, $`L`$ is antisymmetric (this guarantees that energy is conserved, i.e. the First Law) and satisfies the stringent Jacobi property . $`M`$ is a positive definite symmetric matrix (this guarantees that the entropy is a nondecreasing function of time, i.e. the Second Law). If the system presents dynamical invariants $`I(x)`$ different from the total energy, then further restrictions on the form of $`L`$ and $`M`$ are required

$$ILE=0,IMS=0.$$
(3)
These conditions ensure that $`dI/dt=0`$. The deterministic equations (1) are, actually, an approximation in which thermal fluctuations are neglected. If thermal fluctuations are not neglected, the dynamics is described by a Fokker-Planck equation (FPE) that governs the probability distribution function $`\rho =\rho (x,t)`$ . This FPE has as equilibrium solution the Einstein’s distribution function generalized to take into account the presence of dynamical invariants $`E_0,I_0`$ , this is

$$\rho ^{\mathrm{eq}}(x)=\delta (E(x)E_0)\delta (I(x)I_0)\mathrm{exp}\{S(x)/k_B\},$$
(4)
where $`k_B`$ is Boltzmann’s constant.

Our approach for constructing the entropy function for a liquid-vapor system is to compute microscopically $`\rho ^{\mathrm{eq}}(x)`$ and, through Eqn. (4), identify the entropy functional. That is, we compute the joint probability that an extended simple fluid in equilibrium has a particular realization of the mass, momentum and internal energy density fields. We follow a cell-method first used by van Kampen in order to derive the equilibrium probability that a simple fluid has a particular realization of the mass density field . The explicit details of the calculation are given in Appendix VI. Physical space is divided in $``$ cells of volume $`𝒱_\mu `$ and the mass $`M_\mu `$, momentum $`𝐏_\mu `$ and internal energy $`ϵ_\mu `$ of cell $`\mu `$ in terms of microscopic coordinates of the molecules is considered. The internal energy contains the potential energy of interaction between particles within the same cell and the energy of interaction with particles of other cells. Provided that the microscopic coordinates are distributed according to the microcanonical ensemble we compute in Appendix VI the joint probability $`P[M,𝐏,ϵ]`$ following similar steps as in Ref. . An essential step in the derivation is the assumption that the molecular potential has a repulsive hard core $`\varphi ^{\mathrm{hc}}(r)`$ and a long range attractive tail $`\varphi ^l(r)`$. The interaction of particles between different cells is taken in mean field approximation and it involves only $`\varphi ^l(r)`$. The final outcome is

$`P[M,𝐏,ϵ]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }_0}}\delta \left({\displaystyle \underset{\mu }{\overset{}{}}}M_\mu M_0\right)\delta \left({\displaystyle \underset{\mu }{\overset{}{}}}𝐏_\mu \right)`$ (5)
$`\times `$ $`\delta \left({\displaystyle \underset{\mu }{\overset{}{}}}\left({\displaystyle \frac{𝐏_\mu ^2}{2M_\mu }}+ϵ_\mu \right)E_0\right)\mathrm{exp}\left\{{\displaystyle \frac{S[M,ϵ]}{k_B}}\right\}.`$ (6)
We recognize the dynamical invariants (total mass $`M_0`$, momentum $`𝐏_0=\mathrm{𝟎}`$ and energy $`E_0`$) in the conserving delta functions. The entropy functional is given by

$$S[M,ϵ]=\underset{\mu }{}S^{\mathrm{hc}}(ϵ_\mu \overline{\varphi }_\mu ,M_\mu ,𝒱_\mu ).$$
(7)
Here, $`S^{\mathrm{hc}}(ϵ,M,V)`$ is the entropy function of an isolated system of $`N=M/m`$ molecules ($`m`$ is the mass of a molecule) interacting with only the hard core potential $`\varphi ^{\mathrm{hc}}(r)`$ in a volume $`V`$. Note the non-trivial appearance of the mean field interaction energy $`\overline{\varphi }_\mu `$, which is defined by

$$\overline{\varphi }_\mu =\frac{1}{2m^2}\underset{\nu }{\overset{}{}}\varphi ^l(R_{\mu \nu })M_\mu M_\nu ,$$
(8)
where $`R_{\mu \nu }`$ is the distance between cell centers and the prime denotes $`\mu \nu `$. This mean field potential involves only the long ranged part of the molecular potential. Integration of (5) with respect to $`N,ϵ`$ and using a steepest descent approximation (see Appendix VII) leads to van Kampen’s expression,

$$P[M]=\frac{1}{Z_0}\mathrm{exp}\left\{\beta F[M]\right\}\delta \left(\underset{\mu }{\overset{}{}}M_\mu M_0\right),$$
(9)
where the coarse-grained free energy functional is given by

$$F[M]=\underset{\mu }{\overset{}{}}F^{\mathrm{hc}}(\beta ,M_\mu ,𝒱_\mu )+\beta \overline{\varphi }_\mu .$$
(10)
Here, $`F^{\mathrm{hc}}(\beta ,M_\mu ,𝒱_\mu )`$ is the free energy of a system of hard core particles at temperature $`\beta ^1/k_B=2K^{}/DM`$ ($`K^{}`$ is the most probable value of the kinetic energy in the system). Both Eqns. (5) and (9) admit a continuum notation provided that the cells are large enough in order that the entropy and free energy are first order functions of their arguments

$`S[\rho _𝐫,ϵ_𝐫]`$ $`=`$ $`{\displaystyle 𝑑𝐫s^{\mathrm{hc}}(ϵ_𝐫\overline{\varphi }_𝐫,\rho _𝐫)},`$ (11)
$`F[\rho _𝐫]`$ $`=`$ $`{\displaystyle 𝑑𝐫\left(f^{\mathrm{hc}}(\rho _𝐫,\beta )+\overline{\varphi }_𝐫\right)},`$ (12)
where $`ϵ_𝐫=ϵ_\mu /𝒱_\mu `$, $`\overline{\varphi }_𝐫=\overline{\varphi }_\mu /𝒱_\mu `$, $`\rho _𝐫=M_\mu /𝒱_\mu `$, and $`s^{\mathrm{hc}}(ϵ_𝐫\overline{\varphi }_𝐫,\rho _𝐫)`$ is the entropy density of a hard core system, while $`f^{\mathrm{hc}}(\rho _𝐫,\beta )`$ is the corresponding free energy density. The continuum form for the mean field potential energy (8) is

$$\overline{\varphi }_𝐫=\rho _𝐫𝑑𝐫^{}\frac{\varphi ^l(|𝐫𝐫^{}|)}{2m^2}\rho _𝐫^{}$$
(13)
The main benefit of the microscopic derivation presented in Appendix VII is the clear separation at the macroscopic level of the effects of the short and long ranged part of the microscopic potential, in particular the realization that the entropy of the full system is given directly in terms of the entropy of the hard core system, suitably evaluated at the corrected internal energy, Eqns. (7) or (11).

We now consider the particular situation in which the density field varies slowly in the scale of the range of $`\varphi ^l(r)`$. In this case, we can Taylor expand $`\rho _𝐫^{}`$ around $`𝐫`$ and retain up to second order in gradients in order to get

$$\overline{\varphi }_𝐫a\rho _𝐫^2\frac{c}{2}\rho _𝐫^2\rho _𝐫,$$
(14)
where we introduced $`a,c>0`$ through

$$a=𝑑𝐫\frac{\varphi ^l(r)}{2m^2},c=𝑑𝐫r^2\frac{\varphi ^l(r)}{3m^2}.$$
(15)
Substitution of (14) into (12) leads to

$$F[\rho _𝐫]=𝑑𝐫\left(f^{\mathrm{hc}}(\rho _𝐫,\beta )a\rho ^2+\frac{c}{2}(\rho _𝐫)^2\right),$$
(16)
which is the usual local gradient density expression for the free energy functional, also known as the Cahn-Hilliard or Ginzburg-Landau free energy functional for liquid surfaces . The corresponding entropy functional (11), thus, has been computed at the same level of approximation. The coefficient $`c`$ provides the overall magnitude of the surface tension coefficient .

The energy function given by

$$E(M,𝐏,ϵ)=\underset{\mu }{\overset{}{}}\left(\frac{𝐏_\mu ^2}{2M_\mu }+ϵ_\mu \right),$$
(17)
can also be written in continuum form as an energy functional

$$E[\rho _𝐫,𝐠_𝐫,ϵ_𝐫]=\left[\frac{1}{2}\frac{𝐠_𝐫^2}{\rho _𝐫}+ϵ_𝐫\right]𝑑𝐫.$$
(18)
where $`𝐠_𝐫=𝐏_\mu /𝒱_\mu `$ is the momentum density field.

## III Hydrodynamic equations for a van der Waals fluid

In this section, we construct the hydrodynamic equations for a van der Waals fluid following the GENERIC formalism.

Due to the form in which the internal energy appears in the entropy, it proves convenient to use as the proper hydrodynamic variable the “intrinsic” internal energy $`ϵ_𝐫^i=ϵ_𝐫\overline{\varphi }_𝐫`$ which involves the interaction of molecules through the hard core part only. The energy (18) and the entropy (11) expressed in terms of these hydrodynamic variables are then

$`E[\rho ,𝐠,ϵ^i]`$ $`=`$ $`{\displaystyle \left[\frac{1}{2}\frac{𝐠_𝐫^2}{\rho _𝐫}+ϵ_𝐫^i+\overline{\varphi }_𝐫\right]𝑑𝐫},`$ (19)
$`S[\rho ,𝐠,ϵ^i]`$ $`=`$ $`{\displaystyle s^{hc}(\rho _𝐫,ϵ_𝐫^i)𝑑𝐫}.`$ (20)
We have now the two basic building blocks for the GENERIC formulation of the dynamic equations (1). The derivatives of the energy and entropy with respect to the relevant variables become functional derivatives with respect to the hydrodynamic variables, this is,

$$\frac{E}{x}=\left(\begin{array}{c}\frac{\delta E}{\delta \rho _𝐫}\\ \\ \frac{\delta E}{\delta 𝐠_𝐫}\\ \\ \frac{\delta E}{\delta ϵ_𝐫^i}\end{array}\right)=\left(\begin{array}{c}\frac{1}{2}𝐯_𝐫^2\frac{2\overline{\varphi }_𝐫}{\rho _𝐫}\\ \\ 𝐯_𝐫\\ \\ 1\end{array}\right),$$
(21)
and

$$\frac{S}{x}=\left(\begin{array}{c}\frac{\delta S}{\delta \rho _𝐫}\\ \\ \frac{\delta S}{\delta 𝐠_𝐫}\\ \\ \frac{\delta S}{\delta ϵ_𝐫^i}\end{array}\right)=\left(\begin{array}{c}\frac{\mu ^{\mathrm{hc}}(\rho _𝐫,ϵ_𝐫^i)}{T^{\mathrm{hc}}(\rho _𝐫,ϵ_𝐫^i)}\\ \\ 0\\ \\ \frac{1}{T^{\mathrm{hc}}(\rho _𝐫,ϵ_𝐫^i)}\end{array}\right).$$
(22)
The velocity field is $`𝐯_𝐫=𝐠_𝐫/\rho _𝐫`$, and the local temperature $`T^{\mathrm{hc}}`$ and chemical potential $`\mu ^{\mathrm{hc}}`$ per unit mass are

$$T_𝐫^{\mathrm{hc}}=\left(\frac{s^{\mathrm{hc}}}{ϵ_𝐫^i}\right)^1,\frac{\mu _𝐫^{\mathrm{hc}}}{T_𝐫^{\mathrm{hc}}}=\frac{s^{\mathrm{hc}}}{\rho _𝐫}.$$
(23)
Note that these equations of state are those of the hard core system. By focusing on the reversible part of the dynamics, we must construct now the matrix $`L`$. A reasonable proposal is to use the same expression for $`L`$ as in the case of hydrodynamics of a simple fluid. Note that if we neglect the interfacial potential energy $`\overline{\varphi }_𝐫`$ we should recover the hydrodynamics of a simple hard core fluid. We propose, therefore, for the matrix $`LL_{\mathrm{𝐫𝐫}^{}}`$ the following one

$$\left(\begin{array}{ccc}0& \rho _𝐫^{}^{}\delta _{\mathrm{𝐫𝐫}^{}}& 0\\ & & \\ \rho _𝐫^{}\delta _{\mathrm{𝐫𝐫}^{}}& 𝐠_𝐫^{}^{}\delta _{\mathrm{𝐫𝐫}^{}}+^{}\delta _{\mathrm{𝐫𝐫}^{}}𝐠_𝐫& ϵ_𝐫^i^{}\delta _{\mathrm{𝐫𝐫}^{}}+P_𝐫^{}^{\mathrm{hc}}^{}\delta _{\mathrm{𝐫𝐫}^{}}\\ & & \\ 0& ϵ_𝐫^{}^i^{}\delta _{\mathrm{𝐫𝐫}^{}}+P_𝐫^{\mathrm{hc}}^{}\delta _{\mathrm{𝐫𝐫}^{}}& 0\end{array}\right).$$
(24)
Here, $`^{}=/𝐫^{}`$ and $`\delta _{\mathrm{𝐫𝐫}^{}}=\delta (𝐫𝐫^{})`$ is Dirac’s delta function, $`P_𝐫^{\mathrm{hc}}=P^{\mathrm{hc}}(\rho _𝐫,ϵ_𝐫^i)`$ is the pressure field as a function of the mass and intrinsic internal energy densities.

As has been shown in Ref. , the matrix $`L`$ satisfies the condition $`LS=0`$ because of the Gibbs-Duhem relationship (for the hard core system). The fulfillment of $`LS=0`$ is the basic reason for using $`P_𝐫^{\mathrm{hc}}`$ in Eqn. (24) instead of any other possible expression for the pressure. The matrix $`L`$ also satisfies the stringent Jacobi identity . With the explicit form (24) for the $`L`$ matrix and the gradient of the energy in Eqn. (21) the reversible part $`LE`$ of the equations of motion for the hydrodynamic variables are readily constructed

$`_t\rho _𝐫`$ $`=`$ $`\rho _𝐫𝐯_𝐫,`$ (25)
$`_t𝐠_𝐫`$ $`=`$ $`𝐯_𝐫𝐠_𝐫P^{\mathrm{hc}}(ϵ_𝐫^i,\rho _𝐫)\rho _𝐫{\displaystyle \frac{2\overline{\varphi }_𝐫}{\rho _𝐫}},`$ (26)
$`_tϵ_𝐫^i`$ $`=`$ $`ϵ_𝐫^i𝐯_𝐫P^{\mathrm{hc}}(ϵ_𝐫^i,\rho _𝐫)𝐯_𝐫.`$ (27)
Eqns. (27) are identical to the Euler equations of an inviscid (hard cored) fluid except for the additional term in the momentum equation that involves the long range part of the molecular potential.

Concerning the irreversible part of the dynamics $`MS`$, we observe that the same matrix $`MM_{\mathrm{𝐫𝐫}^{}}`$ that corresponds to the usual hydrodynamic equations is perfectly adequate. For simplicity, we assume that there is no extra dissipation due to the presence of interfaces. The matrix $`M`$ is thus given by

$$\left(\begin{array}{ccc}0& 0& 0\\ & & \\ 0& (^{}+\mathrm{𝟏}^{})\eta _𝐫T_𝐫^{\mathrm{hc}}\delta _{\mathrm{𝐫𝐫}^{}}+2^{}\widehat{\kappa }_𝐫T_𝐫^{\mathrm{hc}}\delta _{\mathrm{𝐫𝐫}^{}}& \eta _𝐫T_𝐫^{\mathrm{hc}}\dot{\gamma }\delta _{\mathrm{𝐫𝐫}^{}}+\widehat{\kappa }_𝐫T_𝐫^{\mathrm{hc}}\mathrm{tr}\dot{\gamma }\delta _{\mathrm{𝐫𝐫}^{}}\\ & & \\ 0& ^{}\eta _𝐫T_𝐫^{\mathrm{hc}}\dot{\gamma }\delta _{\mathrm{𝐫𝐫}^{}}+^{}\widehat{\kappa }_𝐫T_𝐫^{\mathrm{hc}}\mathrm{tr}\dot{\gamma }\delta _{\mathrm{𝐫𝐫}^{}}& \frac{1}{2}\eta _𝐫T_𝐫^{\mathrm{hc}}\dot{\gamma }:\dot{\gamma }\delta _{\mathrm{𝐫𝐫}^{}}+^{}\lambda _𝐫T_𝐫^{\mathrm{hc2}}\delta _{\mathrm{𝐫𝐫}^{}}+\frac{1}{2}\widehat{\kappa }_𝐫T_𝐫^{\mathrm{hc}}(\mathrm{tr}\dot{\gamma })^2\delta _{\mathrm{𝐫𝐫}^{}}\end{array}\right),$$
(28)
where $`\dot{\gamma }=𝐯_𝐫+(𝐯_𝐫)^T`$ is the symmetrized velocity gradient tensor. This matrix $`M`$ satisfies the degeneracy $`ME=0`$ as in the usual hydrodynamic case , as the only modification of $`E`$ is in the term in the density component. We remark, however, that the transport coefficients (shear viscosity $`\eta _𝐫`$, bulk viscosity $`\kappa _𝐫=\widehat{\kappa }_𝐫+2\eta _𝐫/3`$ and thermal conductivity $`\lambda _𝐫`$) that appear in $`M`$ should be taken as state dependent and are expected to be strongly varying functions of the density field, in order to encompass the fact that vapor and liquid have very different transport properties. This is the reason for its space dependence displayed as a subindex $`𝐫`$. The presence of $`T_𝐫^{\mathrm{hc}}`$ in Eqn. (28) instead of other possible choices for the temperature, is dictated by the fact that $`MS`$ should produce the usual dissipative terms in the Navier-Stokes equation (note that $`S`$ in Eqn. (22) involves $`T_𝐫^{\mathrm{hc}}`$).

The final set of hydrodynamic equations for the thermohydrodynamics of the van der Waals fluid are

$`_t\rho _𝐫`$ $`=`$ $`\rho _𝐫𝐯_𝐫,`$ (29)
$`_t𝐠_𝐫`$ $`=`$ $`𝐯_𝐫𝐠_𝐫P^{\mathrm{hc}}(ϵ_𝐫^i,\rho _𝐫)\tau +\rho _𝐫\overline{𝐅}_𝐫,`$ (30)
$`_tϵ_𝐫^i`$ $`=`$ $`ϵ_𝐫^i𝐯_𝐫P^{\mathrm{hc}}(ϵ_𝐫^i,\rho _𝐫)𝐯_𝐫𝐪\tau :𝐯_𝐫,`$ (31)
where the constitutive equations for the viscous stress $`\tau `$ and heat flux are given by

$`\tau _𝐫`$ $`=`$ $`\eta _𝐫\left[𝐯_𝐫+𝐯_𝐫^T\right](\kappa _𝐫2\eta _𝐫/3)(𝐯_𝐫)\mathrm{𝟏},`$ (32)
$`𝐪_𝐫`$ $`=`$ $`\lambda _𝐫T_𝐫^{\mathrm{hc}}.`$ (33)
and the mean field force is given by

$$\overline{𝐅}_𝐫=_𝐫\frac{1}{m^2}𝑑𝐫^{}\overline{\varphi }^l(|𝐫𝐫^{}|)\rho _𝐫^{}$$
(34)
Equations (31) are deterministic equations in which thermodynamic fluctuations are neglected. Thermodynamic fluctuations are very easily included in the GENERIC formalism. One simply needs to take the square root in matrix sense of the $`M`$ matrix, and this provides the amplitude of the noises , in accordance with the Fluctuation-Dissipation theorem. Because the matrix $`M`$ in Eqn. (28) has the same structure of that of a one phase fluid, we can advance that the final fluctuating equations are obtained by simply adding to the stress tensor $`\tau _𝐫`$ and heat flux $`𝐪_𝐫`$ a random counterparts given by

$`\stackrel{~}{\tau }_𝐫`$ $`=`$ $`\sqrt{2T_𝐫^{\mathrm{hc}}\eta _𝐫}(\stackrel{~}{\sigma }_𝐫(t)\mathrm{tr}[\stackrel{~}{\sigma }_𝐫(t)]\mathrm{𝟏})+\sqrt{3T_𝐫^{\mathrm{hc}}\kappa _𝐫}\mathrm{tr}[\stackrel{~}{\sigma }_𝐫(t)]\mathrm{𝟏}`$ (35)
$`\stackrel{~}{𝐪}_𝐫`$ $`=`$ $`T_𝐫^{\mathrm{hc}}\sqrt{2\lambda _𝐫}\stackrel{~}{\zeta }_𝐫(t),`$ (36)
where $`\stackrel{~}{\sigma }_𝐫(t)`$ is a matrix of delta correlated white noises and $`\stackrel{~}{\zeta }_𝐫(t)`$ is a vector of delta correlated white noises (in space and time). Note that the form in which thermal fluctuations appear is similar to the approach proposed by Landau and Lifshitz . However, one should note that in counterdistintion to their approach, here the random stress tensor and heat flux are multiplicative noises which are not Gaussian. What is Gaussian are the white noises $`\stackrel{~}{\sigma }_𝐫(t),\stackrel{~}{\zeta }_𝐫(t)`$. The strong dependence of the transport coefficients on the density due to the phase change has a direct effect on the amplitude of the thermal fluctuations.

Associated to the stochastic differential equations describing hydrodynamic fluctuating variables there is a mathematical equivalent functional Fokker-Planck equation . The GENERIC structure ensures that the equilibrium solution of this Fokker-Planck equation is given by the continuum version of Eqn. (5), as it should.

Kawasaki proposed a set of hydrodynamic equations for a van der Waals fluid on intuitive grounds . Both the continuity equation and the momentum balance equation in (27) coincide with Kawasaki’s postulated equations. On the other hand, he uses the entropy per unit mass $`\widehat{s}_𝐫^{\mathrm{hc}}=s_𝐫^{\mathrm{hc}}/\rho _𝐫`$ instead of the intrinsic internal energy $`ϵ_𝐫^i`$. A standard calculation which uses the definitions (23) and the Euler equation $`T^{\mathrm{hc}}s^{\mathrm{hc}}+\mu ^{\mathrm{hc}}\rho P^{\mathrm{hc}}ϵ^i=0`$ shows that our equations (27) imply $`_t\widehat{s}^{\mathrm{hc}}+𝐯_𝐫\widehat{s}^{\mathrm{hc}}=0`$. This is consistent with the fact that the reversible part of the dynamics should not produce an increase of the entropy in the system. In contrast, in Kawasaki treatment the reversible part of the substantial derivative of the entropy per unit mass is equated to $`\rho _𝐫𝐯_𝐫(\overline{\varphi }_𝐫/\rho _𝐫)`$ (in our notation). It is apparent that the addition of this term violates the conservation of the energy (19) and should not be present for thermodynamic consistency.

## IV Gradient approximation

In this section we consider the gradient approximation (14) as applied to the obtained equations (31). The last two terms of the momentum equation can be re-arranged in the form

$$P_𝐫^{\mathrm{hc}}\rho _𝐫\frac{\overline{\varphi }_𝐫}{\rho _𝐫}=\overline{P}_𝐫+c\rho _𝐫^2\rho _𝐫,$$
(37)
where we have introduced

$$\overline{P}_𝐫=P^{\mathrm{hc}}(\rho _𝐫,ϵ_𝐫^i)a\rho _𝐫^2.$$
(38)
The momentum equation thus becomes

$$_t𝐠_𝐫=𝐯_𝐫𝐠_𝐫\overline{P}(ϵ_𝐫^i,\rho _𝐫)+c\rho _𝐫^2\rho _𝐫.$$
(39)
It is illuminating to use as the set of hydrodynamic variables $`\rho _𝐫,𝐠_𝐫,u_𝐫`$ instead of $`\rho _𝐫,𝐠_𝐫,ϵ_𝐫^i`$. Here, the local internal energy is defined by

$$u_𝐫=ϵ_𝐫^ia\rho _𝐫^2,$$
(40)
which, from a microscopic point of view, represents the potential energy due to local interactions, that is, hard core interactions plus the mean field attraction “within” the small volume at $`𝐫`$. By taking the time derivative of (40) and using the hydrodynamic equations (31), a standard calculation leads to

$`_t\rho _𝐫`$ $`=`$ $`\rho _𝐫𝐯_𝐫,`$ (41)
$`_t𝐠_𝐫`$ $`=`$ $`𝐯_𝐫𝐠_𝐫P^{\mathrm{vdW}}(u_𝐫,\rho _𝐫)+c\rho _𝐫^2\rho _𝐫\tau ,`$ (42)
$`_tu_𝐫`$ $`=`$ $`u_𝐫𝐯_𝐫P^{\mathrm{vdW}}(u_𝐫,\rho _𝐫)𝐯_𝐫𝐪\tau :𝐯_𝐫,`$ (43)
where the heat flux is given by

$$𝐪_𝐫=\lambda _𝐫T^{\mathrm{vdW}}(u_𝐫,\rho _𝐫).$$
(44)
Here, $`P^{\mathrm{vdW}},T^{\mathrm{vdW}}`$ are the two equations of state for the van der Waals fluid (see Appendix VIII). In obtaining (43) and (44), use has been made of the identities (167) in the Appendix VIII.

The hydrodynamic equations (43) are equivalent to those proposed by Felderhof . However, in Felderhof’s theory, dissipation was neglected and, consistently, thermal fluctuations where not considered. Thus, equations (43) with (32) and (44) can be taken as the natural generalization of Felderhof theory.

It is not obvious at a first glance that the total momentum of the system is conserved by the term $`c\rho _𝐫^2\rho _𝐫`$ in (43). That this is the case can be seen by noting that this term can be cast into the form $`𝐏_𝐫`$, and the momentum equation has the form of a balance equation. The “surface tension” contribution to the stress tensor $`𝐏_𝐫`$ has the form

$$𝐏_𝐫=2c\left[\frac{1}{6}\rho _𝐫\rho _𝐫+\frac{1}{12}(\rho _𝐫)^2\mathrm{𝟏}\frac{1}{3}\rho _𝐫\rho _𝐫\frac{1}{6}\rho _𝐫^2\rho _𝐫\mathrm{𝟏}\right].$$
(45)
This contribution to the stress tensor is identical to the one derived from molecular considerations for the equilibrium stress tensor  under the gradient approximation. We observe, therefore, that the theory that is being presented relies also on the hypothesis of local equilibrium, in much the same way as in the usual continuum hydrodynamics approach.

## V Summary and discussion

We have computed the entropy functional for a simple fluid starting from microscopic principles under the classic van der Waals assumption that the molecular potential has two well-separated ranges, a hard core repulsive part plus a very long ranged attractive tail. Under the further assumption that the hydrodynamic fields vary slowly on the total range of the potential it is possible to write a local entropy functional that depends on the density gradients. We have also presented the connection between this entropy functional and the local gradient free energy functional which is used in the context of the study of equilibrium of liquid-vapor interfaces. Even though the assumption of the unrealistic molecular potential seems to be too restrictive, the equilibrium gradient theory is highly successful in describing the equilibrium properties of these interfaces, and gives confidence on the dynamical approach taken in this paper.

The entropy functional describes the equilibrium properties of the system but it is also one of the basic inputs for constructing the non-equilibrium evolution equations for the hydrodynamic fields through the GENERIC formalism. By making use of all the well-known information about the $`L`$ and $`M`$ matrices for the hydrodynamics of a single phase fluid flow, we can construct the hydrodynamic equations for a van der Waals fluid in a rather simple way.

The hydrodynamic equations proposed predict unstable forces if the thermodynamic state is in the positive slope part of the van der Waals pressure diagram. When this happens, the fluid spontaneously adopts one of the two possible densities corresponding to liquid and vapor. Density variations appear in the short length scales associated to the interfaces of bubbles or drops. It is precisely in these interfacial regions where the last term in Eqn. (39) is important. This term generates forces in the fluid associated to the surface tension of the interface. Actually, the interfaces try to adopt a spherical shape. The interface is a diffuse object with a finite width. Note that no boundary conditions are required on the interface of a bubble or a drop. Actually, there is no need to known a priori the location of the “interface” in order to have a well-posed hydrodynamic problem.

The set of equations (31) with the gradient approximation (14) and the set of equations (43) are mathematically equivalent. We regard as the main benefit of the microscopic derivation presented in this paper the clear physical interpretation of the variables (either $`ϵ_𝐫^i`$ or $`u_𝐫`$) and the correct equations of state (either $`P^{\mathrm{hc}}`$ or $`P^{\mathrm{vdW}}`$) to be used in the equations. This issue is not always clear in the usual presentations of this subject based on phenomenological arguments that simply add gradient terms either in the energy or the entropy .

From the computational point of view of numerically solving the hydrodynamic equations obtained, we observe that a theory like Kawasaki’s in which the mean field long range forces appear in integral form, Eqns. (31), might represent an advantage in front of a theory like Felderhof’s in which these forces appear in differential form through third derivatives of the density field, Eqns. (43). In those regions where the interfaces are located, the density field changes strongly which means that the third derivatives present a rather complex structure in very short length scales. In order to get stable and convergent results, very fine grids are required. These problems do not arise in the discretization of the integral form of the long range forces.

The theory derived here connects with the well-known phase field theories that deal with systems that present complex boundary conditions like those occurring in melting and dendritic growth  or immiscible fluids . In these systems, the partial differential equations governing the evolution of the relevant variables require the formulation of boundary conditions on moving surfaces. The motion of the surface is, in turn, coupled in a complex way with the dynamics of the relevant variables. The approach taken by the phase field method is to include the boundary conditions into the equations of motion through the introduction of an auxiliary phase field that locates in a somewhat diffuse way the position of the interface. A further dynamics is generated for the phase field which is coupled with the dynamics of the rest of fields in the system. This dynamics is obtained from a phenomenological formulation of a “free energy”. The theory we have derived is then in the class of a phase field theory for bubbles and drops. In our case, the phase field has a definite physical meaning in terms of the density field, and no phenomenological free-energy needs to be introduced ad hoc.

P.E. would like to thank useful discussions with H.C. Öttinger, P. Tarazona, J. Cuesta, and R. Delgado. This work has been partially supported by DGYCIT PB97-0077.

## VI Appendix: Microscopic derivation of the entropy functional

In this section, we compute from first principles the joint probability that an extended simple fluid in equilibrium has a particular realization of the mass, momentum and internal energy density fields. This generalizes van Kampen derivation which only considered the mass density field .

The isolated fluid is assumed to be described microscopically by $`z`$, the set of positions $`𝐪_i`$ and momenta $`𝐩_i`$ of the center of mass of its $`N_0`$ constituent molecules (no rotational or vibrational degrees of freedom are considered, for simplicity). The Hamiltonian of the system is given by

$$H(z)=\underset{i}{\overset{N_0}{}}\frac{𝐩_i^2}{2m}+\frac{1}{2}\underset{ij}{}\varphi _{ij},$$
(46)
where $`m`$ is the mass of a molecule and $`\varphi _{ij}`$ is the pair-wise potential function between molecules $`i`$ and $`j`$, assumed to depend only on its separation $`r_{ij}`$. The total momentum of the system is given by

$$𝐏(z)=\underset{i}{\overset{N_0}{}}𝐩_i.$$
(47)
The microscopic dynamics is such that $`H(z)`$ and $`𝐏(z)`$ are dynamical invariants. Obviously, the total number $`N_0`$ of particles (or the total mass) is also a dynamical invariant.

On a mesoscopic scale, the space containing the fluid is partitioned in $`M`$ cells. The state of the system at that scale is described with the mesoscopic variables $`N_\mu ,𝐏_\mu ,ϵ_\mu `$, $`\mu =1,\mathrm{},M`$, where $`N_\mu `$ is the number of particles in cell $`\mu `$, $`𝐏_\mu `$ is the momentum of the center of mass of the particles of cell $`\mu `$ and $`ϵ_\mu `$ is the energy with respect to the center of mass of the particles in cell $`\mu `$. In order to relate the mesoscopic variables with the microscopic ones, we introduce the characteristic function $`\chi _\mu (𝐫)`$ of cell $`\mu `$ which takes the value 1 if $`𝐫`$ is within cell $`\mu `$ and 0 otherwise. For example, the space can be partitioned into Voronoi cells whose characteristic function is $`\chi _\mu (𝐫)_{\nu \mu }^M\theta (|𝐫𝐑_\nu ||𝐫𝐑_\mu |)`$ where $`\theta (x)`$ is the Heaviside step function and $`𝐑_\mu `$ is the center of the Voronoi cell $`\mu `$. Obviously,

$$\underset{\mu }{\overset{M}{}}\chi _\mu (𝐫)=1,$$
(48)
because the cells cover all space without overlapping. The volume of the cell is given by

$$𝒱_\mu =_{V_0}d^3𝐫\chi _\mu (𝐫),$$
(49)
where $`V_0`$ is the total volume of the system. The mesoscopic variables $`N_\mu `$, $`𝐏_\mu `$ and $`ϵ_\mu `$ can now be written as functions of the microstate $`z`$, this is

$`N_\mu (z)`$ $`=`$ $`{\displaystyle \underset{i}{}}\chi _\mu (𝐪_i),`$ (50)
$`𝐏_\mu (z)`$ $`=`$ $`{\displaystyle \underset{i}{}}𝐩_i\chi _\mu (𝐪_i),`$ (51)
$`ϵ_\mu (z)`$ $`=`$ $`{\displaystyle \underset{i}{}}\left({\displaystyle \frac{1}{2m}}\left(𝐩_i{\displaystyle \frac{𝐏_\mu }{N_\mu }}\right)^2+\varphi _i\right)\chi _\mu (𝐪_i),`$ (52)
where we have introduced the potential energy of particle $`i`$ through

$$\varphi _i=\frac{1}{2}\underset{j}{}\varphi _{ij}.$$
(53)
Note that we use Latin indices to refer to variables at the microscopic level and Greek indices to refer to mesoscopic variables.

The dynamical invariants can be written in terms of the mesoscopic variables by using Eqn. (48), this is

$`H(z)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle \frac{𝐏_\mu ^2(z)}{2M_\mu (z)}}+ϵ_\mu (z),`$ (54)
$`𝐏(z)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}𝐏_\mu (z),`$ (55)
where $`M_\mu =mN_\mu `$ is the mass of cell $`\mu `$.

We will assume that the system is at equilibrium, which means that the effective probability density $`\rho ^{\mathrm{eq}}(z)`$ of finding a particular value of $`z`$ is a function of the dynamical invariants of the system only. If we assume that these invariants are known with precision, the equilibrium distribution function is given by the molecular ensemble

$$\rho ^{\mathrm{eq}}(z)=\frac{1}{N_0!\mathrm{\Omega }_0}\delta (𝐏(z)𝐏_0)\delta (H(z)E_0),$$
(56)
where $`𝐏_0`$ is the total momentum of the system and $`E_0`$ the total energy. The factor $`N_0!`$ is quantum mechanically in origin (it comes from the undistinguishability of the particles) and solves the Gibbs paradox (and makes the macroscopic entropy a first order function of energy and number of particles). The normalization factor $`\mathrm{\Omega }_0`$ is given by

$$\mathrm{\Omega }_0=\mathrm{\Omega }(U_0,N_0,V_0)=\frac{dz}{N_0!}\delta (𝐏(z)𝐏_0)\delta (H(z)E_0).$$
(57)
Here, $`U_0=E_0P_0^2/2M_0`$ is the internal energy of the total system and the total volume $`V_0`$ appears parametrically as the region of integration of positions. The macroscopic entropy for this system is defined as

$$S^{\mathrm{mac}}(U,N,V)=k_B\mathrm{ln}\mathrm{\Omega }_0(U,N,V).$$
(58)
The intensive parameters are defined as the derivatives of the macroscopic entropy, this is,

$`{\displaystyle \frac{1}{T}}`$ $`=`$ $`{\displaystyle \frac{S^{\mathrm{mac}}}{U}},`$ (59)
$`{\displaystyle \frac{\mu }{T}}`$ $`=`$ $`{\displaystyle \frac{S^{\mathrm{mac}}}{N}},`$ (60)
$`{\displaystyle \frac{P}{T}}`$ $`=`$ $`{\displaystyle \frac{S^{\mathrm{mac}}}{V}},`$ (61)
where $`T`$ is the temperature, $`\mu `$ the chemical potential, and $`P`$ the pressure.

### A Calculation of $`P[N,P,ϵ]`$

The probability that a set of functions $`X(z)`$ take particular values $`x`$ when the system is at equilibrium is given by

$$P(x)=𝑑z\delta (X(z)x)\rho ^{\mathrm{eq}}(z),$$
(62)
where $`\rho ^{\mathrm{eq}}(z)`$ is the equilibrium ensemble. We now consider the functions $`X(z)`$ that appear in Eqn. (62) to be the set of mesoscopic variables (52). In this way, the probability that the system adopts a particular set of values $`N_\mu ,𝐏_\mu ,ϵ_\mu `$ for each cell is given by

$`P[N,𝐏,ϵ]`$ $`=`$ $`{\displaystyle 𝑑z\rho ^{\mathrm{eq}}(z)\underset{\mu }{\overset{M}{}}\chi \left(N_\mu (z)N_\mu \right)}`$ (63)
$`\times `$ $`\delta \left(𝐏_\mu (z)𝐏_\mu \right)\delta \left(ϵ_\mu (z)ϵ_\mu \right).`$ (64)
In this expression, we have introduced the characteristic function $`\chi \left(N_\mu (z)N_\mu \right)`$ that takes the value $`1`$ (not infinity) for the region of phase space in which the microstate $`z`$ produces exactly the value $`N_\mu `$ for the number of particles in cell $`\mu `$. The probability (64) is normalized according to

$$𝑑𝐏_1\mathrm{}𝑑𝐏_M𝑑ϵ_1\mathrm{}𝑑ϵ_M\underset{N_1,\mathrm{},N_M}{\overset{N_0}{}}P[N,𝐏,ϵ]=1,$$
(65)
where the numbers $`N_\mu `$ are subject to $`_\mu N_\mu =N_0`$. By using the equilibrium ensemble (56) and Eqns. (55) we obtain

$`P[N,𝐏,ϵ]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }_0}}\delta \left({\displaystyle \underset{\mu }{\overset{M}{}}}𝐏_\mu 𝐏_0\right)`$ (66)
$`\times `$ $`\delta \left({\displaystyle \underset{\mu }{\overset{M}{}}}\left({\displaystyle \frac{𝐏_\mu ^2}{2M_\mu }}+ϵ_\mu \right)E_0\right)`$ (67)
$`\times `$ $`{\displaystyle \frac{dz}{N_0!}\underset{\mu }{\overset{M}{}}\chi (N_\mu (z)N_\mu )\delta (𝐏_\mu (z)𝐏_\mu )}`$ (68)
$`\times `$ $`\delta (ϵ_\mu (z)ϵ_\mu ).`$ (69)
The term $`_\mu ^M\chi (N_\mu (z)N_\mu )`$ within the integral is zero unless the microstate $`z`$ is such that there are exactly $`N_1`$ particles in the first cell, $`N_2`$ in the second, etc. There are $`N_0!/N_1!\mathrm{}N_M!`$ ways of having the $`N_0`$ particles distributed among the $`M`$ cells with the prescribed numbers in each cell. Therefore, we can write

$`{\displaystyle \frac{dz}{N_0!}\underset{\mu }{\overset{M}{}}\chi (N_\mu (z)N_\mu )\delta (𝐏_\mu (z)𝐏_\mu )\delta (ϵ_\mu (z)ϵ_\mu )}`$ (70)
$`={\displaystyle \frac{1}{N_1!\mathrm{}N_M!}}{\displaystyle _{𝒱_1}}\underset{N_1}{\underset{}{d𝐪\mathrm{}d𝐪}}\mathrm{}{\displaystyle _{𝒱_M}}\underset{N_M}{\underset{}{d𝐪\mathrm{}d𝐪}}`$ (71)
$`\times {\displaystyle }d𝐩_1\mathrm{}d𝐩_M{\displaystyle \underset{\mu }{\overset{M}{}}}\delta ({\displaystyle \underset{i_\mu }{\overset{N_\mu }{}}}𝐩_{i_\mu }𝐏_\mu )`$ (72)
$`\times \delta \left({\displaystyle \underset{i_\mu }{\overset{N_\mu }{}}}\left({\displaystyle \frac{1}{2m}}\left(𝐩_{i_\mu }{\displaystyle \frac{𝐏_\mu }{N_\mu }}\right)^2+\varphi _{i_\mu }\right)ϵ_\mu \right)`$ (73)
$`={\displaystyle \frac{1}{N_1!\mathrm{}N_M!}}{\displaystyle _{𝒱_1}}\underset{N_1}{\underset{}{d𝐪_1\mathrm{}d𝐪_{N_1}}}\mathrm{}{\displaystyle _{𝒱_M}}\underset{N_M}{\underset{}{d𝐪\mathrm{}d𝐪_M}}`$ (74)
$`\times {\displaystyle \underset{\mu }{\overset{M}{}}}\mathrm{\Phi }(\mathrm{𝟎},ϵ_\mu {\displaystyle \underset{i_\mu }{\overset{N_\mu }{}}}\varphi _{i_\mu },N_\mu ),`$ (75)
where the function $`\mathrm{\Phi }(𝐏,E,N)`$ is introduced and computed in the Appendix VII. The numbers $`N_1,\mathrm{},N_M`$ are subject to $`_\mu N_\mu =N_0`$ and, therefore, we can extend the above result to arbitrary $`N_1,\mathrm{},N_M`$ provided that we multiply Eqn. (75) with the factor $`\chi (_\mu N_\mu N_0)`$, which takes the value 1 if its argument is zero.

Let us write now the potential energy of cell $`\mu `$ as follows,

$$\underset{i_\mu }{\overset{N_\mu }{}}\varphi _{i_\mu }=\frac{1}{2}\underset{i_\mu }{\overset{N_\mu }{}}\underset{j_\mu }{\overset{N_\mu }{}}\varphi _{i_\mu j_\mu }+\frac{1}{2}\underset{\nu }{\overset{}{}}\underset{i_\mu }{\overset{N_\mu }{}}\underset{j_\nu }{\overset{N_\nu }{}}\varphi _{i_\mu j_\nu }.$$
(76)
The first term is the potential energy due to the particles of cell $`\mu `$ and it is a function solely of the coordinates of the particles that are in cell $`\mu `$. The second term, the potential energy due to the interaction of particles in different cells, involve coordinates of particles in cells $`\nu `$ different from $`\mu `$. It is this last term which hinders the decoupling of the integrals in positions in Eqn. (75). For this reason, it is necessary to make approximations to this last term. In what follows, we will approximate the potential energy due to different cells by a mean field constant (independent of coordinates of particles), this is,

$$\frac{1}{2}\underset{\nu }{\overset{}{}}\underset{i_\mu }{\overset{N_\mu }{}}\underset{j_\nu }{\overset{N_\nu }{}}\varphi _{i_\mu j_\nu }\underset{\nu }{\overset{}{}}\frac{1}{2}N_\mu N_\nu C_{\mu \nu }.$$
(77)
This approximation is fully justified when the potential $`\varphi (r)`$ has two typical ranges, a short length scale $`r_s`$ much smaller than the typical size of Voronoi cells and a long length scale $`r_l`$ much larger than the size of the cells. In this case, we can separate $`\varphi (r)=\varphi ^s(r)+\varphi ^l(r)`$ and write

$$\frac{1}{2}\underset{i_\mu }{\overset{N_\mu }{}}\underset{j_\nu }{\overset{N_\nu }{}}\varphi _{i_\mu j_\nu }\frac{1}{2}\underset{i_\mu }{\overset{N_\mu }{}}\underset{j_\nu }{\overset{N_\nu }{}}\varphi _{i_\mu j_\nu }^l,$$
(78)
because the short range part is negligible for different cells $`\mu \nu `$. Also, because the range $`r_l`$ is much larger than the size of a cell, we can further approximate Eqn. (78) as if all particles in each cell where located at the center of the cell and, therefore, we can write

$$\frac{1}{2}\underset{i_\mu }{\overset{N_\mu }{}}\underset{j_\nu }{\overset{N_\nu }{}}\varphi _{i_\mu j_\nu }^l\frac{1}{2}\varphi ^l(R_{\mu \nu })N_\mu N_\nu ,$$
(79)
and the potential energy due to different cells has the structure of Eqn. (77).

Under the mean field approximation (77), the different position integrals in (75) decouple, so that Eqn. (75) can be written as

$$\underset{\mu }{\overset{M}{}}\left[\frac{1}{N_\mu !}_{𝒱_\mu }𝑑𝐪_{1_\mu }\mathrm{}𝑑𝐪_{N_\mu }\mathrm{\Phi }(\mathrm{𝟎},ϵ_\mu \varphi _\mu ,N_\mu )\right],$$
(80)
where

$`\varphi _\mu `$ $`=`$ $`\varphi _\mu (𝐪_{1_\mu },\mathrm{},𝐪_{N_\mu })={\displaystyle \frac{1}{2}}{\displaystyle \underset{i_\mu }{\overset{N_\mu }{}}}{\displaystyle \underset{j_\mu }{\overset{N_\mu }{}}}\varphi _{i_\mu j_\mu }+\overline{\varphi }_\mu `$ (81)
$`\overline{\varphi }_\mu `$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\nu }{\overset{}{}}}\varphi ^l(R_{\mu \nu })N_\mu N_\nu ,`$ (82)
where $`\overline{\varphi }_\mu `$ is the mean field inter-cell potential energy of cell $`\mu `$.

We can write Eqn. (80) as

$$\underset{\mu }{\overset{M}{}}\mathrm{\Omega }^{\mathrm{hc}}(ϵ_\mu \overline{\varphi }_\mu ,N_\mu ,𝒱_\mu ),$$
(83)
where $`\mathrm{\Omega }^{\mathrm{hc}}(U,N,V)`$ is defined in Eqn. (57). Note that the Hamiltonian that should appear in Eqn. (57) is that of a system of particles interacting through the short range part of the potential $`\varphi ^s(r)`$. Finally, Eqn. (69) can be written as

$`P[N,𝐏,ϵ]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }_0}}\chi \left({\displaystyle \underset{\mu }{\overset{M}{}}}N_\mu N_0\right)\delta \left({\displaystyle \underset{\mu }{\overset{M}{}}}𝐏_\mu 𝐏_0\right)`$ (84)
$`\times `$ $`\delta \left({\displaystyle \underset{\mu }{\overset{M}{}}}\left({\displaystyle \frac{𝐏_\mu ^2}{2M_\mu }}+ϵ_\mu \right)E_0\right)`$ (85)
$`\times `$ $`\mathrm{exp}\left\{S[N,ϵ]/k_B\right\},`$ (86)
where the entropy functional at the mesoscopic level is defined by

$$S[N,ϵ]=\underset{\mu }{}S^{\mathrm{hc}}(ϵ_\mu \overline{\varphi }_\mu ,N_\mu ,𝒱_\mu ),$$
(87)
where we have used Eqn. (58). The mesoscopic entropy is given in terms of the sum of the macroscopic entropies of each cell, as if they were isolated systems in a volume $`𝒱_\mu `$, with number of particles $`N_\mu `$ and with an energy $`ϵ_\mu \overline{\varphi }_\mu `$. The internal energy $`ϵ_\mu `$ in (52) includes in its microscopic definition the potential energy of particles that interact with particles of neighboring cells, that is, $`ϵ_\mu `$ includes the inter-cell energy. In this way, $`ϵ_\mu \overline{\varphi }_\mu `$ represents the purely internal energy due to the particles within the cell (in mean field). Even though we have made progress in writing the probability for the mesoscopic variables, Eqn. (86), the entropy of each cell $`S^{\mathrm{hc}}(ϵ_\mu ,𝒱_\mu ,N_\mu )`$ is still an unknown quantity. We remark that this entropy is the macroscopic entropy of a system of molecules interacting with the short range part of the potential. In applications of the model, this entropy will be given by simple models (see Appendix VIII) or by very accurate expressions like the Carnahan-Starling equation of state .

Apparently, Eqn. (87) simply says that “the entropy is additive”. However, this sentence is imprecise: The entropy $`S[N,ϵ]`$ in the lhs of Eqn. (87) is a different object from the entropy $`S^{\mathrm{hc}}(ϵ_\mu ,N_\mu ,𝒱_\mu )`$ appearing in the rhs. They depend on a different number of variables. Therefore, one should rather say that “the mesoscopic entropy at the level of hydrodynamic variables is the sum of the entropy at the level of dynamical invariants (the macroscopic entropy) of each cell as if they were isolated”.

The macroscopic entropy is a first order function of its variables and therefore,

$$S[N,ϵ]=\underset{\mu }{}𝒱_\mu S^{\mathrm{hc}}(\frac{ϵ_\mu }{𝒱_\mu }\frac{\overline{\varphi }_\mu }{𝒱_\mu },\frac{N_\mu }{𝒱_\mu },1)$$
(88)
which admits the continuum notation

$$S[n_𝐫,ϵ_𝐫]=𝑑𝐫s^{\mathrm{hc}}(ϵ_𝐫\overline{\varphi }_𝐫,n_𝐫),$$
(89)
where we have introduced the continuum version of the effective interaction energy between molecules in different points of space (cells),

$`𝐑_\mu `$ $``$ $`𝐫,`$ (90)
$`{\displaystyle \underset{\mu }{}}𝒱_\mu `$ $``$ $`{\displaystyle 𝑑𝐫},`$ (91)
$`{\displaystyle \frac{ϵ_\mu }{𝒱_\mu }}`$ $``$ $`ϵ_𝐫,`$ (92)
$`{\displaystyle \frac{N_\mu }{𝒱_\mu }}`$ $``$ $`n_𝐫,`$ (93)
$`{\displaystyle \frac{\overline{\varphi }_\mu }{𝒱_\mu }}`$ $``$ $`\overline{\varphi }_𝐫.`$ (94)
Note that in this continuum notation we can write the dynamical invariants (55) of the system as

$`E`$ $`=`$ $`{\displaystyle 𝑑𝐫\left(\frac{1}{2}\frac{𝐠_𝐫}{mn_𝐫}+ϵ_𝐫\right)},`$ (95)
$`P`$ $`=`$ $`{\displaystyle 𝑑\mathrm{𝐫𝐠}_𝐫},`$ (96)
where $`𝐏_\mu /𝒱_\mu 𝐠_𝐫`$ is the momentum density field. We stress that we regard Eqn. (89) not as an strict mathematical limit in which the volume of the cells tend to zero, but rather as a convenient notational tool.

In spite of its apparent form, the mesoscopic entropy (89) is not local in space. This is, it is not given by a sum of a function of the variables at a given point. This is due to the presence of the term $`\overline{\varphi }_𝐫`$ which can be written in continuous notation as

$$\overline{\varphi }_𝐫=𝑑𝐫^{}\varphi ^l(|𝐫𝐫^{}|)n_𝐫n_𝐫^{}.$$
(97)
We have implicitly associated the potential energy of interaction with other cells $`\overline{\varphi }_𝐫`$ with the phenomena of surface tension. In order to make this connection more explicit, in the following subsections we consider the marginal probability distribution functions and will make contact with the Density Functional Theory.

### B The marginal distribution $`P[N,ϵ]`$

By integrating the distribution function $`P[N,𝐏,ϵ]`$ over momenta we will have the probability of a realization of the “fields” $`N,ϵ`$ irrespective of the values of the momenta in each cell, this is

$`P[N,ϵ]`$ $`=`$ $`\mathrm{exp}\left\{S[N,ϵ]/k_B\right\}{\displaystyle \frac{1}{\mathrm{\Omega }_0}}\chi \left({\displaystyle \underset{\mu }{\overset{M}{}}}N_\mu N_0\right)`$ (98)
$`\times `$ $`{\displaystyle 𝑑𝐏_1\mathrm{}𝑑𝐏_M\delta \left(\underset{\mu }{\overset{M}{}}𝐏_\mu 𝐏_0\right)}`$ (99)
$`\times `$ $`\delta \left({\displaystyle \underset{\mu }{\overset{M}{}}}\left({\displaystyle \frac{𝐏_\mu ^2}{2M_\mu }}+ϵ_\mu \right)E_0\right)`$ (100)
$`=`$ $`\mathrm{exp}\left\{S[N,ϵ]/k_B\right\}{\displaystyle \frac{1}{\mathrm{\Omega }_0}}\chi \left({\displaystyle \underset{\mu }{\overset{M}{}}}N_\mu N_0\right)`$ (101)
$`\times `$ $`{\displaystyle \underset{\mu }{\overset{M}{}}}(2mN_\mu )^{D/2}{\displaystyle \frac{\omega _{D(M1)}}{2}}\left[U_0{\displaystyle \underset{\mu }{}}ϵ_\mu \right]^{\frac{D(M1)}{2}1},`$ (102)
where we have used once more Eqn. (150) in the appendix.

The most probable realization $`N_1^{},\mathrm{},N_M^{},ϵ_1^{},\mathrm{},ϵ_M^{}`$ of the fields is the one which maximizes the functional

$`k_B^1S[N,ϵ]+\left({\displaystyle \frac{D(M1)}{2}}1\right)\mathrm{log}\left(U_0{\displaystyle \underset{\mu }{}}ϵ_\mu \right)`$ (104)
$`+{\displaystyle \frac{D}{2}}{\displaystyle \underset{\mu }{}}\mathrm{log}(2mN_\mu )+\beta \lambda {\displaystyle \underset{\mu }{}}N_\mu ,`$ (105)
where we have introduce the Lagrange multiplier $`\beta \lambda `$ that takes into account the restriction $`_\mu N_\mu =N_0`$ The maximum $`N^{},ϵ^{}`$ of the functional (105) is the solution of the following set of equations

$`{\displaystyle \frac{S}{N_\mu }}[N^{},ϵ^{}]`$ $`=`$ $`k_B\beta \lambda ,`$ (106)
$`{\displaystyle \frac{S}{ϵ_\mu }}[N^{},ϵ^{}]`$ $`=`$ $`{\displaystyle \frac{(\frac{D(M1)}{2}1)k_B}{U_0_\nu ϵ_\nu ^{}}},`$ (107)
$`{\displaystyle \underset{\mu }{}}N_\nu ^{}`$ $`=`$ $`N_0.`$ (108)
We have neglected a term of order $`N_\mu ^1`$ in the first equation, which is reasonable if typically there are many particles in each cell. Eqns. (108) is a set of $`2M+1`$ equations for the unknowns $`N_\mu ^{},ϵ_\mu ^{},\lambda `$. Note that the solution $`N_\mu ^{},ϵ_\mu ^{},\beta \lambda `$ depends parametrically on $`N_0,U_0`$.

We find now a convenient approximation to Eqn. (LABEL:pne) by noting that this probability is expected to be highly peaked around the most probable state. Therefore, for those values of the field $`ϵ`$ for which $`P[N,ϵ]`$ is appreciably different from zero (that is, around $`ϵ^{}`$), we can approximate

$`\left[U_0{\displaystyle \underset{\mu }{}}ϵ_\mu \right]^P=\left[U_0{\displaystyle \underset{\mu }{}}ϵ_\mu ^{}\right]^P`$ (109)
$`\times `$ $`\left[1+{\displaystyle \frac{\beta }{P}}{\displaystyle \underset{\mu }{}}(ϵ_\mu ^{}ϵ_\mu )\right]^P`$ (110)
$``$ $`\left[U_0{\displaystyle \underset{\mu }{}}ϵ_\mu ^{}\right]^P\mathrm{exp}\{\beta {\displaystyle \underset{\mu }{}}\left(ϵ_\mu ^{}ϵ_\mu \right)\}`$ (111)
$`=`$ $`\mathrm{ctn}.\mathrm{exp}\{\beta {\displaystyle \underset{\mu }{}}ϵ_\mu \},`$ (112)
where $`P=D(M1)/21`$ is a very large number and we have introduced

$$\beta =\frac{D(M1)/21}{U_0_\nu ϵ_\nu ^{}}.$$
(113)
Note that $`\beta `$ is proportional to the inverse of the most probable value of the kinetic energy per cell.

Finally, we can write Eqn. (LABEL:pne) as

$`P[N,ϵ]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }_0^{}}}\mathrm{exp}\left\{S[N,ϵ]/k_B\beta {\displaystyle \underset{\mu }{}}ϵ_\mu \right\}`$ (114)
$`\times `$ $`\chi \left({\displaystyle \underset{\mu }{}}N_\mu N_0\right),`$ (115)
where $`\mathrm{\Omega }_0^{}`$ is the corresponding normalization function. We see, therefore, that by integrating the momenta the “microcanonical” form Eqn. (86) becomes the “canonical” form (115).

### C The marginal distribution $`P[N]`$

There are two different routes to compute the probability of a certain distribution $`N_1,\mathrm{},N_M`$ of the particles in the cells. We could simply start the calculation from Eqn. (64) without the momentum and energy conserving delta functions. This route is essentially the one taken by van Kampen . The alternative is to integrate out the energy field in Eqn. (115). The same results are obtained in both approaches and we illustrate here this second one, this is,

$`P[N]`$ $`=`$ $`{\displaystyle 𝑑ϵ_1\mathrm{}𝑑ϵ_MP[N,ϵ]}`$ (116)
$`=`$ $`{\displaystyle \frac{1}{Z_0}}\mathrm{exp}\left\{\beta {\displaystyle \underset{\mu }{}}F^{\mathrm{hc}}(\beta ,N_\mu ,𝒱_\mu )\beta \overline{\varphi }_\mu \right\}`$ (117)
$`\times `$ $`\chi \left({\displaystyle \underset{\mu }{\overset{M}{}}}N_\mu N_0\right),`$ (118)
where we have used the form (87) for the entropy function, we have changed variables to $`ϵ_\mu ^{}=ϵ_\mu \overline{\varphi }_\mu `$ and have introduced the macroscopic free energy $`F^{\mathrm{hc}}(\beta ,N,V)`$ of a system of $`N`$ particles at inverse temperature $`\beta `$ in a volume V, and interacting with a short range potential. In deriving Eqn. (118) we have made use of the following identity

$$_0^{\mathrm{}}\mathrm{exp}\{S(U)/k_B\beta U\}𝑑U=\mathrm{exp}\{\beta F(\beta )\}.$$
(119)
This expression is the statistical mechanics link between the macroscopic entropy and the free energy and can be proved by computing the Laplace transform of $`\mathrm{\Omega }_0`$ in Eqn. (57). One has

$`{\displaystyle _0^{\mathrm{}}}\mathrm{\Omega }_0(U,N,V)\mathrm{exp}\{\beta U\}dU`$ (120)
$`={\displaystyle 𝑑z\delta (\underset{i}{}𝐩_i)\mathrm{exp}}\{\beta H(z)\}=Z(\beta ,N1,V),`$ (121)
where the partition function $`Z(\beta )`$ is defined through this equation. The presence or absence of the momentum conservation delta function is irrelevant when the number of particles is large so we may very well drop it. The free energy $`F(\beta )`$ is defined as

$$\beta F(\beta )=\mathrm{ln}Z(\beta ),$$
(122)
The usual thermodynamic link between entropy and free energy,

$$F(T,N,V)=UTS(U,N,V),$$
(123)
can be obtained under the assumption that the integrand in (119) is highly peaked.

## VII Appendix: Molecular ensemble

In this appendix we compute explicitly the following integral

$`\mathrm{\Phi }(𝐏_0,E_0,M)`$ $`=`$ $`{\displaystyle d^{DM}𝐩\delta \left(\underset{i}{\overset{M}{}}\frac{𝐩_i^2}{2m_i}E_0\right)}`$ (124)
$`\times `$ $`\delta ^D\left({\displaystyle \underset{i}{\overset{M}{}}}𝐩_i𝐏_0\right),`$ (125)
which appears repeatedly when computing molecular averages. By a simple change of variables we obtain

$`\mathrm{\Phi }(𝐏_0,E_0,M)`$ $`=`$ $`{\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{D/2}{\displaystyle d^{DM}𝐩\delta \left(\underset{i}{\overset{M}{}}𝐩_i^2E_0\right)}`$ (126)
$`\times `$ $`\delta ^D\left({\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{1/2}𝐩_i𝐏_0\right),`$ (127)
The equation $`_i^M(2m_i)^{1/2}𝐩_i=𝐏_0`$ are actually D equations (one for each component of the momentum) which define D planes in $`R^{DM}`$. The integral in (125) is actually over a submanifold which is the intersection of the D planes with the surface of a $`DM`$ dimensional sphere of radius $`E_0^{1/2}`$. This intersection will be also a sphere, which will be now of smaller radius and also of smaller dimension, $`D(M1)`$.

In order to compute (125), we change to the following notation

$`𝓟`$ $`=`$ $`(p_1^x,\mathrm{},p_M^x,p_1^y,\mathrm{},p_M^y,p_1^z,\mathrm{},p_M^z)`$ (128)
$`𝓒_x`$ $`=`$ $`((2m_1)^{1/2},\mathrm{},(2m_M)^{1/2},0,\mathrm{},0,0,\mathrm{},0)`$ (129)
$`𝓒_y`$ $`=`$ $`(0,\mathrm{},0,(2m_1)^{1/2},\mathrm{},(2m_M)^{1/2},0,\mathrm{},0)`$ (130)
$`𝓒_z`$ $`=`$ $`(0,\mathrm{},0,0,\mathrm{},0,(2m_1)^{1/2},\mathrm{},(2m_M)^{1/2}).`$ (131)
Note that these vectors satisfy $`𝓒_x𝓒_y=0`$ $`𝓒_y𝓒_z=0`$ $`𝓒_z𝓒_x=0`$. With these vectors so defined, Eqn. (125) becomes

$`\mathrm{\Phi }(𝐏_0,E_0,M)`$ $`=`$ $`{\displaystyle d^{DM}𝓟\delta \left(𝓟^2E_0\right)\delta \left(𝓒_x𝓟P_0^x\right)}`$ (132)
$`\times `$ $`\delta \left(𝓒_y𝓟P_0^y\right)\delta \left(𝓒_z𝓟P_0^z\right)`$ (133)
$`\times `$ $`{\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{D/2}.`$ (134)
Now we consider the following change of variables

$$𝓟^{}=𝓟\left(\frac{𝓒_x}{|𝓒_x|^2}P_0^x+\frac{𝓒_y}{|𝓒_y|^2}P_0^y+\frac{𝓒_z}{|𝓒_z|^2}P_0^z\right),$$
(135)
which is simply a translation and has unit Jacobian. Simple algebra leads to

$`\mathrm{\Phi }(𝐏_0,E_0,M)`$ $`=`$ $`{\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{D/2}{\displaystyle d^{DM}𝓟^{}\delta \left(𝓟\text{}^2U_0\right)}`$ (136)
$`\times `$ $`\delta \left(𝓒_x𝓟^{}\right)\delta \left(𝓒_y𝓟^{}\right)\delta \left(𝓒_z𝓟^{}\right),`$ (137)
where we have introduced the total internal energy

$$U_0=\left(E_0\frac{𝐏_0^2}{2_0}\right),$$
(138)
where $`_0=_im_i`$ is the total mass.

We now consider a second change of variables $`𝓟^{\prime \prime }=𝚲𝓟^{}`$ through a rotation $`𝚲`$ such that

$`𝚲𝓒_x`$ $`=`$ $`(1,\mathrm{},0,0,\mathrm{},0,0,\mathrm{},0)`$ (139)
$`𝚲𝓒_y`$ $`=`$ $`(0,\mathrm{},0,1,\mathrm{},0,0,\mathrm{},0)`$ (140)
$`𝚲𝓒_z`$ $`=`$ $`(0,\mathrm{},0,0,\mathrm{},0,1,\mathrm{},0).`$ (141)
It is always possible to find a matrix $`𝚲`$ that satisfies Eqns. (141). For example, consider a block diagonal matrix made of three identical blocks of size $`M\times M`$. Then assume that each block is the same orthogonal matrix which transforms the vector $`((2m_1)^{1/2},(2m_2)^{1/2},\mathrm{},(2m_M)^{1/2})`$ into $`(2_0)^{1/2}(1,0,\mathrm{},0)`$. After the rotation (which has unit Jacobian and leaves the modulus of a vector invariant) the integral (137) becomes

$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{D/2}{\displaystyle d^{DM}𝓟^{\prime \prime }\delta \left(𝓟^{\prime \prime 2}U_0\right)}`$ (142)
$`\times `$ $`\delta (p_1^{\prime \prime x})\delta (p_1^{\prime \prime y})\delta (p_1^{\prime \prime z})`$ (143)
$`=`$ $`{\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{D/2}{\displaystyle d^{D(M1)}𝓟^{\prime \prime }\delta \left(𝓟^{\prime \prime 2}U_0\right)}.`$ (144)
We compute now the integral over the sphere in Eqn. (144) by using that the integral of an arbitrary function $`F(𝐱)=f(|𝐱|)`$ that depends on $`𝐱`$ only through its modulus $`|𝐱|`$ can be computed by changing to polar coordinates

$$F(𝐱)d^M𝐱=\omega _M_0^{\mathrm{}}f(r)r^{M1}𝑑r.$$
(145)
The numerical factor $`\omega _M`$, which comes from the integration of the angles, can be computed by considering the special case when $`f(r)`$ is a Gaussian. The result is

$$\omega _M=2\frac{\pi ^{M/2}}{\mathrm{\Gamma }(M/2)}$$
(146)
By using Eqn. (145), Eqn. (144) becomes

$`\mathrm{\Phi }(𝐏_0,E_0,M)`$ $`=`$ $`{\displaystyle \underset{i}{\overset{M}{}}}(2m_i)^{D/2}\omega _{D(M1)}`$ (147)
$`\times `$ $`{\displaystyle 𝑑pp^{D(M1)1}\delta \left(p^2U\right)}.`$ (148)
We need now the property

$$\delta (f(x))=\underset{i}{}\frac{\delta (xx_i)}{|f^{}(x_i)|}$$
(149)
where $`x_i`$ are the zeros of $`f(x_i)=0`$. For the case of Eqn. (148) we have $`f(x)=x^2U`$, $`x_i=\pm (U)^{1/2}`$ and $`f^{}(x)=2x`$. Therefore,

$$\mathrm{\Phi }(𝐏_0,E_0,M)=\frac{1}{2}\omega _{D(M1)}U_0^{\frac{D(M1)2}{2}}\underset{i}{\overset{M}{}}(2m_i)^{D/2}.$$
(150)
## VIII Appendix: van der Waals and hard core models

The particular functional form of $`s^{\mathrm{hc}}(ϵ,\rho )`$ cannot be computed from the microscopic analysis presented in the previous appendix. The approach taken in this paper is that this fundamental thermodynamic equation (in the sense of Callen ) is known either from empirical sources or by suitable modelling. In this appendix we summarize the results for the fundamental equation of van der Waals model and for its corresponding hard core model.

The van der Waals model can be defined through the following fundamental equation, which relates the entropy density with the internal energy density $`ϵ`$ and number density $`n`$,

$$s^{\mathrm{vdW}}(ϵ,n)=\frac{D+2}{2}k_Bnk_Bn\mathrm{ln}\left(\frac{\left[\mathrm{\Lambda }^{\mathrm{vdW}}(ϵ,n)\right]^Dn}{1nb}\right).$$
(151)
The number density is $`n=\rho /m_0`$ where $`\rho `$ is the mass density and $`m_0`$ is the mass of a molecule. The thermal wavelength is defined by

$$\mathrm{\Lambda }^{\mathrm{vdW}}(ϵ,n)=\frac{h}{(2\pi m_0k_B𝒯^{\mathrm{vdW}}(ϵ,n))^{1/2}},$$
(152)
and we have introduced the function

$$𝒯^{\mathrm{vdW}}(ϵ,n)=\frac{2}{D}\frac{ϵ+an^2}{k_Bn}.$$
(153)
The two equations of state are the derivatives of the entropy with respect to each variable

$`{\displaystyle \frac{s^{\mathrm{vdW}}}{ϵ}}`$ $`=`$ $`{\displaystyle \frac{1}{T^{\mathrm{vdW}}(ϵ,n)}},`$ (154)
$`{\displaystyle \frac{s^{\mathrm{vdW}}}{n}}`$ $`=`$ $`{\displaystyle \frac{\mu ^{\mathrm{vdW}}(ϵ,n)}{T^{\mathrm{vdW}}(ϵ,n)}},`$ (155)
where $`T^{\mathrm{vdW}}`$ is the temperature and $`\mu ^{\mathrm{vdW}}`$ is the chemical potential. Straightforward calculations lead to

$`T^{\mathrm{vdW}}(ϵ,n)`$ $`=`$ $`𝒯^{\mathrm{vdW}}(ϵ,n),`$ (156)
$`\mu ^{\mathrm{vdW}}(ϵ,n)`$ $`=`$ $`k_BT^{\mathrm{vdW}}\left(\mathrm{ln}\left({\displaystyle \frac{n[\mathrm{\Lambda }^{\mathrm{vdW}}]^D}{1nb}}\right)+{\displaystyle \frac{nb}{1nb}}\right)2an.`$ (157)
The third equation of state can be obtained from the Euler equation $`P=Tsϵ+\mu n`$, which leads to

$$P^{\mathrm{vdW}}(ϵ,n)=n\frac{k_BT^{\mathrm{vdW}}(ϵ,n)}{1nb}an^2.$$
(158)
The hard core model that corresponds to the van der Waals model can be defined as the model that results from taking the attractive parameter $`a=0`$ in the van der Waals model. This model is essentially the gas ideal model but with excluded volume effects. In this way, one obtains

$$s^{\mathrm{hc}}(ϵ,n)=\frac{D+2}{2}k_Bnk_Bn\mathrm{ln}\left(\frac{\mathrm{\Lambda }^{\mathrm{hc}D}(ϵ,n)n}{1nb}\right),$$
(159)
where now

$`\mathrm{\Lambda }^{\mathrm{hc}}(ϵ,n)`$ $`=`$ $`{\displaystyle \frac{h}{(2\pi m_0k_B𝒯^{\mathrm{hc}}(ϵ,n))^{1/2}}},`$ (160)
$`𝒯^{\mathrm{hc}}(ϵ,n)`$ $`=`$ $`{\displaystyle \frac{2}{D}}{\displaystyle \frac{ϵ}{k_Bn}}.`$ (161)
The equations of state are now

$`T^{\mathrm{hc}}(ϵ,n)`$ $`=`$ $`𝒯^{\mathrm{hc}}(ϵ,n),`$ (162)
$`\mu ^{\mathrm{hc}}(ϵ,n)`$ $`=`$ $`k_BT^{\mathrm{vdW}}\left(\mathrm{ln}\left({\displaystyle \frac{n[\mathrm{\Lambda }^{\mathrm{hc}}]^D}{1nb}}\right)+{\displaystyle \frac{nb}{1nb}}\right),`$ (163)
$`P^{\mathrm{hc}}(ϵ,n)`$ $`=`$ $`n{\displaystyle \frac{k_BT^{\mathrm{hc}}(ϵ,n)}{1nb}}.`$ (164)
Note that the following functional relations hold

$`T^{\mathrm{hc}}(u+an^2,n)`$ $`=`$ $`T^{\mathrm{vdW}}(u,n),`$ (165)
$`P^{\mathrm{hc}}(u+an^2,n)an^2`$ $`=`$ $`P^{\mathrm{vdW}}(u,n),`$ (166)
$`\mu ^{\mathrm{hc}}(u+an^2,n)2an`$ $`=`$ $`\mu ^{\mathrm{vdW}}(u,n).`$ (167)


text:
# Order and nFl Behavior in UCu4Pd

## I Introduction

Since the first observation of non–Fermi liquid (nFl) behavior in 1991 , many intermetallic compounds have been found with low temperature thermodynamic properties that cannot be described within Landau’s Fermi liquid theory ($`\gamma =const.,\chi =const.,\rho \rho _rT^2`$). Typically, such alloys show a logarithmic or a weak power law dependence in C/T and in the magnetic susceptibility $`\chi `$ at low temperatures, as well as a non–T2 dependence of the electrical resistivity $`\rho `$ (for example $`\rho \rho _r`$ T). All of these compounds are close to magnetism. In searching for further examples for nFl–behavior, Andraka and Stewart  characterized UCu4Pd as a non-Fermi liquid. They found a power–law temperature dependence of the Sommerfeld coefficient $`\gamma =C/T=\mathrm{\Delta }T^\delta `$ ($`\mathrm{\Delta }`$ = 450 $`\frac{mJ}{moleK^{2\delta }}`$, $`\delta `$ = 0.32) over a decade in temperature between 1 and 10 K: the magnetic susceptibility $`\chi `$ also followed a power law with the same $`\delta `$ between 1.8 K and 10 K. Furthermore, the electrical resistivity increased linearly below 10 K with a slope six times larger than in U0.2Y0.8Pd3 . A possible explanation for this nFl–behavior is, according to ref. , the proximity to antiferromagnetism suppressed to TN = 0 in this compound. This is called the ”quantum critical point” (QCP) scenario.

Following this seminal work, much research has been done on UCu4Pd and its physical properties have been explained by various other theories. Castro Neto et al.  suggested that a Griffiths phase causes the nFl–behavior, whereupon de Andrade et al.  fit their UCu4Pd–data down to 0.6 K to this theory. Scheidt et al. , however, showed that their C/T data down to 60 mK matches this model only in a small temperature range around 1 K. This was also confirmed by Vollmer et al. .

Bernal et al.  interpreted the nFl–behavior of UCu4Pd within the framework of the Kondo disorder model. This model presupposes a distribution of microscopic exchange couplings and was argued to successfully explain the observed nFl–behavior in the magnetic susceptibility and the specific heat.

UCu4Pd is particularly well–suited for studying the influence of disorder because of the possibility of perfect chemical ordering in this compound which crystallizes in the cubic AuBe5 structure. In this ideal case the Pd atoms would occupy the 4c sites and the Cu atoms would occupy the smaller 16e sites. XAFS experiments , however, indicate that in an unannealed sample of UCu4Pd only 76 % of the Pd atoms populate the 4c sites, the remaining 24 % are on the 16e sites. Another hint for disorder in this compound is, according to ref. , that the change of the slope of the lattice parameter of UCu5-xPdx versus the Pd concentration x occurs at x $``$ 0.85 and not as would be expected at x = 1 for an ideal, ordered compound.

In this work we study annealed UCu4Pd–samples, where the Pd atoms are rearranged from the 16e sites to the 4c sites. This allows us to study UCu4Pd with less (but not zero) intrinsic disorder in order to analyze the importance of disorder for its nFl–behavior. The original discovery work  reported the effect of annealing on UCu3.5Pd1.5, but not on UCu4Pd. The only other annealing study of which we are aware reported that , within their error bar, no difference in disorder in ball milled (which introduces strain that broadens diffraction peak line widths) unannealed and annealed (at 925 C for one week) samples of UCu4Pd could be found in elastic neutron diffraction measurements.

We show by X-ray, ac-susceptibility, resistivity and specific heat data that the physical properties of UCu4Pd are extremely sensitive to annealing. We conclude that disorder may not be the underlying mechanism for the nFl–behavior in the specific heat and resistivity.

## II Sample Preparation

The samples of UCu4Pd presented in this work were all obtained from the same batch. This batch was arc–melted under a highly purified argon atmosphere. To obtain the highest possible homogeneity, the sample was flipped over four times and remelted. The loss in weight was 0.2 %. The sample was then cut in several pieces. These pieces were treated in different ways: one was measured as cast, one was annealed in an evacuated quartz glass tube for seven days at 750 C, one for 14 days at 750 C and another one was splat cooled (cooling rate more than 104 K/s) as a way to introduce more disorder than the as cast sample. All these samples were identified as single-phase (AuBe5 structure) in extra long counting time powder diffraction measurements. Also in specific heat and ac–susceptibility measurements no impurity phases like UCu5 or UPd3 were detected within $`1\%`$ accuracy. The annealed samples showed no weight loss within $`0.1\%`$ accuracy.

## III Results and Discussion

In the inset of figure 1) the lattice parameter a of unannealed UCu5-xPdx is plotted versus the Pd concentration x. The concentration dependence of the lattice parameter a can be described by two straight lines with an intersection at x $``$ 0.85. This points to the beginning of augmented (i.e. above random disorder present for x $`<`$ 0.85) occupation of the smaller 16e sites by the Pd atoms at this x = 0.85 concentration . The expanded plot around x = 1 (Figure 1a) shows in addition the lattice parameters of the two annealed and the one splat cooled samples. One observes that the lattice parameters of the annealed samples nearly lie on the straight line which fits the data points from unannealed samples with Pd concentrations lower than x = 0.8. The decrease of the lattice parameter with annealing can be attributed to a thermally activated movement of the Pd atoms from the smaller 16e sites to the bigger 4c sites. The lattice parameter of the splat cooled sample, however, is clearly much larger than the unannealed one, which is indicative of increased disorder. The lattice parameter of the splat cooled UCu4Pd sample is almost the same as that of unannealed UCu3.8Pd1.2.

A further test for increasing order upon annealing is the full width at half maximum of an X–ray peak. As an example, figure 1b) shows the \[3 1 1\] peak from X–ray powder diffraction measurements of the splat cooled, the unannealed, the 7 and the 14 days annealed UCu4Pd samples. The full width at half maximum clearly gets smaller in this sequence, which is a strong indication of increasing order with annealing the samples. The resistivity ratio of the 7 days and 14 days annealed samples between T=1.8 K and T=400 K decreases by $`11\%`$ (see figure 3), which also supports increasing order for extended annealing.

With these samples with different degrees of disorder, we are now able to make a systematic study of the role that disorder plays for the physical behavior of UCu4Pd.

In figure 2a) $`\mathrm{\Delta }`$C/T is displayed versus log T for all these samples except the splat cooled one. It can be seen that the antiferromagnetic transition  of the unannealed UCu4Pd sample (T$`{}_{N}{}^{}`$ 170 mK) is suppressed by annealing. The 7 days annealed sample of UCu4Pd still shows an inkling of such a transition ($``$ 140 mK), whereas UCu4Pd ”14 days annealed” follows a logarithmic dependence down to the lowest measured temperature of 0.08 K. This logarithmic temperature dependence expands upon annealing to cover more than two decades as can be seen in the inset of figure 2a).

To compare our specific heat data to the Griffiths phase model  that predicts a power–law behavior $`C/TT^{1+\lambda }`$ with $`\lambda <1`$, we have performed a careful $`\chi ^2`$–analysis of our data: i) A two–parameter logarithmic fit $`a_1+a_2\mathrm{ln}T`$ leads to $`\chi ^2=9.0`$ and $`\chi ^2=6.1`$ for fits in the temperature regions $`T<10K`$ and $`T<1K`$, resp. ii) A three–parameter power–law fit $`a_1^{}+a_2^{}T^{1+\lambda }`$ converges to the same $`\chi ^2`$–values with the fit parameter $`\lambda 1.0`$ and $`a_1^{}`$, $`a_2^{}`$ large and with opposite signs . One easily sees that this is effectively equivalent to our two–parameter logarithmic fit i). From this statistical analysis we conclude that our specific heat data below 10 K has a clear logarithmic temperature behavior.

It is remarkable that increasing order in UCu4Pd does not destroy its nFl–behavior in the specific heat, but rather expands the temperature range of its logarithmic behavior. Thus this work reports a new tuning parameter for nFl–behavior besides pressure , doping  or magnetic field : crystallographic order.

The splat cooled UCu4Pd, as expected for the most disordered sample, shows however a completely different behavior. The temperature dependence is more like that of UCu3.8Pd1.2 (which has almost the same lattice parameter as our ”splat cooled” sample, see fig. 1) as can be seen in figure 2b).

The electrical resistivity shows a particularly strong sensitivity to annealing. Figure 3 shows the electrical resistivity of the unannealed and the annealed samples. The residual resistivity drops strongly by a factor of about 2.5. Also an enormous change in the qualitative behavior of the curve can be observed. While the resistivity of the unannealed sample increases continuously below 400 K, the resistivity of the annealed samples has a Kondo–like minimum at about 35 K. In the temperature region between 2 K and 8 K our resistivity data for the 14 days annealed sample can be well–described by a Fermi liquid expansion $`\rho \rho _r=AT^2+BT^4`$ (see the inset of figure 3). The small upturn of the resistivity by about 0.5% below 1 K with an inkling of a maximum at 70 mK seen also in the inset of figure 3 is so far not understood. It is also observed in the unannealed sample (not shown here). Except for this small upturn below 1 K we conclude that the behavior of the resistivity of the annealed sample is clearly no longer non–Fermi liquid like with $`\rho \rho _rT`$ below 10 K as first reported by Andraka and Stewart  for the unannealed case. In fact, the observed behavior of the resistivity below 10 K is neither consistent with Fermi liquid nor with non–Fermi liquid behavior with a unique power–law dependence. Also notice that the sign of $`A`$ in our fit is unusual for f–moments forming a Kondo lattice.

Furthermore, a strong indication for a change of the ground state in this system is the change from a negative magnetoresistance at 1.8 K of the unannealed sample to a positive magnetoresistance of the annealed samples (not shown). This indicates a transition from a spin disordered compound to a ”normal” metal.

Finally, we observe that also the spin glass behavior found in unannealed UCu4Pd  is sensitive upon annealing. The spin glass temperature TSG determined from ac–susceptibility measurements shifts to lower temperature with decreasing frequency (95Hz, 995Hz) in the unannealed sample . In the annealed samples we could still observe a maximum in the ac–susceptibility at 65 mK (7 days annealed) and at 55 mK (14 days annealed), but no temperature shift with frequency (not shown). This temperature shift of the maximum might be due to a suppression of the antiferromagnetic transition of the unannealed sample as recently discussed by Körner et al. , supporting increasing order between the 7 days and 14 days annealed samples. Further work on the effects of annealing upon spin glass behavior is in progress.

## IV Conclusions

Summarizing, we have shown that in UCu4Pd intrinsic order yields a new tuning parameter for nFl behavior in the specific heat in addition to the usual parameters doping, pressure, and magnetic field. Therefore samples must also be characterized by their intrinsic order: X–ray line width and lattice parameter measurements indicate a rearrangement of the Pd atoms from the 16e sites to the 4c sites in UCu4Pd leading to increased sub-lattice order upon annealing.

Increasing order has a markedly different effect on the behavior of the specific heat and the resistivity: i) A higher degree of crystallographic order suppresses the antiferromagnetic transition temperature TN to $`<`$ 0.08 K. The best ordered sample then shows a logarithmic behavior in C/T over more than two decades in temperature from 0.08 K to 10 K.

Therefore annealing expands the range of nFl–behavior in the specific heat. ii) On the other hand, the behavior of the resistivity in the annealed sample shows a $`T^2`$–dependence with higher order corrections between 2 K and 8 K. In addition, the magnetoresistance changes sign: The non–Fermi liquid behavior ($`\rho \rho _rT`$) of the unannealed samples observed by Andraka and Stewart  in the same temperature range thus vanishes for increased crystallographic order.

Since the Kondo disorder model would predict nFl–behavior for the specific heat (logarithmic dependence of the Sommerfeld coefficient on temperature ) and for the resistivity ($`\rho \rho _rT`$ ) in about the same temperature range, this model is therefore clearly not applicable in annealed UCu4Pd. With respect to unannealed UCu4Pd, it seems furthermore plausible from figure 2a) to expect that the low temperature behavior of its Sommerfeld coefficient is due to the same physical mechanism as in annealed UCu4Pd. Thus, since disorder (at least as understood in the Kondo disorder model) is not the underlying mechanism in C/T for nFl–behavior in annealed UCu4Pd as shown above, it seems possible that it is not the origin for nFl–behavior in unannealed UCu4Pd either. With respect to the Griffiths phase model  we have found that our specific heat data is not consistent with the theoretically predicted power–law behavior $`C/TT^{1+\lambda }`$ with $`\lambda <1`$. Therefore these two theoretical models which presuppose disorder for explaining the nFl–behavior in UCu4Pd are not applicable. Our observations lead support to a quantum critical point scenario that might reconcile our specific heat, susceptibility and resistivity data. Further theoretical and experimental work, in particular $`\mu `$SR and NMR–studies on the annealed samples, is in progress in order to establish such a scenario.

Acknowledgements: We thank A. Castro Neto for valuable discussions and suggestions. S. K. has been supported by SFB 484 of the Deutsche Forschungsgemeinschaft (DFG).

## Figure Captions