format: text

text:
# On the mean field treatment of attractive interactions in nonuniform simple fluids

## 1 Introduction

In this paper we are concerned with the interplay between oscillatory excluded volume correlations and the formation of smooth liquid-vapor interfaces in simple liquids. Broadly speaking, the former arise because the harshly repulsive molecular cores cannot overlap while the latter are associated with the longer ranged and more slowly varying attractive intermolecular interactions. The theoretical challenge is to develop a general approach that can naturally explain the complete range of behavior. Early in his career Bruce Berne and coworkers  examined some of these issues in a pioneering simulation study of the effects of confining a Lennard-Jones (LJ) liquid-vapor system in a slit with a varying width. This is a long-standing and very general problem, with work extending back to the time of van der Waals, and some interesting questions remain even now.

Here we consider a simpler geometry where the different correlations manifest themselves in a particularly clear way. We study the density response and associated thermodynamic properties of a LJ liquid at a state very close to the triple point as the radius of a hard sphere solute is varied. Let $`S`$ denote the radius of the *solute cavity* within which the centers of fluid particles are excluded. As shown in recent computer simulations by Huang and Chandler  (HC), for small molecular sized $`S`$ of order unity (with the usual LJ parameter $`\sigma `$ as the unit of length) the induced density profile is highly oscillatory and resembles the radial distribution function of the uniform LJ fluid, with a large density maximum near the solute. However for much larger $`S`$ the solute approaches a hard wall, and the density response is very different: a “drying layer” of low vapor density forms near the wall, and the transition to the bulk liquid density occurs through an essentially unperturbed and smooth vapor-liquid interface. In general, as $`S`$ increases from molecular size, the contact value of the density profile decreases and mixed states with increasingly damped oscillations are seen.

We will develop a quantitative description of this process, generalizing ideas based on mean field theory that we have developed in a series of recent papers . Consider first the qualitative physics. As pointed out by Widom , in most typical configurations in the uniform LJ fluid the attractive intermolecular forces on a given particle tend to cancel in pairs between oppositely situated particles. This leaves only the excluded volume correlations induced by the harshly repulsive forces, which are well approximated by those in a uniform *hard sphere reference system*. This cancellation argument should also apply, though somewhat less accurately, to a small hard sphere solute of the same molecular size in the LJ fluid. The density response to such a fixed particle should thus closely resemble the oscillatory radial distribution in the associated hard sphere reference model, with a density *maximum* at contact.

However, since the hard core solute provides no balancing attractive forces, as its radius increases, LJ particles near the solute will increasingly experience a net “uncancelled” or “unbalanced” attractive force directed *away* from the solute arising from neighboring particles situated farther away . This unbalanced attractive interaction induces a second more slowly varying “interfacial component” in the density response that competes with the oscillatory excluded volume correlations and tends to produce a *lower* density near the solute. The final density profile results from the interplay between these components, and depends on how close the thermodynamic state of the initial uniform LJ fluid is to coexistence and on the magnitude of the unbalanced force, clearly related in this case to the size of the hard core solute.

The quantitative theory behind these ideas makes use of a fundamental approximation: the attractive interactions are described by a general *mean* or *molecular field* approach in which the locally averaged effects of the attractive intermolecular interactions are replaced by an effective single particle potential. Since the attractive interactions are relatively slowly varying, it seems plausible that such an averaged description could often provide a useful simplification. Thus, in the present case the structure of the nonuniform LJ system is approximated by that of a (hard sphere) reference fluid in an appropriately chosen *effective reference field* (ERF) comprised of two parts: the bare external field from the hard core solute and a much more slowly varying part describing the unbalanced attractive interactions.

Of course, there are inherent errors in the theory arising from the use of an ERF to describe the effects of attractive interactions. In addition, other errors can arise from the numerical determination of the properties of the reference fluid in the presence of the ERF. However, we have previously introduced a new and generally very accurate way of calculating the structure and thermodynamics of the reference fluid in the presence of a *general* external field, based on a locally optimized application of linear response theory . Thus, the latter source of error is very small (as will be verified later in this paper) and the quantitative validity of the theory for the LJ-hard core solute system depends mainly on the choice of the ERF.

In our earlier studies of nonuniform fluids in a variety of different geometries, we obtained very good qualitative and often even quantitative accuracy  using the simplest possible mean or molecular field description of the ERF, given in Eq. (3) below. However in the present case, where an essentially unperturbed vapor-liquid interface with a very large change in density can form at large $`S`$, the usual molecular field theory will yield shifted (mean field) values for the coexisting vapor and liquid densities. The main problem with the theory is not so much its description of the local density gradients, but its predictions for the thermodynamic properties of the coexisting *bulk* phases. In this paper we introduce a simple modification of the mean field expression for the ERF that incorporates information from an accurate equation of state for the *uniform* fluid. The modified expression for the ERF insures that *exact* thermodynamic results (consistent with the given uniform fluid equation of state) are found in the “hydrostatic limit” of a very slowly varying ERF, and it gives in particular accurate results for the coexisting vapor and liquid densities.

Using this modified expression for the ERF, we obtain excellent agreement with the simulations of HC for the density response to the hard sphere solute for all values of $`S`$ tested. HC have already shown that good qualitative agreement can be obtained using a simplified version of the mean field theory introduced by Lum, Chandler, and Weeks  (LCW). The LCW theory does not require detailed information about the intermolecular potential (and hence can be applied to drying transitions in more complicated liquids like water) but empirically fits a few key quantities using experimental data for properties of the liquid-vapor system. The present results suggest that the full theory can be used for quantitative calculations as well in simple liquids, requiring only the intermolecular potential and an accurate bulk equation of state along with known properties of the hard sphere reference system.

## 2 Simple molecular field approximation for the ERF

Fluid particles interact with a known external field $`\varphi (𝐫)`$ from the hard core solute. We consider a grand ensemble with fixed chemical potential $`\mu ^B`$, which determines $`\rho ^B`$, the uniform fluid density far from the solute where $`\varphi (𝐫)=0.`$ The LJ pair potential $`w(r)u_0(r)+u_1(r)`$ is separated into rapidly and slowly varying parts associated with the intermolecular *forces*  so that all the harshly repulsive forces arise from $`u_0`$ and all the attractive forces from $`u_1`$.

The structure of the nonuniform LJ system is related to that of a simpler *nonuniform* *reference fluid* , with only repulsive intermolecular pair interactions $`u_0(r)`$ (equal to the LJ repulsions) and a chemical potential $`\mu _0^B`$ corresponding to the same bulk density $`\rho ^B,`$ but in a different renormalized or *effective reference field* $`\varphi _R(𝐫)`$. (In the numerical calculations that follow, the soft-sphere reference fluid is approximated by a hard sphere fluid of appropriately chosen diameter, with the finite softness of $`u_0`$ taken into account by the usual “blip function” expansion , as described in detail in earlier work . The errors introduced by this treatment of the reference fluid are very small, and for most purposes the reference system can be thought of simply as a hard sphere system.)

What is the best choice for $`\varphi _R(𝐫)\mathrm{?}`$ Since we want the reference fluid structure to approximate that of the full fluid to the extent possible, it seems reasonable to determine $`\varphi _R(𝐫)`$ formally by the requirement that it produces a *local* (singlet) density at every point $`𝐫`$ in the reference fluid equal to that of the full nonuniform LJ fluid :

$$\rho _0(𝐫;[\varphi _R])=\rho (𝐫;[\varphi ]).$$
(1)
In practice we will make approximate choices for $`\varphi _R`$ motivated by mean field ideas. The subscript $`0`$ in Eq. (1) denotes the reference fluid and the notation $`[\varphi _R]`$ indicates that all distribution functions are functionals of the appropriate external field.

We can derive a formally exact equation for such a $`\varphi _R`$ by subtracting the first equations of the YBG hierarchy  for the full and reference systems with $`\varphi _R`$ chosen so that Eq. (1) is satisfied . The result appropriately focuses on *forces*  and can be written exactly as

$`_1[\varphi _R(𝐫_1)\varphi (𝐫_1)]`$ $`=`$ $`{\displaystyle 𝑑𝐫_2\rho _0(𝐫_2|𝐫_1;[\varphi _R])_1u_1(r_{12})}`$
$`{\displaystyle 𝑑𝐫_2\{\rho (𝐫_2|𝐫_1;[\varphi ])\rho _0(𝐫_2|𝐫_1;[\varphi _R])\}_1w(r_{12})}.`$
Here $`\rho _0(𝐫_2|𝐫_1;[\varphi _R])\rho _0^{(2)}(𝐫_1,𝐫_2;[\varphi _R])/\rho _0(𝐫_1;[\varphi _R])`$ is the *conditional* singlet density, i.e., the density at $`𝐫_2`$ given that a particle is fixed at $`𝐫_1`$.

To integrate this equation and obtain a simple expression for the effective field $`\varphi _R`$, we make some physically motivated approximations. If we assume that Eq. (1) produces similar local environments for the (identical) repulsive cores in the two fluids, which then mainly determine density correlations through excluded volume effects, then when Eq. (1) is satisfied, the conditional singlet densities in the two fluids should also be very similar. This suggests that the last term on the right side in Eq. (2) should be very small, at least at high densities where excluded volume correlations dominate. Moreover the last term clearly vanishes at very low density where $`\varphi _R=\varphi .`$

If we ignore the last term entirely, we obtain an approximate equation for $`\varphi _R`$ suggested by Weeks, Selinger and Broughton . We focus here on the even simpler *molecular field equation* that arises from further approximating the conditional singlet density $`\rho _0(𝐫_2|𝐫_1;[\varphi _R])`$ in the first term on the right side by the ordinary singlet density $`\rho _0(𝐫_2;[\varphi _R]).`$ This approximation is much better than one might at first suppose, since the main difference in these two functions occurs when $`𝐫_2`$ is close to $`𝐫_1,`$ but then for small $`r_{12}`$ the multiplicative factor $`_1u_1(r_{12})`$ (the *attractive* part of the LJ force) vanishes identically. The gradient $`_1`$ can then be taken outside the integral and the equation can be integrated. Choosing the constant of integration so that $`\varphi _R^{MF}`$ vanishes far from the solute where the density equals $`\rho ^B`$, we obtain the simple *molecular field* (MF) equation  for the ERF:

$$\varphi _R^{MF}(𝐫_1)\varphi (𝐫_1)=𝑑𝐫_2\rho _0(𝐫_2;[\varphi _R^{MF}])u_1(r_{12})+2\rho ^Ba,$$
(3)
where

$$a\frac{1}{2}𝑑𝐫_2u_1(r_{12})$$
(4)
corresponds to the attractive interaction parameter $`a`$ in the uniform fluid van der Waals equation, as discussed below.

A self-consistent solution of Eq. (3) must be found, since the ERF $`\varphi _R^{MF}`$ appears explicitly on the left side and implicitly on the right side through the dependence of the reference density $`\rho _0(r_2;[\varphi _R^{MF}])`$ on $`\varphi _R^{MF}`$. As mentioned above, we have developed new and generally very accurate ways of calculating $`\rho _0(r;[\varphi _R])`$ for a given $`\varphi _R`$. Since these methods have been described in some detail previously , we will not review them here. Using these methods, it is straightforward to solve the molecular field equation (by iteration, for example) to determine the self-consistent $`\varphi _R^{MF}`$ and the associated density $`\rho _0(r;[\varphi _R^{MF}])`$, which from Eq. (1) approximates that of the original LJ-hard core solute system.

While the simple molecular field expression (3) often gives quite satisfactory results , particularly for local density gradients, for quantitative accuracy in the present application we require a better treatment of two phase coexistence. In the limit of a *uniform* system, Eq. (3) describes all effects of attractive interactions in terms of the *constant* parameter $`a`$ as in the van der Waals equation. Indeed the theory then reduces to the generalized van der Waals theory of Longuet-Higgins and Widom , where one combines an accurate description of the uniform (hard sphere) reference system with the simple treatment of the attractive interactions in terms of the constant $`a.`$ This very simple approximation captures much essential physics and gives a good qualitative description of the uniform fluid thermodynamic properties. To achieve quantitative agreement with bulk thermodynamic properties one can replace the constant $`a`$ by a *function* $`\alpha `$ that depends (hopefully weakly, to the extent the van der Waals theory is reasonably accurate) on temperature and density . We will adopt such a strategy here, modifying Eq. (3) so that it incorporates information from a given accurate bulk equation of state and gives exact results when $`\varphi `$ and $`\varphi _R`$ are very slowly varying.

While the details differ, our method is similar in spirit to the procedure introduced in important work by van Swol and Henderson . They modified the usual mean field description of attractive interactions in density functional theory to ensure that accurate coexistence behavior is obtained in a density functional treatment of a square well fluid at a square well wall. The main idea behind both approaches is that when repulsive forces are accurately described, relatively simple modifications of the van der Waals treatment of attractive forces can produce accurate thermodynamic and structural data in the limit of very slowly varying fields or in bulk phases. The modified theory combines the simplicity of mean field theory with an accurate description of two phase coexistence.

## 3 Constant fields and chemical potential changes

Consider first the important special case where the external field $`\varphi `$ is a *constant* everywhere. Since in the grand ensemble, one encounters only the difference between the chemical potential and the external field, the effect of such a field just amounts to a change in the chemical potential of the uniform LJ fluid from $`\mu ^B`$ to $`\mu ^B\varphi ,`$ which then changes the bulk density . Equation (3) then reduces to the simple molecular field *approximation* for the analogous field $`\varphi _R`$ (or the reference chemical potential change) that from Eq. (1) is supposed to produce the same density change in the uniform reference fluid. In general there will be errors arising from the molecular field approximation even in this simple limit. However, if an accurate equation of state for the uniform fluid is known, we can calculate these density changes exactly and this will suggest how to modify Eq. (3) so that it gives exact results for very slowly varying fields.

To that end, let $`\mu (\rho )`$ denote the chemical potential as a function of density $`\rho `$ for the uniform LJ fluid. (This also depends on the temperature, but we are interested in density variations along particular isotherms, so we will not indicate the temperature dependence explicitly.) We determine $`\mu (\rho )`$ from the accurate 33-parameter equation of state  given by Johnson, et al. This provides a very good global description of the stable liquid and vapor phases in the (cutoff and shifted) LJ fluid used in the simulations of HC and provides a smooth interpolation in between by using analytic fitting functions. Thus it naturally produces a modified “van der Waals loop” in the two phase region and seems quite appropriate for our use here in improving the simplest mean field or van der Waals description of the uniform fluid. One could alternately use the method of van Swol and Henderson  and achieve essentially identical results. If $`\varphi (\rho )`$ denotes the magnitude of the uniform external field that yields a particular density $`\rho `$ when applied to the LJ system with chemical potential $`\mu ^B`$, then its value is exactly given by

$$\varphi (\rho )=\mu ^B\mu (\rho ).$$
(5)
Using known properties of the hard sphere fluid, we also have essentially exact expressions for $`\mu _0(\rho ),`$ the chemical potential as a function of density in the uniform reference fluid, and for the inverse function $`\rho _0(\mu ),`$ giving the uniform reference density as a function of chemical potential. The analogous uniform reference field yielding the same density $`\rho `$ as in Eq. (5) is then exactly given by

$$\varphi _R(\rho )=\mu _0^B\mu _0(\rho ).$$
(6)
Finally the exact $`\mu (\rho )`$ and $`\mu _0(\rho )`$ can be related through a function $`\alpha (\rho )`$ defined so that

$$\mu (\rho )=\mu _0(\rho )2\rho \alpha (\rho ).$$
(7)
Thus the exact chemical potentials in the uniform LJ and reference systems are related in the same way as predicted by the simple mean field approximation of Eq. (3) except that the constant $`a`$ in Eq. (4) is replaced by a (temperature and density dependent) function $`\alpha (\rho )`$ chosen so that Eq. (7) holds. Since even the simplest mean field theory is qualitatively accurate we expect that the ratio $`\alpha (\rho )/a`$ will be of order unity and rather weakly dependent on density and temperature.

In Fig. 1 we give the calculated effective $`\alpha (\rho )`$ from the equation of state of Johnson et al. , together with the value $`a`$ of the van der Waals constant in Eq. (4) for a number of different isotherms. The value of the constant $`a=2\pi \underset{0}{\overset{\mathrm{}}{}}𝑑rr^2u_1(r)`$ seems to be a good overall compromise and may suffice for many qualitative purposes. Our derivation resulting from approximately integrating the force balance in Eq. (2) indicates why this expression with $`u_1`$ giving all the attractive *forces* and integrated from $`0`$ to infinity is appropriate. This gives much better results than the earlier intuitive suggestion $`a^{}=2\pi \underset{\sigma }{\overset{\mathrm{}}{}}𝑑rr^2w(r),`$ which integrates only the *negative part* of the LJ potential from the point $`\sigma `$ where the potential changes sign . However the true $`\alpha (\rho )`$ exhibits some noticeable variations in density and temperature, illustrating the need for an accurate equation of state for quantitative accuracy. We note that some of the thermodynamic states plotted in Fig. 1 are in the two phase region of the bulk LJ fluid and thus are not thermodynamically or sometimes even mechanically stable. The analytic equation of state provides a smooth interpolation through these states and produces behavior analogous to the van der Waals loop in the simplest mean field description. By introducing a function $`\alpha (\rho )`$ that reduces to the *constant* $`a`$ in the van der Waals theory we hopefully have removed much of the (temperature and) density dependence otherwise arising from the loops, leaving the relatively smooth $`\alpha (\rho )`$ suitable for interpolation .

## 4 Slowly varying fields and the hydrostatic density

These results for a strictly constant field can be used to determine simple and accurate approximations for very slowly varying fields. In particular when the field $`\varphi _R(𝐫)`$ varies so slowly that it is essentially *constant* over the range of a correlation length in the reference fluid, then the density $`\rho _0(𝐫_1;[\varphi _R])`$ at a given point $`𝐫_1`$ is very accurately approximated  by the *local hydrostatic density* $`\rho ^{𝐫_1}\rho _0(\mu _0^{𝐫_1})`$, the density of the *uniform* reference fluid at the shifted chemical potential $`\mu _0^{𝐫_1}\mu _0^B\varphi _R(𝐫_1).`$ Since in the hydrostatic approximation the density at any given point responds only to the *local* value of the field, all the equations defined above for constant fields and uniform densities also hold for very slowly varying fields and the corresponding hydrostatic densities.

By subtracting Eq. (5) from Eq. (6) and using Eq. (7) we then arrive at an equation analogous to the simple mean field equation (3) that gives essentially exact results in the hydrostatic limit of very slowly varying fields:

$$\varphi _R(𝐫_1)=\varphi (𝐫_1)2\rho ^{𝐫_1}\alpha (\rho ^{𝐫_1})+2\rho ^B\alpha (\rho ^B).$$
(8)
## 5 Modified molecular field equation

We now compare this exact result to the hydrostatic limit of the simple mean field equation (3). It is easy to see that this takes the same form as Eq. (8) except that the function $`\alpha (\rho )`$ is replaced by the constant $`a.`$ Our goal is to modify Eq. (3) so that in the hydrostatic limit it reduces exactly to Eq. (8), while still giving reasonable results for more rapidly varying fields.

There is no unique way to do this, but the following simple prescription seems most natural, and gives our final result, which we will call the  *modified molecular field* (MMF) approximation for the ERF:

$$\varphi _R(𝐫_1)\varphi (𝐫_1)=\frac{\alpha (\rho ^{𝐫_1})}{a}𝑑𝐫_2\rho _0(𝐫_2;[\varphi _R])u_1(r_{12})+2\rho ^B\alpha (\rho ^B).$$
(9)
Thus the molecular field integral in Eq. (3) is multiplied by a factor $`\alpha (\rho ^{𝐫_1})/a`$ of order unity that depends on $`𝐫_1`$ through the dependence of the hydrostatic density $`\rho ^{𝐫_1}`$ on the local value of the field $`\varphi _R(𝐫_1),`$ and the constant of integration $`2\rho ^Ba`$ is replaced by the appropriate limiting value of the modified integral. Note that the hydrostatic density $`\rho ^{𝐫_1}`$ remains smooth and relatively slowly varying even when $`\varphi _R(𝐫_1)`$ contains a hard core. The nonlocal oscillatory excluded volume correlations that can exist in the full density $`\rho _0(𝐫_1;[\varphi _R])`$do not appear in $`\rho ^{𝐫_1}`$ because of the strictly local response to $`\varphi _R`$.

The next sections report results for the structure and thermodynamics of the the LJ-hard core solute system studied by HC with the ERF determined from Eq. (9). We use methods based on generalized linear response theory  to determine the structure of the nonuniform reference fluid.

## 6 Results

### 6.1 Structure

To make contact with the work of HC, we study the LJ liquid at a state near the triple point with $`\rho ^B=0.70`$ and $`T=0.85`$ in the presence of a hard sphere solute. (We use the standard LJ reduced units.) By definition the solute centered at the origin interacts with the LJ particles through the hard core potential:

$$\varphi (r;S)=\{\begin{array}{cc}\mathrm{},\hfill & rS,\hfill \\ 0,\hfill & r>S.\hfill \end{array}$$
(10)
The MMF theory discussed above allows us to reduce this problem to that of the *reference* fluid in the presence of the effective field $`\varphi _R(r;S)`$ satisfying Eq. (9). By construction from Eq. (1), the density profile of the full LJ fluid subjected to the “bare” external field $`\varphi (r;S)`$ and the profile of the reference LJ fluid in the presence of the external field $`\varphi _R(r;S)`$ “dressed” by the attractive interactions are supposed to be identical to each other.

We have calculated *self-consistently* the ERF $`\varphi _R(r;S)`$ and the associated density response $`\rho _0(r;S)`$ of the reference fluid, solving Eq. (9) by iteration. In Fig. 2 we compare these results for the density profiles in the presence of the hard sphere solutes with $`S`$ equal to $`1.0`$, $`2.0`$, $`3.0`$ and $`4.0`$ in reduced units with the simulation results  of the same LJ system by HC. There is very good agreement between theory and simulation.

Figure 3 shows the corresponding ERF’s obtained in these calculations. For small solutes with $`S`$ less than about $`0.7,`$ attractive interactions do not give rise to any substantial modification of the bare external field, as can be seen from the plot of $`\varphi _R(r;S)`$ for $`S=0.5`$. (Clearly, for $`S=0`$ there are no solute induced iteractions of any kind and the profile reduces to the constant $`\rho ^B`$). However, the effects due to unbalanced attractions begin to become important even for $`S=1.0`$, which is about the same size as the LJ core, and all larger sizes give rise to a very strong and relatively soft repulsion in $`\varphi _R(r;S)`$. The corresponding density profiles show pronounced depletion near the surface of the solute, characteristic of surface induced drying.

To demonstrate the importance of the renormalization of the effective field we show in Fig. 4 the density profiles of the reference fluid in the presence of the *bare* external field $`\varphi (r;S)`$, thus neglecting all effects due to the attractions. These density profiles reflect the characteristic packing effects associated with excluded volume correlations, which result in a local ordering (or layering) of fluid particles next to the solute. We would also expect such behavior from the LJ system at very high temperatures, but for the state near the triple point considered here the strong ordering is absent even though the density is high. The unbalanced attractive interactions generate an additional strong repulsive component in the ERF as shown in Fig. 3. This is capable of destroying the ordering and effectively places a layer of the metastable low density vapor next to the solute’s surface for states near coexistence.

In the self-consistent theory discussed above we make approximations both in determining the ERF and in calculating the structure of the reference fluid the presence of a given ERF. To make sure that in the present case the agreement with the simulation results for the LJ fluid given in Fig. 2 is not the result of some fortuitous cancellation of errors associated with our treatment of the reference fluid, we have carried out computer simulations  of the *reference* fluid structure in the presence of the self-consistently determined ERF’s. This directly tests our theory  for reference system correlations in the presence of a given external field. Results of these calculations are shown in Fig. 5 and demonstrate a very good agreement between our theory for the reference fluid structure  and the simulations.

### 6.2 Solvation free energy

Another quantity of great interest is the free energy of the nonuniform system. This is the main focus of attention in density functional theory , where one tries to approximate the “intrinsic” free energy as a functional of the nonuniform singlet density. The equilibrium density and value of the free energy for a given external field and set of intermolecular interaction potentials is determined by minimizing an appropriate free energy functional of the density and those potentials . By starting from the free energy, certain exact *sum rules* relating integrals of those correlation functions to the thermodynamic properties are automatically and consistently satisfied . Of course consistency does not necessarily imply correctness — a poor functional will give consistently poor values for the thermodynamic properties — and the physical implications of particular approximations made in DFT can sometimes be hard to understand.

In contrast, our approach focuses first on the liquid *structure*, as do most integral equation methods. We believe this permits physical insight to play a more direct role. However, since we can determine the density response to an arbitrary external field, the free energy can be easily calculated from a coupling parameter type integration that connects some initial state (e.g., the bulk fluid) whose free energy is known to the final state as the field is varied. Of course, our structural predictions, though generally accurate, are not exact and different routes to the same thermodynamic property can give different answers. Fortunately, in the present case there is a very simple route to the free energy of the nonuniform LJ system that uses structural features that we know from simulations are accurately determined.

### 6.3 Virial route

The basis for our calculation lies in the definition of the grand canonical free energy (thermodynamic potential) in the LJ fluid, which is related to the partition function $`\mathrm{\Xi }`$ through the following equation:

$$\mathrm{\Omega }=k_BT\mathrm{ln}\mathrm{\Xi }.$$
(11)
We also use the property that the partition function is the generating function for the density correlation functions; in particular, the density profile in the presence of a general field $`\varphi `$ is given by

$$\rho (𝐫;\mu ^B,[\varphi ])=k_BT\frac{\delta \mathrm{ln}\mathrm{\Xi }(\mu ^B,[\varphi ])}{\delta [\mu ^B\varphi (𝐫)]}.$$
(12)
We now specialize to the hard core solute $`\varphi (r;S)`$ and introduce a coupling-parameter dependent external field $`\varphi _\lambda (𝐫)\varphi (𝐫/\lambda )`$ with corresponding density $`\rho _\lambda (𝐫)\rho (𝐫;\mu ^B,[\varphi _\lambda ])`$. This $`\lambda `$-dependence corresponds to growing the range of the external field from zero at $`\lambda =0`$ to the full external field at $`\lambda =1`$. This construction is the basis of scaled particle theory . As is well-known, the free energy change for the hard core field then takes the particularly simple form:

$$\beta \mathrm{\Delta }\mathrm{\Omega }_S=4\pi S^3\underset{0}{\overset{1}{}}𝑑\lambda \lambda ^2\rho _\lambda (\lambda S^+),$$
(13)
which requires only the *contact value* $`\rho _\lambda (\lambda S^+)`$ of the density profile. This is very accurately given by the theory described above. To use this “virial route”, we can discretize the $`\lambda `$-integration and calculate the density profile for several values of $`\lambda `$ at the fixed bulk chemical potential $`\mu ^B`$.

In Fig. 6 we plot the contact values as a function of the size $`S`$. For solvent fluids with purely repulsive interactions, such as the hard sphere fluid, the contact value monotonically increases because the fluid particles increasingly tend to order next to the surface of the solute and the density profile always has a sharp peak at contact. A nonmonotonic dependence of the contact values is a very characteristic feature of surface induced drying, and has been predicted for the case of water in contact with hydrophobic objects by Stillinger  and in the present case by the LCW theory.

Using Eq. (13) we obtain the dependence of solvation free energy on the size of the hard sphere solute. The free energy per unit surface area of the solute $`\mathrm{\Delta }\mathrm{\Omega }_S/4\pi S^2`$ we obtain is shown in Fig. 7. For small solutes attractive interactions do not play an important role and the solvation free energy agrees well with a pure hard sphere model, which completely neglects attractions by using the bare solute potential, as shown by the dotted line. At the solute size of about $`0.7`$ the behavior changes drastically and the reduced free energy rapidly crosses over to the practically constant plateau in agreement with the simulation results. The small slope of the curves in Fig. 7 for large $`S`$ can be understood by separating the free energy into volume ($`V_S=4\pi S^3/3`$) and surface ($`A_S=4\pi S^2`$) contributions  as discussed by HC:

$$\mathrm{\Delta }\mathrm{\Omega }_SV_Sp^B+A_S\gamma _S.$$
(14)
The first term in this expression corresponds to the work required to remove liquid particles from the volume occupied the solute, where $`p^B`$ is the bulk liquid pressure, and is very small for the values of $`S`$ considered here. The second term determines the cost of forming the liquid-solute interface and is proportional to the interface tension $`\gamma _S`$, which is essentially independent of the solute size for large solutes.

## 7 Conclusions

The theory described here combines an accurate treatment of correlations induced in the nonuniform hard sphere reference system by a general external field with a mean field description of the effects of attractive interactions in terms of an appropriately chosen ERF. We have shown here how relatively simple modification of the simple molecular field expression for the ERF can ensure that accurate results are found in the hydrostatic limit of very slowly varying fields and for properties of coexisting bulk phases. Of course there are important limitations inherent in any mean field treatment of long wavelength fluctuations, such as those seen near the critical point or the capillary wave fluctuations of the liquid-vapor interface . However in many other applications where such fluctuations do not play an important role, mean field theory provides a simple and often very accurate starting point. In future work we will discuss properties of the liquid-vapor interface and of fluids confined in slits and cylinders from this perspective.

It is a pleasure to dedicate this paper to Bruce Berne on the happy occasion of his 60th birthday. This work was supported in part by NSF Grant No. CHE9528915.


text:
# X-ray emission from middle-aged pulsars

## 1 Introduction

X-ray emission from cooling neutron stars (NSs) had been predicted and studied in detail already in the early 1960’s (Chiu & Salpeter 1964 ; Tsuruta 1964 ). However, it was not until the late 70’s that the Einstein Observatory first detected radiation from isolated NSs, and until the ROSAT era (the 1990’s) that more detailed information on the spectral characteristics and pulsed fraction ($`Pf`$) could be obtained.

Thermal radiation due to cooling can be best observed in a small subset of pulsars which are middle-aged and have characteristic lifetimes $`\tau 10^5`$ yr. In younger pulsars, such as the Crab, X-rays of magnetospheric origin often dominate over thermal emission, while older ones are too cold to have any observable radiation due to cooling. So far, the three best candidates for which thermal emission has likely been detected are PSR 1055-52, PSR 0656+14 and PSR 0630+18.

Apart from a non-thermal, power-law component dominating the high-energy tail (this is particularly evident in PSR 0656+14; Greiveldinger et al. 1996 ) and which is thought to be due to magnetospheric processes, the bulk of the emission between $`0.11`$ keV is well fitted by a double blackbody at two different temperatures. The ratio of the area of the hotter component to the colder one is typically very small ($``$ a few $`\times 10^5`$ to a few $`\times 10^3`$). The colder component has been interpretated as being due to thermal cooling, while the hotter one most likely arises from heated polar caps, with the extra heating being due to bombardment by high-energy particles (Greiveldinger et al. 1996 ).

An important property of the radiation from all three pulsars (for a review see Becker & Trumper 1997) is a phase shift of $`80^{}100^{}`$ at a transition point which is around 0.5-0.6 keV. The pulsed fraction also shows transitions in its behavior in this region for both PSR 1055-52 and PSR 0656+14. Below $`0.5`$ keV, it is roughly constant at the 10% level, while it rapidly increases after that point to about 20-30% for PSR 0656+14 and to about 80-90% for PSR 1055-52. On the other hand, the pulsed fraction from Geminga shows an intriguing feature: in the PSPC channels 8-28 (i.e. roughly at energies below 0.3 keV) the $`Pf`$ is much larger than in channels 28-53 (roughly corresponding to the energy range 0.3-0.5 keV): 33% versus 20%. However, as discussed by Page 1995  and Page & Sarmiento 1996 , no matter what the surface temperature distribution is, as long as it is not uniform, blackbody emission always gives an increase of the $`Pf`$ with energy. Detailed modeling by Shibanov et al. 1995  shows that realistic atmospheres are able to produce a slight decrease of the $`Pf`$ with increasing energy, but it is still much smaller than what is observed. A model which could account for the large observed decrease has been presented by Page, Shibanov & Zavlin 1995 . They assume the presence of warm, magnetized plates on the surface of Geminga, surrounded by cold, unmagnetized regions, with the warmer plates emitting a softer spectrum than the surrounding areas. In this model, Geminga would therefore have different characteristics than the two other objects in the same class.

In this Letter, we present a simple model for the emission in the $`0.11`$ keV range that is able to account for the properties of all three pulsars within a unified framework. We model the cooling component as a blackbody modified by a light-element atmosphere; as discussed in the literature (Heyl & Hernquist 1998b ; Zavlin, Pavlov & Shibanov 1996 ; Pavlov et al. 1994 ), this radiation, which emerges through an atmosphere heated from below, is most likely emitted in a “pencil” beam (Pavlov et al. 1994 ; Zavlin, Pavlov & Shibanov 1996 ). The hotter component, due to particle bombardment of the polar caps, is produced in an atmosphere that is hotter on top, i.e. at lower densities. This is due to the fact that the tenuous upper reaches of the neutron-star atmospheres are extremely inefficient emitters; therefore, as the bombarding particles deposit energy in this region, it heats up dramatically to maintain radiative balance. Meanwhile, the denser layers of the atmosphere which have a larger opacity heat up relatively little. The problem of atmospheres heated from above has been discussed in the literature in the context of slow accretion onto neutron stars (Zane, Turolla & Treves 2000 ) and the X-ray illumination of normal stars in X-ray binaries (e.g.  Milgrom & Salpeter 1975 ). It has been shown that the radiation from an atmosphere illuminated from above may be emitted in a “fan” pattern (e.g.  Milgrom & Salpeter 1975 ). In fact, as one moves away from the normal, the effective photosphere of the atmosphere moves to lower densities due to geometrical effects. In a bombarded atmosphere, these lower densities are hotter, so one sees a hotter and brighter spectrum as one approaches grazing incidence.

We show that a combination of pencil beaming from the cooling component and fan beaming from the hotter polar caps is able to account not only for the magnitude of the pulsations in PSR 1055-52 and PSR 0656+14, but also for the phase shift between the soft and the hard components. Moreover, this same model is able to account also for the decrease of the pulsed fraction observed in the soft X-ray emission of Geminga, without requiring any special model for the composition of the surface of this object.

The paper is organized as follows: in §2, we describe in detail the X-ray spectral model; specific applications to PSR 1055-52, PSR 0656+14 and Geminga are discussed in §3; our results are summarized and discussed in §4.

## 2 X-ray spectrum

We model the local emission from the surface of the star, $`n(E,T)`$, as a blackbody spectrum modified by the presence of an atmosphere, for which we adopt the semianalytical model developed by Heyl & Hernquist 1998c  with a geometric generalization to study the spectral intensity away from the normal. The dependence of our results on the atmospheric composition will be discussed in §3. We assume that the neutron star has a dipolar magnetic field. With this field geometry and a sufficiently strong intensity, the flux transmitted through the envelope can be approximated as $`F\mathrm{cos}^2\psi `$, where $`\mathrm{cos}^2\psi =4\mathrm{cos}^2\theta _p/(3\mathrm{cos}^2\theta _p+1)`$ is the angle between the radial direction and the magnetic field (Greenstein & Hartke 1983 ; Heyl & Hernquist 1998a , 2000); here

$$\theta _p=\mathrm{arccos}(\mathrm{cos}\theta \mathrm{cos}\alpha +\mathrm{sin}\theta \mathrm{sin}\alpha \mathrm{cos}\varphi )$$
(1)
is the angle that the magnetic pole makes with the normal to the star at position $`(\theta ,\varphi )`$, while $`\alpha `$ is the angle that it makes with the line of sight. If $`\xi `$ is the angle between the magnetic dipole and the rotation axis, and $`\chi `$ the angle between the observer’s direction and the rotation axis, then the angle $`\alpha `$ is given by

$$\alpha =\mathrm{arccos}(\mathrm{cos}\chi \mathrm{cos}\xi +\mathrm{sin}\chi \mathrm{sin}\xi \mathrm{cos}\gamma ),$$
(2)
with $`\gamma `$ being the phase angle.

Outside of the heated polar caps, the local temperature on the star due to thermal cooling is given by

$$T_{\mathrm{th}}(\theta ,\varphi )=T_p\left[\frac{4\mathrm{cos}^2\theta _p}{3\mathrm{cos}^2\theta _p+1}(0.75\mathrm{cos}^2\theta _p+0.25)^{0.2}\right],$$
(3)
where, following Heyl & Hernquist 1998a , we have assumed a further dependence of the flux on magnetic field strength scaling as $`B^{0.4}`$. In Equation (3), $`T_p`$ is the temperature that the pole would have if there were no reheating.

Now, let $`\beta `$ be the angular size of the radius of the polar caps, which are centered around the poles. The hot spot region (in the upper hemisphere) is defined by

$$\theta \beta ,\text{if}\alpha =0$$
(4)
and

$$\{\begin{array}{cc}\alpha \beta \theta \alpha +\beta \hfill & \\ \varphi \varphi _p\text{or}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}2}\pi \varphi \varphi _p,\text{if}\alpha 0\text{and}\beta \alpha \hfill & \end{array}$$
(5)
where

$$\varphi _p=\mathrm{arccos}\left[\frac{\mathrm{cos}\beta \mathrm{cos}\alpha \mathrm{cos}\theta }{\mathrm{sin}\alpha \mathrm{sin}\theta }\right].$$
(6)
Here, we have chosen a system of coordinates such that the prime meridian passes through the North pole of the star, at the center of the upper polar cap. The hot spot in the lower hemisphere is defined by the condition $`\theta \pi \beta `$, if $`\alpha =0`$, and by Equation (5) and Equation (6) with the substitutions $`\alpha \pi \alpha `$ and $`\varphi \pi \varphi `$, if $`\alpha 0`$.

Inside the hot spot, we assume the temperature to be constant, $`T_{\text{hs}}(\theta ,\varphi )=\text{c}onstT_{\mathrm{hs}}`$. As already discussed, we adopt pencil beaming for the thermal component and fan beaming for the radiation produced in the polar caps, and, for simplicity, we parameterize them with the functions

$$\{\begin{array}{cc}f_{\text{th}}(\delta )\mathrm{cos}^{n_1}(\delta )\hfill & \\ f_{\text{hs}}(\delta )\mathrm{sin}^{n_2}(\delta )\hfill & \end{array},$$
(7)
where $`\delta `$ is the angle that a photon emitted at a colatitude $`\theta `$ on the star makes with the normal to the surface at the moment of emission. The relation between $`\theta `$ and $`\delta `$ (which is a consequence of the general relativistic effects of light deflection) is given by the ray-tracing function (Page 1995)

$$\theta (\delta )=_0^{R_s/2R}x𝑑u/\sqrt{\left(1\frac{R_s}{R}\right)\left(\frac{R_s}{2R}\right)^2(12u)u^2x^2},$$
(8)
having defined $`x\mathrm{sin}\delta `$. The value of $`n_1`$ in the above equation is found to be roughly unity for an unmagnetized light-element atmosphere (Zavlin, Pavlov & Shibanov 1996 ) and slightly larger for a magnetized light-element atmosphere (Pavlov et al. 1994 ), while it is about 0.5 for an atmosphere made of heavy elements (Rajagopal, Romani & Miller 1997 ). In our calculations we restrict ourselves to this range of values. The exponent $`n_2`$, on the other hand, is treated as a fit parameter and is allowed to vary from source to source, as one would expect the extent of fan beaming in the bombarded atmosphere to depend on the spatial and energy distribution of the bombarding particles which are unknown (and likely to vary among different objects).

The calculation of the emitted spectrum fully accounts for the consequences of gravitational bending of light and gravitational redshift. Let us define $`e^{\mathrm{\Lambda }_s}=\sqrt{1R/R_s}`$, where $`R`$ is the radius of the neutron star, and $`R_s=2GM/c^2`$ is its Schwarzschild radius. A range of radii compatible with the currently available models for the NS equation of state requires $`2(R/R_s)4`$. Here we take $`R`$ in this range, with $`M=1.4M_{}`$. Let $`D`$ be the distance from the star to the observer, and $`N_\text{H}`$ the intervening column density. The flux measured by an observer at infinity is then given by the sum of the contributions from the thermal component and from the polar caps

$`F(E;\gamma )`$ $`=`$ $`{\displaystyle \frac{\pi R^2e^{\mathrm{\Lambda }_s}}{4\pi D^2}}e^{\sigma (E)N_\text{H}}{\displaystyle _0^1}2x𝑑x{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\varphi }{2\pi }}`$ (9)
$`\times `$ $`\left\{{\displaystyle \frac{1}{kT_p}}\sigma T_{\text{th}}^4(\theta ,\varphi )n[Ee^{\mathrm{\Lambda }_s};T_{\text{th}}(\theta ,\varphi )]+{\displaystyle \frac{1}{kT_{\text{hs}}}}\sigma T_{\text{hs}}^4(\theta ,\varphi )n[Ee^{\mathrm{\Lambda }_s};T_{\text{hs}}(\theta ,\varphi )]\right\},`$
in units of phot cm-2 s-1 keV-1. Note that the dependence on $`\gamma `$ comes in through Eqs. (1) – (3).

Finally, the energy dependent pulsed fraction is defined by

$$Pf(E)=\frac{F^{\text{max}}(E)F^{\text{min}}(E)}{F^{\text{max}}(E)+F^{\text{min}}(E)},$$
(10)
where $`F^{\text{max}}(E)`$ and $`F^{\text{min}}(E)`$ are, respectively, the maximum and minimum flux over a rotation period of the star. The phase of the pulsation, $`Ph(E)`$, is defined as the phase angle $`\gamma `$ at which the flux is maximum at energy $`E`$.

## 3 Application to PSR 1055-52, PSR 0656+14 and PSR 0630+18

For each pulsar, we considered three types of models: a two-component blackbody (of the type expected from iron atmospheres), or a two-component light-element atmosphere, and for the latter we considered two different types of opacities. The emergent spectrum from iron atmospheres is similar to a blackbody in its gross properties (Rajagopal, Romani & Miller 1997 ) and we assume that it is moderately beamed as discussed in §2. The two-component light-element atmosphere is substantially bluer than a blackbody and also exhibits a hard tail. We have used two different semianalytic atmospheres of Heyl & Hernquist 1998c  with opacities decreasing as $`\nu ^3`$ and $`\nu ^1`$. The former model is much bluer than a blackbody and exhibits a strong hard tail. The latter is most appropriate to model magnetized light-element atmospheres (e.g.  Pavlov et al. 1994 ); it is intermediate between the $`\nu ^3`$-model and the blackbody.

The effective temperatures (parameterized through $`T_p`$ and $`T_{\mathrm{hs}}`$) and the size of the hot region (called $`\beta `$ in our paper) correponding to the same effective blackbody temperatures are different in the three models, and, to precisely determine them in each case, we generated a simulated spectrum (using the software XSPEC, Arnaud 1996 , and the $`ROSAT`$ PSPC detector) of each object, using the temperature values of the double-blackbody fits and the ratio between the areas of the two components given in the literature. This spectrum was then fitted with the phase averaged spectrum of our model111This is given by $`\frac{1}{2\pi }_0^{2\pi }𝑑\gamma F(E;\gamma )`$. to determine $`T_p`$, $`T_{\mathrm{hs}}`$ and $`\beta `$. These values were then used to compute pulsed fractions and phases with some choices of $`n_1`$ and $`n_2`$.

We found that an iron-type atmosphere, which has a rather mild degree of beaming, was unable to reproduce the degree of modulation observed in the softest part of the spectrum (dominated by the thermal component222Note, however, that this result is dependent on the assumption of a dipolar field. The value of the pulsed fractions can be increased if a quadrupole component is added (Page & Sarmiento 1996 ).). A light element atmosphere with an opacity $`\nu ^3`$ was able to reproduce reasonably well the pulsed fractions in the soft component of the spectrum. However, the tail of the thermal component is generally so strong that a second hotter component was not statistically required by the data in the case of PSR 0656 and PSR 1055, and it was much smaller than the one required by the double blackbody fit in the case of Geminga333 Note that in the case of Geminga, similar results had been found by Meyer et al. 1990 .. Such a model was unable to produce the sudden transition in the behaviour of the pulsed fractions accompanied by a phase shift around 0.5 keV which is observed in PSR 0656 and PSR 1055. Moreover, such models implied surface areas much larger than what would be consistent with any reasonable equation of state for the neutron star. A two-component atmosphere model which is less blue ($`\nu ^1`$) than the former required a second hotter component in the spectral fitting, and was able to reproduce the behaviour of the observed pulsed fractions and phase shifts. More details on each of the object that we considered are given in the following.

### 3.1 PSR 1055-52

This pulsar has a characteristic age $`\tau =P/2\dot{P}5\times 10^5`$ yr, and a rotational energy loss rate of $`\dot{E}3\times 10^{34}`$ erg/sec. ROSAT observations (Ögelman & Finley 1993 ) clearly show that the pulses have an energy dependent phase and pulsed fraction, with a transition at $`0.5`$ keV. The pulsed fraction below the transition is roughly constant at a value on the order of 7-8%, while after the transition it rapidly increases to a maximum of 0.85 around 1 keV. The hard photons lead the soft ones by an angle $`120^{}`$. A blackbody fit to the spectrum requires two components, with temperatures $`T_{\text{soft}}8\times 10^5`$ K and $`T_{\text{hard}}3.7\times 10^5`$ K, and a ratio between the areas of $`A_{\text{soft}}/A_{\text{hard}}3\times 10^5`$ (Greiveldinger et al. 1996 ). The angles $`\xi `$ and $`\chi `$ have been inferred by Malov 1990  to be both about $`30^{}`$. With these values, and $`T_\text{p}`$, $`T_{\text{hs}}`$ and $`\beta `$ calibrated on the phase averaged spectrum, we show the derived pulsed fractions and phase for this pulsar in Figure 1, having assumed $`n_1=1`$ and $`n_2=2`$.

### 3.2 PSR 0656+14

PSR 0656+14 was discovered in an Einstein satellite survey of ultrasoft sources (Cordova et al. 1989 ). It has a relatively young spin down age of $`1.1\times 10^5`$ yr and a rotational energy loss rate $`\dot{E}3.8\times 10^{34}`$ erg/sec. The combined ASCA and ROSAT spectrum reveals two blackbody components with $`T_{\text{soft}}8\times 10^5`$ K, and $`T_{\text{hard}}1.5\times 10^6`$ K, and shows evidence that a power-law component is needed to account for higher energy photons. The ratio of the hot polar cap area to the neutron star surface area is $`7\times 10^3`$ (Greiveldinger et al. 1996 ).

As for PSR 1055-52, both the phase and the pulsed fraction are energy dependent. Below $`0.5`$ keV, the pulsed fraction stays a little under $`10\%`$, and it increases modestly in the range $`0.51`$ keV to around 20-30 % at 1 keV. The soft and the hard components are shifted in phase by about $`85^{}`$ with respect to each other. The angles $`\xi `$ and $`\chi `$ have been estimated by Malov 1990  to be both about $`35^{}`$. With these values, and $`T_\text{p}`$, $`T_{\text{hs}}`$ and $`\beta `$ derived from spectral fitting, we show the pulsed fraction and phase that our model predicts with the choice $`n_1=0.8`$ and $`n_2=0.2`$ in Figure 2.

### 3.3 PSR 0630+18 (Geminga)

Geminga was first observed as a high-energy $`\gamma `$-ray source with the SAS-2 satellite in the 100 MeV band (Fichtel et al. 1975 ). It was only in 1992 that it was observed in the X-ray band with ROSAT (Halpern & Holt 1992 ). Its period and spin-down rate (Bertsch et al. 1992 ; Hermsen et al. 1992 ) yield a dynamic age of $`3.4\times 10^5`$ yr, and a rotational energy loss rate $`\dot{E}3.5\times 10^{35}`$ erg/sec. A spectral analysis of the ROSAT data made by Halpern & Ruderman 1993  shows that the X-ray spectrum consists of two blackbody components with $`T_{\text{soft}}5\times 10^5`$ K, and $`T_{\text{hard}}3\times 10^6`$ K. Both components are modulated at the pulsar rotation period, but the harder X-ray pulse leads the soft pulse by about $`105^{}`$ in phase. Geminga is believed to be an orthogonal rotator, with $`\xi =\chi 90^{}`$ (Malov 1990 ).

Unlike PSR 1055-52 and PSR 0656+14, the pulsed fraction observed in Geminga shows a decrease with energy at low energies: the amplitude of the pulsations in the PSPC channels 8-28 (i.e. roughly at energies below 0.3 keV) is much larger than in channels 28-53 (roughly corresponding to the energy range 0.3-0.5 keV): 33% versus 20%. Figure 3 shows that a significant decrease of the pulsed fractions in the 0.3-0.5 keV energy range can be reproduced with the model described in §2: here we have taken $`n_1=1.5`$ and $`n_2=3`$, while all other model parameters have again been calibrated from a fit to the phase-averaged spectrum444Note that a similar behaviour for the $`Pfs`$ in the softest part of the spectrum could be obtained with a more modest beaming (i.e. smaller $`n_1`$), but a quadrupole component of the $`B`$ field added to the dipole (Page96).. The decrease in the pulsed fraction occurs in the region where the harder component starts to overtake the softer one. The hard component brings more photons at phase angles in which the soft component has less photons, therefore leading to a decrease of the intensity fluctuations produced by the thermal component alone. This effect happens in the same way also in the other two examples that we showed, but due to the less intense beaming of the harder component assumed there, the decrease is not so pronounced. Interestingly, a slight decrease in the pulsed fraction before its rapid increase can be observed also in the data for these other two pulsars (e.g. Ögelman 1995 ).

Our model for the phase shifts and pulsed fraction in Geminga is similar to that proposed by Halpern & Ruderman 1993 . Page, Shibanov & Zavlin 1995  later argued that this picture could not work because the flux in the hard component is more than an order of magnitude less than that in the soft component in the $`0.30.5`$ keV energy range, where the decrease in the $`Pf`$ is observed. However, the pencil-beaming reduces the flux of the soft component at the same phase of the rotation where the fan-beaming of the hard component increases its flux; therefore, although the phase averaged fluxes in the two components are not comparable, the instantaneous fluxes at $`\gamma 90^{}`$ can be, and the hard component can begin to influence the soft one well before it becomes dominant in the phase-averaged spectrum.

## 4 Summary and Discussion

We have proposed a simple model that is able to account for properties of the X-ray emission observed in the three middle aged pulsars PSR 1055-52, PSR 0656+14 and Geminga within a unified framework. The emission in the $`0.11`$ keV energy range is believed to be due to the combination of thermal radiation from the entire star and emission from heated polar caps. We have argued that, while pencil beaming is expected for the thermal, soft component, the hard component from the polar caps is expected to be beamed into a fan, as discussed in the literature for similar problems in other contexts. We have shown that such a model is able to account not only for the magnitude of the pulsations in PSR 1055-52 and PSR 0656+14, but also for the phase shifts between the soft and the hard component observed in all the three objects. We have shown that this same model is able to reproduce also the decrease of the pulsed fraction observed in the soft X-ray emission of Geminga, without requiring any special model for the composition of the surface of this object. However, we have found that the type and composition of the atmosphere plays an important role. An atmosphere made of heavy elements is not able to account for significant $`Pfs`$ (if only a dipolar field is assumed for the thermal component), while a light element atmosphere with opacity $`\nu ^3`$ has a very pronounced tail which is able to account for the all spectrum (up to about 1.5 - 2 keV) for PSR 0656 and PSR 1055, without requiring a second, hotter component. Such a model is not able to account for the sudden increase of the $`Pfs`$ (accompanied by a phase shift) that is observed for these two pulsars around 0.5 keV. In the case of Geminga, the implied size of the hot spots is too small to produce any significant effect. On the other hand, a light element atmosphere with opacity $`\nu ^1`$ requires a second, hotter component to account for the overall spectrum, and, as we have shown, such a model is able to account for the observed pulsed fractions and phase shifts, if one allows the emission from the hot spots to have a variable degree of beaming for the various objects. We have argued that this this is plausible, as it depends on the spatial and energy distribution of the bombarding particles, which is likely to vary among the various objects.

In conclusion, we need to stress that the model we have presented here was aimed at explaining features observed in the X-ray emission from middle aged pulsars from a qualitative point of view. A more quantitative analysis, which includes detector and absorption effects 555For a perfect (diagonal) detector, absorption does not affect the energy-dependent $`Pfs`$ (but it does affect the $`Pf`$ over a finite energy bandwidth; Perna et al. 2000 ); however, for a real detector, absorption can affect also the energy-dependent $`Pfs`$ to a certain extent (Page 1995 ). would not change any of the qualitative features reproduced here. On the other hand, a thorough analysis, inclusive of detailed predictions for the shape of the pulse profiles, would require very detailed models for the atmosphere and for the temperature distribution on the star, as well as much better data with which to compare the model. We anticipate that forthcoming observational data from the Chandra and XMM missions will soon make this endeavor feasible.

We thank Jonathan McDowell for support with use of the XSPEC software.


text:
# The Evolution of Dark-Matter Dominated Cosmological Halos

## Abstract

Adaptive SPH and N-body simulations were carried out to study the evolution of the equilibrium structure of dark matter halos that result from the gravitational instability and fragmentation of cosmological pancakes. Such halos resemble those formed by hierarchical clustering from realistic initial conditions in a CDM universe and, therefore, serve as a test-bed model for studying halo dynamics. The dark matter density profile is close to the universal halo profile identified previously from N-body simulations of structure formation in CDM, with a total mass and concentration parameter which grow linearly with scale factor $`a`$. When gas is included, this concentration parameter is slightly larger than the pure N-body result. We also find that the dark matter velocity distribution is less isotropic and more radial than found by N-body simulations of CDM.

CDM Simulations vs. Observed Halos N-body simulations of structure formation from Gaussian-random-noise density fluctuations in a cold dark matter (CDM) universe have revealed that dark matter halos possess a universal density profile that diverges as $`r^\gamma `$ near the center, with $`1\gamma 2`$ Moore00 NFW97 . There exists a discrepancy between these singular density profiles found in N-body simulations and current observations of the rotation curves of nearby dwarf galaxies Moore00  and of strong lensing of background galaxies by the galaxy cluster CL0024+1654 Tyson98 Shapiro00 , which suggest that dark-matter dominated halos of all scales have flat density cores, instead.

Halo Formation by Pancake Instability and Fragmentation The model we use to examine the formation of dark-matter dominated halos is that of cosmological pancake instability and fragmentation, previously discussed in detail by Valinia97 . Halos formed by such a pancake instability have density profiles very similar to those formed hierarchically in CDM models (e.g. NFW profile NFW97 ), providing a convenient alternative to more complicated simulations with more realistic initial conditions Alvarez00 Martel01 ValPhd . The ASPH/P3M simulations considered here were described by Alvarez00 ; this paper extends that analysis to evolutionary trends in the dark matter halo structure.

Main Results:

* For $`a/a_c`$ between 3 and 7, the halo can be fit by an NFW profile, with mass within $`r_{200}`$ growing linearly with scale factor $`a`$, when simulated either with or without gas: $`M_{200}(x)0.07x`$, where $`xa/a_c`$, and $`a_c`$ is the scale factor at primary pancake collapse (see Figs. 1 & 2). This mass evolution resembles that of self-similar spherical infall Bert85 , despite the anisotropy associated with pancake collapse and filamentation and periodic boundary conditions.
* After $`a/a_c=3`$, the concentration parameter $`c_{NFW}r_{200}/r_s`$, determined by best-fitting an NFW density profile NFW97  to our simulation halos, grows roughly linearly with scale factor $`a`$: $`c_{NFW}(x)1.33x0.18`$ (without gas), $`c_{NFW}(x)1.49x0.37`$ (with gas) (i.e. the linear slope is steeper in the case with gas included). Fluctuations in $`c_{NFW}`$ around this trend are smaller when gas is included. This evolution we find for $`c_{NFW}`$ is reminiscent of that reported for halos in CDM N-body simulations Bullock00 . However, the latter applies to halos of a given mass which are observed at different epochs, and, therefore, reflects the statistical correlation of halo mass with collapse epoch in the CDM model, while our result follows an individual halo.
* The anisotropy parameter $`\beta 1v_t^2/(2v_r^2)`$, where $`v_t`$($`v_r`$) are tangential (radial) velocities, is shown in Figure 2. Pancake halos are somewhat radially biased, with $`\beta 0.6`$, about twice the value reported for halos in CDM N-body simulations ENF98 Thomas98 . With no gas included, the average anisotropy in the halo does not change very much with time, while the inclusion of gas leads to a slight drop after $`a/a_c=5`$.

This work was supported by NASA ATP grants NAG5-7363 and NAG5-7821, NSF grant ASC-9504046, and Texas Advanced Research Program grant 3658-0624-1999.


text:
# Dynamical Initial Conditions in Quantum Cosmology

## Isotropic loop quantum cosmology.

In the triad representation of isotropic loop quantum cosmology  the scale factor $`a^+`$ is replaced by a discrete label $`n`$ which parameterizes eigenvalues of the triad operator. An orthonormal basis of the kinematical Hilbert space is given by quantum states $`|n`$ labeled by the triad eigenvalue $`n`$ which also determines volume eigenvalues: $`\widehat{V}|n=V_{(|n|1)/2}|n`$ with

$$V_j=(\gamma l_\mathrm{P}^2)^{\frac{3}{2}}\sqrt{\frac{1}{27}j(j+\frac{1}{2})(j+1)}$$
(1)
($`\gamma ^+`$, which is of order one, is the Barbero–Immirzi parameter labeling inequivalent representations of the classical Poisson algebra, and $`l_\mathrm{P}=\sqrt{\kappa \mathrm{}}`$ with $`\kappa =8\pi G`$ is the Planck length). The volume operator has eigenvalue zero with threefold degeneracy (for the states $`|\pm 1`$ and $`|0`$), but only one of them, $`|0`$, has degenerate triad and so corresponds to the classically singular state. The wave function $`\psi (a,\varphi )`$ of standard quantum cosmology is replaced by the coefficients $`s_n(\varphi )`$ of a state $`|s=_ns_n|n`$. For large $`|n|`$ the correspondence between $`a`$ and $`n`$ is $`|n(a)|6a^2\gamma ^1l_\mathrm{P}^2`$ which follows from the volume spectrum ($`|n|=2j+1`$).

The Hamiltonian constraint equation for spatially flat models takes the form of a discrete evolution equation (see  for the case of models with positive spatial curvature):

$`{\displaystyle \frac{1}{4}}(1+\gamma ^2)A_n^{(8)}s_{n+8}(\varphi )A_n^{(4)}s_{n+4}(\varphi )2A_n^{(0)}s_n(\varphi )`$ (2)
$`A_n^{(4)}s_{n4}(\varphi )+{\displaystyle \frac{1}{4}}(1+\gamma ^2)A_n^{(8)}s_{n8}(\varphi )`$ (3)
$`={\displaystyle \frac{1}{3}}\gamma \kappa ł_\mathrm{P}^2\widehat{H}_\varphi (n)s_n(\varphi )`$ (4)
with

$`A_n^{(\pm 8)}`$ $`:=`$ $`\left(V_{|n\pm 8|/2}V_{|n\pm 8|/21}\right)k_{n\pm 8}^\pm k_{n\pm 4}^\pm `$ (5)
$`A_n^{(\pm 4)}`$ $`:=`$ $`\left(V_{|n\pm 4|/2}V_{|n\pm 4|/21}\right)`$ (6)
$`A_n^{(0)}`$ $`:=`$ $`\left(V_{|n|/2}V_{|n|/21}\right)`$ (8)
$`\times \left({\displaystyle \frac{1}{8}}(1+\gamma ^2)(k_n^{}k_{n+4}^++k_n^+k_{n4}^{})1\right)`$
where $`V_j`$ are the eigenvalues (1) of the volume operator with $`V_1=0`$, and the coefficients $`k_n^\pm `$ can be chosen to be non-vanishing by a suitable ordering of the extrinsic curvature operator and are approximately $`\mathrm{sgn}(n)`$ for large $`|n|`$ (see  for explicit expressions in terms of the volume eigenvalues). Here we introduced a matter Hamiltonian $`\widehat{H}_\varphi `$ whose particular form is irrelevant. It only matters that it acts diagonally in the triad degrees of freedom which is always the case in the absence of curvature couplings.

## Pre-Classicality.

Compared to the standard Wheeler–DeWitt equation of second order the discrete evolution equation is of order sixteen. So the problem of a unique solution seems to be more severe at first sight, but many of the additional solutions can easily be seen to not correspond to a semiclassical solution. Let us call a wave function $`s_n`$ pre-classical if and only if, at large volume ($`n1`$), it is not strongly varying at the Planck scale (increasing the large label $`n`$ by one), although it may oscillate on much larger scales (increasing $`n`$ by a macroscopic amount). Note that a difference equation with fixed step size, as is always the case here, may have solutions which are very different from those of an approximating differential equation even though all solutions of the differential equation are well approximated when the step size goes to zero . A common possibility is a solution with alternating sign between successive $`n`$, which cannot correspond to a continuous solution of a differential equation. Due to instabilities there can also be solutions with exponentially increasing absolute value, even in regimes where the solutions of the differential equation are purely oscillating (i.e. in the classically allowed range in a WKB approximation). A precise formulation of the phrase “not strongly varying” can be given in the following way. Note first that the Barbero–Immirzi parameter $`\gamma `$ enters $`a(n)=\gamma l_\mathrm{P}^2n/6`$, which is used here as internal time. Although the physical value of $`\gamma `$ is fixed and of order one , we can use the $`\gamma 0`$ limit, together with $`n\mathrm{}`$ such that $`a(n)`$ is fixed, to decide whether a wave function is pre-classical. In this limit the difference $`a(n+1)a(n)`$ becomes infinitesimal implying a continuum limit. A wave function $`s_n`$ is pre-classical if and only if its limit $`\gamma 0`$, $`n\mathrm{}`$ exists, providing a rigorous check of the pre-classicality condition. Note that $`\kappa `$ and $`\mathrm{}`$, and so $`l_\mathrm{P}`$, are fixed in this limit and we are still dealing with quantum cosmology. In fact, standard quantum cosmology can be shown to be the above limit of loop quantum cosmology.

Our condition picks out only those solutions which are oscillatory on large scales but almost constant on the Planck scale. Since this is a pre-requisite for a subsequent WKB-approximation, we call it pre-classicality. Whenever it is fulfilled, a discrete wave function $`s_n(\varphi )`$ can be approximated at large $`n`$ by a standard continuous wave function $`\psi (a):=s_{n(a)}`$ with $`n(a)=6a^2\gamma ^1l_\mathrm{P}^2`$ as above, which approximately solves the standard Wheeler–DeWitt equation up to corrections of order $`\sqrt{\gamma }l_\mathrm{P}/a`$ . Thus, standard quantum cosmology is realized only as an approximation valid at large volume where the discreteness of quantum geometry is irrelevant (see Fig. 1).

Since the Wheeler–DeWitt equation is of second order and so has two independent solutions, there can be at most two independent pre-classical solutions $`s^\pm `$ of our discrete evolution equations, such that any pre-classical solution can be written as $`s=as^++bs^{}`$ with $`a,b`$.

To demonstrate this explicitly, we introduce

$`t_m`$ $`:=`$ $`\gamma ^1l_\mathrm{P}^2(V_{2|m|}V_{2|m|1})s_{4m}`$
$`P(m)`$ $`:=`$ $`{\displaystyle \frac{1}{3}}\gamma \kappa l_\mathrm{P}^2H_\varphi (m)(V_{2|m|}V_{2|m|1})^1,`$
using the expectation value $`H_\varphi (n)`$ of $`\widehat{H}_\varphi (n)`$ in a matter state, such that for $`|n|1`$, where $`k_n^\pm \mathrm{sgn}(n)`$, the evolution equation (2) takes the form

$`{\displaystyle \frac{1}{4}}(1+\gamma ^2)t_{m+2}t_{m+1}+({\displaystyle \frac{1}{2}}(3\gamma ^2)+P(m))t_m`$ (9)
$`t_{m1}+{\displaystyle \frac{1}{4}}(1+\gamma ^2)t_{m2}=0.`$ (10)
In a classical regime $`|m|`$ is large and $`P(m)\frac{2}{3}\kappa H_\varphi /a`$ is approximately constant on a range small compared to $`|m|`$. In this case we have a linear difference equation with constant coefficients whose solutions can be found by an ansatz $`t_me^{im\theta }`$ with $`\theta `$ which in (10) yields the quadratic equation

$`(1+\gamma ^2)\mathrm{cos}^2\theta 2\mathrm{cos}\theta +1\gamma ^2+P=0`$

which has solutions

$`\mathrm{cos}\theta =(1+\gamma ^2)^1\left(1\pm \sqrt{\gamma ^4(1+\gamma ^2)P}\right)`$

being real with modulus smaller than one such that $`\theta `$ is real when $`\gamma `$ is of the order one and $`P`$ is small.

If the matter does not contribute a Planck size energy, $`P`$ is small and we have $`\mathrm{cos}\theta _0=1ϵ+O(ϵ^2)`$ or $`\mathrm{cos}\theta _1=(1+\gamma ^2)^1(1\gamma ^2)+ϵ+O(ϵ^2)`$ with $`0<ϵ:=\frac{1}{2}\gamma ^2P1`$. The first possibility, expanding $`\mathrm{cos}\theta _0=1\frac{1}{2}\theta _0^2+O(\theta _0^4)`$, leads to two solutions $`\theta _0=\gamma \sqrt{P}+O(P)`$ and $`\theta _0`$ with $`|\theta _0|1`$ both of which imply pre-classical $`t_m^\pm =e^{\pm im\theta _0}`$. Because $`\gamma `$ is not large compared to one, the second possibility $`\mathrm{cos}\theta _1`$ leads to $`\theta _1`$ which violates pre-classicality (e.g., for $`\gamma =1`$ we have $`\theta _1=\pm \pi /2`$ and $`t_m(\pm i)^m`$).

All 16 independent solutions $`t_{n/4}=\gamma ^1l_\mathrm{P}^2(V_{|n|/2}V_{|n|/21})s_n`$ of (2) can be obtained as

$`t_{n/4}`$ $`=`$ $`e^{\pm in\theta _0/4},e^{\pm in\theta _1/4},`$
$`(1)^ne^{\pm in\theta _0/4},(1)^ne^{\pm in\theta _1/4},`$
$`\sigma ^ne^{\pm in\theta _0/4}\text{ or }\sigma ^ne^{\pm in\theta _1/4}`$
where $`\sigma `$ can be $`+i`$ or $`i`$. Obviously, only the first two are pre-classical. The definition using $`\gamma 0`$ ($`a`$ finite) is applied as follows: with $`n=6a^2\gamma ^1l_\mathrm{P}^2`$ we have

$`\underset{\gamma 0}{lim}\gamma ^1l_\mathrm{P}^2(V_{|n|/2}V_{|n|/21})=a/2,`$
$`\underset{\gamma 0}{lim}e^{\pm in\theta _0/4}=\mathrm{exp}\left(\pm {\displaystyle \frac{3}{2}}i\sqrt{P}a^2/l_\mathrm{P}^2\right),`$
whereas $`lim_{\gamma 0}\theta _1=\pi `$ and so the limit for $`\gamma 0`$, $`n\mathrm{}`$ does not exist for the remaining 14 solutions.

Thus, by using the pre-classicality requirement, which is a prerequisite for any semiclassical analysis, we arrive at the same situation as in standard quantum cosmology: there are two independent solutions from which we have to select a linear combination up to norm.

## Dynamical Initial Conditions.

Up to now we considered only the semiclassical regime, but there is an additional feature of loop quantum cosmology which emerges right at the classical singularity, deeply in the Planck regime where the approximation by standard quantum cosmology breaks down: the highest order (or lowest order when we evolve backwards) coefficient vanishes when we try to determine $`s_0`$. At first sight, it seems that this is a breakdown of the evolution similar to the classical situation. However, as demonstrated in , this is not the case in the particular factor ordering of the constraint chosen above because $`s_0`$ completely drops out of the evolution equation. (This observation depends crucially on the fact that $`\widehat{H}_\varphi s_0(\varphi )=0`$ which is always true in quantum geometry  but would be impossible without space-time discreteness.)

Instead of determining $`s_0`$ the evolution equation leads to a consistency condition for the initial data: starting from a general pre-classical solution $`s_n=as_n^++bs_n^{}`$ for large $`n`$ and evolving backwards, we eventually arrive at a point where we have to apply (2) for $`n=8`$. At this value of $`n`$ the lowest order coefficient $`A_8^{(8)}`$ vanishes as noted above causing $`s_0`$ to drop out. Since by assumption we have already determined all $`s_n`$ for $`n>0`$ (for which there is no vanishing coefficient in the evolution equation), the would-be equation for $`s_0`$ leads to a further condition for higher $`s_n`$ ($`s_4`$, $`s_8`$, $`s_{12}`$ and $`s_{16}`$) which upon inserting the general pre-classical solution $`s_n=as_n^++bs_n^{}`$ implicitly yields a linear relation between the two free parameters $`a`$ and $`b`$. This leaves us with a unique solution (up to norm).

## Conclusions.

We have shown that loop quantum cosmology implies a discrete evolution equation which uniquely determines a state (up to norm) behaving semiclassically at large volume. It is important to adapt the standard condition for semiclassicality in a WKB approximation taking the discreteness of time into account. This leads already to a strong reduction of the allowed solutions, but the crucial condition for the uniqueness arises only from the particular structure of the evolution equation in quantum geometry. We remark that in general it is only possible to require pre-classicality at one connected domain of large volume. If one evolves through a classical singularity, the wave function may pick up components which oscillate at the Planck scale (see Fig. 1). The precise form of these oscillations depends on factor ordering ambiguities (in the coefficients $`k_n^\pm `$ entering the constraint) and the use of the Lorentzian (versus Euclidean) theory.

Such a unique wave function generally differs from those obtained with boundary proposals of standard quantum cosmology. By choosing a real prefactor it is always real (for flat spatial slices the evolution equation has real coefficients; this no longer holds true for spatially curved models) and so cannot coincide with the “tunneling” wave function . While the “no-boundary” proposal  also leads to a real wave function, it is imposed on the standard Wheeler–DeWitt equation at the Planck scale where large deviations to loop quantum cosmology occur. Thus, in general its wave function of a universe is different from the unique pre-classical solution found here. The consistency condition for the initial data in loop quantum cosmology may be expressed as $`s_0=0`$ which is reminiscent of DeWitt’s $`\psi (0)=0`$  (to achieve this, an ad hoc Planck potential has been introduced in ). However, since these two conditions are imposed on completely different evolution equations, the selected solutions in general differ. As Fig. 1 shows, there may be a good coincidence in certain models, but only if the curvature is small at all times which can happen only in the absence of matter.

Contrary to all other proposals for boundary conditions in quantum cosmology, our dynamical initial conditions are not chosen to fulfill an a priori intuition about the “creation” of a universe but derived from the evolution equation which, in turn, is derived from quantum geometry, a candidate for a complete theory of quantum gravity. Therefore, one equation provides both the dynamical law and initial conditions. As we have seen, the critical condition, which crucially depends on quantum geometry, emerges from evaluating the evolution equation at the state which corresponds to the classical singularity. So in contrast to the classical situation where a singularity leads to unpredictability, in quantum geometry the regime of the classical singularity fixes ambiguities in the wave function of a universe.

## Acknowledgements.

The author is grateful to A. Ashtekar for discussions and a careful reading of the manuscript. This work was supported in part by NSF grant PHY00-90091 and the Eberly research funds of Penn State.


text:
# Experimental Constraints on the Axion Dark Matter Halo Density

## 1 Introduction

Measurements such as the BOOMERANG and MAXIMA cosmic microwave background radiation power spectra (de Bernardis et al., 2000; Hanany et al., 2000) and Type Ia supernovae data (Perlmutter et al., 1999; Reiss et al., 1998) suggest that dark matter makes up the bulk of the matter content in the universe. The agreement between expectations from Big Bang nucleosynthesis (Burles et al., 1999) and the primordial abundance of the lightest elements (Olive et al., 2000) strongly constrains the total baryonic content to be a small value, implying that the majority of the dark matter must be non-baryonic. The results of other diverse astrophysical measurements based on the Sunyaev-Zel’dovich effect (Grego et al., 2001), strong gravitational lensing (Cheng and Krauss, 2000) and galactic flows (Zehavi and Dekel, 1999) lend powerful support to this picture.

Although the nature of dark matter remains unknown, its gravitational effect is pronounced: much of the dynamics of spiral galaxies cannot be understood without there being a massive dark matter halo (Sofue and Rubin, 2000). Two categories of particle cold dark matter (CDM) candidates have survived experimental and theoretical scrutiny over time: the lightest supersymmetric particle (Goldberg, 1983; Ellis et al., 1984) and the axion (Weinberg, 1978; Wilczek, 1978). Neutrinos and other forms of hot dark matter are thought to contribute insignificantly to closure density and, in any case, cannot explain structure formation. Galaxy formation requires CDM, i.e., dark matter which is already non-relativistic at the time of decoupling. Recently, experiments have begun with the sensitivity either to detect or exclude possible CDM halo candidates. In this paper we present upper limits on the local axion halo density derived from a search for cold dark matter axions.

## 2 Axion Physics

The axion is the pseudo Nambu-Goldstone boson (Weinberg, 1978; Wilczek, 1978) associated with a new spontaneously broken global UPQ(1) symmetry invented to suppress strong CP-violation (Peccei and Quinn, 1977). There is some model dependence in assigning UPQ(1) charges to particles: in the KSVZ scheme (Kim, 1979; Shifman et al., 1980) the axion only couples to quarks at tree level, while in the GUT-inspired DFSZ model (Zhitnitsky, 1980; Dine et al., 1981) it couples to both quarks and leptons. The axion acquires a mass that scales inversely with the (unknown) energy scale fa at which the UPQ(1) symmetry breaking occurs. Initially fa was presumed to be the electroweak energy scale, but such massive axions were quickly ruled out in, e.g., beam dump experiments (Asano et al., 1981). Subsequently, it was proposed that axions possess such small couplings to matter and radiation that for all practical purposes they would remain forever “invisible”. Shortly thereafter, an experiment was proposed that could make even very light axions detectable (Sikivie, 1983). The U.S. axion search experiment is predicated on this approach whereby the axion converts into a single photon via the inverse Primakoff effect. We use a resonant cavity permeated by a strong static magnetic field, where the large number density of virtual photons from the field enhances axion decay.

The allowed axion mass is constrained to 10-2 to 10-6 eV. Axions with a mass less than a few eV would have cooled the core of supernova 1987a to such an extent that the distribution of neutrino arrival times would be inconsistent with observation (Turner, 1988). Even heavier axions have been ruled out by a variety of astrophysical and terrestrial searches (Turner, 1990; Raffelt, 1990; Ressell, 1991; Gnedin, 1999). Conversely, if the axion mass is less than some value they would have been overproduced in the early universe. This lower mass limit has been calculated for various axion production mechanisms under different early-universe scenarios, e.g., “vacuum realignment” (Preskill et al., 1983; Abbott and P. Sikivie, 1983; Dine and Fischler, 1983), “string decay” (Battye and Shellard, 1994; Yamaguchi et al., 1999; Hagmann et al., 2001) and “wall decay” (Chang et al., 1999). The vacuum misalignment mechanism provides a lower mass limit of $``$ 10-6 eV which we adopt for our search strategy; the other mechanisms produce a value that is in fairly close agreement. Common to all of these mechanisms is the misalignment of the axion field with respect to the CP-conserving (minimum energy) direction when the axion mass turns on during the QCD phase transition. Axions produced in this way are very cold. Their typical momentum was of order the inverse of the horizon scale at the QCD phase transition ($`p_a10^8`$ eV) when the temperature was of order 1 GeV. The cosmological energy density in these cold axions is of order: (Kolb and Turner, 1990)

$$\mathrm{\Omega }_a0.5(\frac{\mu \mathrm{eV}}{m_a})^{7/6}.$$
(1)
Hence, if $`m_a`$ is of order a few $`\mu `$eV, the mass range where we search, axions contribute significantly to the energy density of the universe. Studies of large scale structure formation support the view that the dominant fraction of matter is in the form of CDM. Since CDM necessarily contributes to galactic halos by falling into the gravitational wells of galaxies (halo axions in our galaxy possess a virial velocity $``$ 10-3c), there is excellent motivation to search for axions as constituents of our galactic halo.

There is a substantial body of evidence that our own galaxy is surrounded by a massive dark halo, though its exact properties are not well-constrained (Zaritsky, 1998; Alcock, 2000). Of particular interest to this paper is the local dark matter halo density, whose value depends on the degree of halo flattening as well as the core radii of the various dark matter components. To derive a reliable mass density, one can turn to parameterizations of the density distribution, rejecting distributions which fail to match observational constraints, such as reproducing the local rotation speed of 200-240 kmsec-1 (Gates et al., 1995). A key element of this approach is the use of microlensing data to estimate the fraction of local dark matter that is in the form of compact objects. Employing this methodology, one arrives at a halo density of 9.2$`{}_{3.1}{}^{}{}_{}{}^{+3.8}`$$`\times `$10-25gm/cm3. If massive compact halo objects (MACHOs) comprise a negligible fraction of the local halo density, then the above number is likely an underestimate. The local halo density may also be enhanced because of our proximity to a possible dark matter caustic (Sikivie, 1998). Our experimental analysis directly constrains the local density of the axionic component of the halo, and as such is independent of astronomical observations and assumptions.

## 3 Experimental Technique

The interaction between axions and photons can be written as

$$=g_{a\gamma \gamma }a\stackrel{}{E}\stackrel{}{B},$$
(2)
where $`g_{a\gamma \gamma }`$ is the relevant coupling, $`a`$ the axion field and $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ the electric and magnetic fields, respectively. Since $`g_{a\gamma \gamma }`$ is very small in the mass range of interest, the spontaneous decay lifetime of an axion to two real photons is vastly greater than the age of the universe. In our experiment, located at Lawrence Livermore National Laboratory, a high-$`Q`$ resonant cavity and superconducting magnet stimulate axion conversion into a single real photon. Resonant conversion occurs when the cavity resonant frequency equals the axion rest mass. Because this mass is, a priori, unknown, resonant frequencies are changed by moving either ceramic or metallic tuning rods from the wall to the center of the cavity. For a resonant cavity with a loaded quality factor $`Q_L`$, the axion-to-photon conversion power is

$$P=410^{26}\text{W}\left(\frac{V}{0.22\text{m}^3}\right)\left(\frac{B_0}{8.5\text{ T}}\right)^2C_{nl}\left(\frac{g_\gamma }{0.97}\right)^2\left(\frac{\rho _a}{\frac{1}{2}10^{24}\frac{g}{\text{cm}^3}}\right)\left(\frac{m_a}{2\pi (\text{GHz})}\right)\text{min}(Q_L,Q_a),$$

where $`V`$ is the cavity volume, $`B_0`$ the magnetic field strength, $`C_{nl}`$ the mode-dependent cavity form factor, $`g_\gamma `$ the reduced coupling constant (equal to $`g_{a\gamma \gamma }\pi f_a/\alpha `$), $`\rho _a`$ the axion halo density and min$`(Q_L,Q_a)`$ the smaller of either the cavity or axion quality factors. Typical values for the first four parameters are 0.2 $`m^3`$, 7.5 $`T`$, 0.6 and 0.97, respectively. The copper cavity has a loaded (critically coupled) $`Q_L`$ $``$ 105, whereas $`Q_a`$, the ratio of the energy to the energy dispersion of the axion, is a factor of ten or so larger over the present frequency range. The total power that results from Eq. 3 is of order $``$ $`10^{22}`$ W; our cavity and amplifiers are cooled to a few degrees Kelvin to minimize thermal noise. Fig. 1 is a schematic of the axion receiver chain. A microwave signal centered at the cavity resonant frequency and approximately 30 kHz wide is coupled out of the cavity by an electric field probe and subsequently mixed down (in two stages) to near audio frequencies. At any given frequency 104 spectra are averaged by fast-Fourier-transform (FFT) hardware, with each spectrum sampled for 8 ms. The corresponding Nyquist resolution of 125 Hz is well matched to the width ($``$ 750 Hz) of axions thermalized by interactions with the galactic gravitational potential.

The Dicke radiometer equation (Dicke et al., 1946) dictates the integration time necessary to achieve a specified signal-to-noise ratio (SNR)

$$\mathrm{SNR}=\frac{\mathrm{P}}{\mathrm{T}_\mathrm{n}}\sqrt{\frac{\mathrm{t}}{\mathrm{B}}},$$
(3)
where $`P`$ is the axion power from Eq. 3 (times a factor which accounts for the external coupling), $`T_n`$ is the noise temperature, $`t`$ is the integration time and $`B`$ is the axion bandwidth (defined as $`f/Q_a`$). Since the desired SNR is not attained in a single pass over a given frequency interval, data in a single frequency bin are the result of combined data from numerous overlapping spectra. The SNR for this multiple-pass data in the interval 550 $``$ f $``$ 810 MHz is shown in Fig. 2. A second data set was formed by co-adding six neighboring bins into single 750 Hz bins suitable for the virialized axion signal. To determine the number of candidate peaks that must be rescanned to obtain an overall confidence limit of $``$90%, artificial peaks are injected into the data via software. As the cut threshold is lowered, the number of candidate peaks increases rapidly. A typical cut threshold of 2.3 $`\sigma `$ (6-bin) applied to our data yields numerous candidates which are rescanned to a SNR commensurate with the original data. These data are subsequently added to the original data and from these combined data sets a reduced set of candidates is generated and scanned at the corresponding frequencies. A final round of data-combining produces a persistent-candidate list. Candidates above a threshold of 3.5 $`\sigma `$ in these data are manually inspected. A detailed description of the experiment and analyses may be found in (Peng et al., 2000; Asztalos et al., 2001).

## 4 Results

We have examined these data for candidates in each of 2.08$`\times `$106 125 Hz and 750 Hz bins in the region 550 $``$ f $``$ 810 MHz. A total of 13712, 1369 and 34 candidates survived each stage of 6-bin data cuts, respectively. All 34 persistent candidates have been identified with strong external radio peaks. To derive an upper limit on the axion contribution to the local halo density, we fix the axion-to-photon coupling $`g_{a\gamma \gamma }`$ at the KSVZ level and invert Eq. 3 to calculate $`\rho _a`$ as a function power deposited in the cavity and axion mass. The absence of a persistent signal in these data over this range permits us to impose the limits shown in Fig. 3, where we plot the excluded axion dark matter halo densities for both KSVZ (lower curve) and DFSZ (upper curve) axions as a function of axion mass and frequency over the interval 550 $``$ f $``$ 810 MHz. The small variations in these density limits represent effective integration times somewhat longer or shorter than that prescribed by Eq. 3. The nominal excluded mass density lies near 0.45 GeV/cm3 for KSVZ axions and 3.0 GeV/cm3 for DFSZ axions. The former is comparable to the best estimate of the local dark matter halo density.

## 5 Conclusions

There is abundant evidence that our own galaxy, like other spiral galaxies, contains a vast dark matter halo. Observation can neither differentiate the various candidates, nor well constrain other parameters that describe the halo, e.g., the local dark matter density. Since 1995 we have been using a single resonant cavity to search for axions which may constitute the local dark matter halo over the frequency interval 550 $``$ f $``$ 810 MHz. The lack of a persistent signal allows us to exclude the axion from contributing more than 0.45 GeV/cm3 to the halo dark matter mass density over the mass range of 2.3$`\times `$10-6 $``$ $`m_a`$ $``$ 3.4$`\times `$10-6 eV, should axions couple only to hadrons according to the KSVZ prescription, with 90% confidence. This restriction is relaxed to around 3.0 GeV/cm3 in the DFSZ model, also with 90% confidence. It should be noted that other KSVZ- and DFSZ-like implementations exist. Some of these models (including some DFSZ-like models) give rise to coupling constants that are larger than the benchmark KSVZ $`g_{a\gamma \gamma }`$ used in this paper (Kim, 1998). These, too, are ruled out by our results over the mass range quoted above.

## 6 Acknowledgments

This work was supported by the U.S. Department of Energy under Contract Nos. DE-FC02-94ER40818, W-7405- ENG-48, DE-FG02-97ER41029, DE-AC02-76CH0300, DE-FG02-90ER40560, DE-AC-03-76SF00098, and the National Science Foundation under Award No. PHY-9501959.


text:
# Comment on “Requirement of optical coherence for continuous–variable quantum teleportation” by Terry Rudolph and Barry C. Sanders

## Abstract

The argument of Rudolph and Sanders, while technically correct, raises conceptual problems. In particular, if carried to its logical conclusion, it would disallow the use in our theories of any time $`t`$ with implied resolution beyond that of direct human experience.

In a recent paper , Rudolph and Sanders have argued that, contrary to Ref. , continuous–variable quantum teleportation has not been and, in fact, cannot be, achieved using a laser as a source of coherent radiation. They base their argument on the fact that a laser is not a source of coherent radiation, in the sense that the output of a laser is not a coherent state, but a mixture of coherent states over all possible phases. Although the formal analysis of Rudolph and Sanders is indisputably correct, there are deep conceptual issues raised by this analysis that they seem not to have considered.

Having established, following Mølmer , the lack of absolute phase of a laser beam, Rudolph and Sanders say that they are not asserting that the production of coherent states of light is impossible, and that “basic quantum electrodynamics shows that a classical oscillating current can produce coherent states.” At first sight this seems unassailable, but let us examine it more closely.

First, how would one obtain a current oscillating at optical frequencies? The natural oscillators at optical frequencies are the electrons in atoms and molecules. But how can one create coherent excitations in such oscillators if one cannot start with coherent light? One answer would be to “strike” the atom (with a free electron, for example), to set it “ringing”. Even assuming that the dynamics of the collision can be fully determined, to produce a coherent oscillation, the time of the collision would have to be known to an accuracy less than an optical cycle, of order $`10^{15}`$s. Otherwise one would have to average over all possible phases and one would be left with exactly the same problem as with the laser.

Let us say for arguments sake that it is possible to know the time of collision to an accuracy of $`10^{15}`$s. The question is, with respect to what? What clock ticks $`10^{15}`$ times per second? If we ignore that problem and allow Alice and Bob such clocks, ticking in phase, then surely this solves the problem? Well not really, because how do we know that the clock really has a definite phase? How do we know that, relative to the absolute time of the universe, Alice’s clock does not have a random phase? Certainly Alice cannot simply look at her clock and verify that it does not have a random phase, because she cannot perceive anything in $`10^{15}`$s.

The point is that the whole idea of an absolute time standard for the universe is highly questionable, even ignoring any issues to do with relativity. As conscious beings we feel that we experience time directly, but experiments show that the limit of our time resolution is in the range of tens or even hundreds of milliseconds . On this basis, it is impossible to establish the absolute phase of any oscillator of frequency greater than a few tens of Hertz. At higher frequencies we can only talk about the phase of one oscillator relative to another oscillator. This conclusion is not altered by oscillations obtained by frequency $`2^n`$-upling, because the timing of the zeros of the highest harmonic can be no more accurately defined than that of the fundamental.

It follows then that, if we accept the arguments of Rudolph and Sanders, we must conclude that it is impossible to teleport the state of a high-frequency oscillator by any means. Indeed, we should conclude that the state of any high frequency oscillator is always mixed (with all phases equally weighted) so there is probably little point in trying to teleport it.

While this stance is a logically consistent one, it is far more useful to acknowledge that there is no absolute time standard. All we can ever do, for experiments involving time resolution beyond direct human experience, is to use an agreed time standard. In this context, a laser field “ticking” at $`10^{15}`$ Hz is as good a “clock” as anything. There are no better clocks, even in principle.

To conclude, treating a laser “clock” as if it had a fixed phase relative to some absolute standard (that is, ascribing to it a coherent state) may be a “convenient fiction” . It may even be committing the “partition ensemble fallacy” . However, the alternative, if carried to its logical conclusion, would be never to write down a time $`t`$ or a phase $`\varphi `$ in our theories if its implied resolution would be beyond that of direct human experience. To scientists and engineers, this would be unacceptable pedantry.