text:
# 1 Introduction

## 1 Introduction

The density functional approach is one of the most powerful and widely applicable approaches to nonuniform fluids . Its idea is to express the free energy as a functional of locally varying one-particle densities $`\{\rho (𝐫)\}`$, with an ideal gas contribution and an excess free energy which accounts for the interparticle interactions. This allows one to calculate the structure and properties of fluids with various kinds of inhomogeneities. Density functionals for simple fluids have reached a considerable degree of sophistication and are able to describe such complex phenomena as freezing, wetting, and surface melting .

One of the most established density functionals for hard sphere fluids is the Rosenfeld functional . It reproduces by construction the Percus-Yevick solution for the correlation function, which is known to be very good . As a consequence, it provides highly accurate predictions for the structure of inhomogeneous hard sphere fluids . It has been used successfully to study mixtures of hard spheres  and polydisperse fluids . After a simple anisotropic mapping procedure it can be applied to fluids of fully aligned ellipsoids . Exploiting the Gauss-Bonnet theorem, it has been generalized for isotropic molecular fluids . This generalized functional has been used as a starting point to calculate direct correlation functions of isotropic multicomponent fluids . A Rosenfeld type approach has been developed for mixtures of rods and needles, assuming that the needles are too thin to interact with each other directly , and for systems containing only one single rod immersed in a fluid of spheres . However, to our knowledge, the Rosenfeld functional has not yet been extended to general anisotropic fluids.

The simplest density functional for anisotropic particles is the Onsager functional , which truncates the virial expansion after the leading coefficient. Onsager showed that this functional produces a nematic-isotropic phase transition  in fluids of sufficiently elongated particles. Its predictions are in good agreement with experimental results on systems of tobacco mosaic viruses . One can even show that it describes the transition *exactly* in homogeneous systems of infinitely elongated rods . Thus any functional for nematic liquid crystals should reduce to the Onsager functional in this limit.

In the present paper, we propose a density functional for hard anisotropic particles which interpolates between the Rosenfeld functional and the Onsager functional. It reduces to the Rosenfeld functional in the case of particles with spherical symmetry, and to the Onsager functional for homogeneous fluids of infinitely elongated particles at density close to zero. As a first test of the functional, we have calculated the nematic-isotropic transition for hard spherocylinders and hard ellipsoids, and obtained reasonable results. We believe that our functional might provide a useful new approach to the study of inhomogeneous liquid crystals, e.g. the study of interfacial phenomena such as surface anchoring.

## 2 Background

We consider a fluid of hard anisotropic particles with positions $`𝐫`$ and orientations $`\varpi `$. The grand canonical free energy is the minimum of the free energy functional

$$\beta \mathrm{\Omega }[\rho ]=\beta ^{id}[\rho ]+\beta ^{ex}[\rho ]+𝑑𝐫𝑑\varpi \rho (𝐫,\varpi )[\mu V(𝐫,\varpi )]$$
(1)
with respect to the one-particle density $`\rho (𝐫,\varpi )`$, where $`\beta =\frac{1}{k_BT}`$ is the Boltzmann factor with the temperature $`T`$, $`\mu `$ the chemical potential and $`V(𝐫,\varpi )`$ summarises all external potentials. The first term describes the ideal gas contribution

$$\beta ^{id}=𝑑𝐫𝑑\varpi \rho (𝐫,\varpi )\left(\mathrm{ln}[\rho (𝐫,\varpi )\lambda _T^3]1\right)$$
(2)
with the de Broglie wavelength $`\lambda _T`$. The second term $`^{ex}`$ is the excess free energy, the central quantity in density functional theories .

In the Onsager theory, the excess free energy of two hard particles is given by

$$\beta _O^{ex}=\frac{1}{2}𝑑𝐫𝑑\varpi 𝑑𝐫^{}𝑑\varpi ^{}\rho (𝐫,\varpi )\rho (𝐫^{},\varpi ^{})\text{M}(𝐫𝐫^{},\varpi ,\varpi ^{}).$$
(3)
The Mayer-function M takes the value $`(1)`$ if two particles at $`(𝐫,\varpi )`$ and $`(𝐫^{},\varpi ^{})`$ overlap and vanishes otherwise. In the homogeneous case, the one-particle density $`\rho `$ does not depend on $`𝐫`$ and can be written as $`\rho (𝐫,\varpi )=\rho f(\varpi )`$. The expression 3 then reduces to

$$\frac{\beta _O^{ex}}{N}=\frac{\rho }{2}𝑑\varpi 𝑑\varpi ^{}f(\varpi )f(\varpi ^{})\nu (\varpi ,\varpi ^{}),$$
(4)
where $`\nu (\varpi ,\varpi ^{})`$ is the covolume of two particles with orientations $`\varpi `$ and $`\varpi ^{}`$, and $`N`$ the total number of particles. Eqn. (4) corresponds to a virial expansion up to second order. As mentioned in the introduction, the Onsager theory is exact in the limit of infinitely elongated particles . Compared to real thermotropic liquid crystals, it tends to overestimate the shape anisotropy required to observe a nematic phase at a given finite density.

A natural extension of the Onsager model would be to include higher order terms in the virial expansion. On principle, this is feasible, and an extension of the Onsager functional up to third order was actually carried out for ellipsoids . However , the calculations are very cumbersome due to the increasing complexity of the integrals. As long as one is interested in homogeneous systems, one way to overcome the problem is the *decoupling approximation.* It consists of a resummation of the virial expansion, where the first virial coefficients are calculated exactly, and the remaining ones are approximated by a mapping onto a reference system. Given a virial expansion of the form:

$$\frac{\beta ^{ex}}{N}=\underset{l=1}{\overset{\mathrm{}}{}}\frac{\rho ^l}{l}B_{l+1}\left[f\left(\varpi \right)\right],$$
(5)
the decoupling approximation reads -:

$$\frac{\beta ^{ex}}{N}=\underset{l=1}{\overset{r}{}}\frac{\rho ^l}{l}B_{l+1}\left[f\left(\varpi \right)\right]+\left[\frac{\beta F^{ex,ref}}{N}\underset{l=1}{\overset{r}{}}\frac{\rho ^l}{l}B_{l+1}^{ref}\right]\frac{B_{r+1}\left[f\left(\varpi \right)\right]}{B_{r+1}^{ref}},$$
(6)
where $`B_{l+1}\left[f\left(\varpi \right)\right]`$ are the virial coefficients of the anisotropic fluids, calculated exactly up to the $`r+1`$ order and ref  refers to the reference system. Two possible reference systems have been proposed: the hard sphere fluid and the isotropic fluid of the particles under investigation. If $`r=1`$ and the reference system is the hard sphere fluid with the Carnahan-Starling equation of state , we recover the Parsons-Lee functional

$$\beta _{PL}^{ex}=\frac{1}{8}\frac{(43\eta )}{(1\eta )^2}𝑑𝐫𝑑\varpi 𝑑𝐫^{}𝑑\varpi ^{}\rho (𝐫,\varpi )\rho (𝐫^{},\varpi ^{})\text{M}(𝐫𝐫^{},\varpi ,\varpi ^{}),$$
(7)
where $`\eta `$ is the packing fraction. The Parsons-Lee functional predicts accurately the location of the nematic-isotropic phase transition for a wide range of shape anisotropies . It can also be applied to inhomogeneous systems. However, it does not describe the microscopic structure very well: it yields very crude correlation functions  and produces rather unrealistic density profiles in inhomogeneous systems as a consequence .

The success of the Parsons-Lee theory in predicting the isotropic-nematic coexistence densities motivated Poniewierski and Holyst  and Somoza and Tarazona  to combine the decoupling approximation with a “weighted density approximation” (WDA) scheme . Both are constructed such that one automatically recovers the Onsager functional in the low density limit. In the hard sphere limit, Somoza’s and Tarazona’s  and related  versions reduce to the well-known Tarazona functional . The latter incorporates the Carnahan-Starling equation of state and the Percus-Yevick direct correlation function for homogeneous hard sphere fluids, and has been used with great success to study inhomogeneous hard sphere fluids and solids. Applied to systems of spherocylinders, the extensions  of the Somoza functional generate a nematic phase and several smectic phases. Graf and Löwen have put forward a simplified “modified weighted density approximation” (MWDA), and used it to reproduce complex phase diagrams of spherocylinder fluids .

Besides the Tarazona functional, another equally successful hard sphere functional has established itself in recent years: the fundamental measure theory or Rosenfeld functional . It has the advantage of being based on somewhat more fundamental considerations: it does not require the explicit input of the equation of state and the direct correlation functions, instead they pop out automatically. One obtains by construction the Percus-Yevick solution. Moreover, the functional is formulated a priori for general convex particles of arbitrary geometry. It thus seems to call for a generalization to anisotropic fluids. Unfortunately, there is one problem: the original Rosenfeld functional does not reproduce the Onsager functional in the low density limit.

Here, we propose a simple modification of the Rosenfeld functional which solves that problem for anisotropic particles, and reduces to the original Rosenfeld functional for isotropic particles. Like Somoza’s and Tarazona’s functional, our functional interpolates between an established hard sphere functional and the Onsager functional. Thus we hope that it will be equally useful.

Before introducing our approach, we sketch briefly the basic equations of the original Rosenfeld functional . It can be formulated for general multicomponent fluids of convex hard particles. In our case, different components $`i`$ may be identified with different particle orientations. The functional then reads:

$$\beta _R^{ex}=𝑑𝐫\left(\mathrm{\Phi }_1(𝐫)+\mathrm{\Phi }_2(𝐫)+\mathrm{\Phi }_3(𝐫)\right),$$
(8)
$$\mathrm{\Phi }_1(𝐫)=n_0(𝐫)\mathrm{ln}(1n_3(𝐫))$$
(9)
$$\mathrm{\Phi }_2(𝐫)=\frac{n_1(𝐫)n_2(𝐫)𝐧_1(𝐫)𝐧_2(𝐫)}{1n_3(𝐫)}$$
(10)
$$\mathrm{\Phi }_3(𝐫)=\frac{\frac{1}{3}n_2^3(𝐫)n_2(𝐫)(𝐧_2(𝐫)𝐧_2(𝐫))}{8\pi (1n_3(𝐫))^2}\lambda $$
(11)
The $`\mathrm{\Phi }`$s depend on the weighted densities

$$n_\alpha (𝐫)=\underset{i}{}𝑑𝐫^{^{}}\rho _i(𝐫^{^{}})w_i^{(\alpha )}(𝐫𝐫^{^{}}),$$
(12)
with the number density of the $`i`$th component $`\rho _i(𝐫)`$ and the weight functions

$$\begin{array}{cccccc}\hfill w_i^{(0)}(𝐫)& =& C_i(\widehat{𝐫})\delta (R_i(\widehat{𝐫})r)K_i(\widehat{𝐫})/4\pi \hfill & \hfill w_i^{(2)}(𝐫)& =& C_i(\widehat{𝐫})\delta (R_i(\widehat{𝐫})r)\hfill \\ \hfill w_i^{(1)}(𝐫)& =& C_i(\widehat{𝐫})\delta (R_i(\widehat{𝐫})r)H_i(\widehat{𝐫})/4\pi \hfill & \hfill w_i^{(3)}(𝐫)& =& \mathrm{\Theta }(R_i(\widehat{𝐫})r)\hfill \\ \hfill 𝐰_i^{(1)}(𝐫)& =& C_i(\widehat{𝐫})\delta (R_i(\widehat{𝐫})r)\widehat{𝐧}_i(\widehat{𝐫})H_i(\widehat{𝐫})/4\pi \hfill & \hfill 𝐰_i^{(2)}(𝐫)& =& C_i(\widehat{𝐫})\delta (R_i(\widehat{𝐫})r)\widehat{𝐧}_i(\widehat{𝐫}).\hfill \end{array}$$
(13)
Here $`\widehat{𝐫}`$ denotes the unit vector in the direction of $`𝐫`$, $`R_i(\widehat{𝐫})`$ the radius from the centre of a particle of type $`i`$ to the surface in the direction $`\widehat{𝐫}`$, $`H_i(\widehat{𝐫})`$ and $`K_i(\widehat{𝐫})`$ the local mean and Gaussian curvature on that specific surface point , $`\widehat{𝐧}_i(\widehat{𝐫})`$ the outward unit normal (see Figure 1), and $`\delta `$ and $`\theta `$ are the usual delta and step function. The factor $`C_i(\widehat{𝐫})`$ ensures that integrals over the weight functions $`w_i^{(\alpha )}(𝐫)`$ and $`𝐰_i^{(\alpha )}(𝐫)`$ $`(\alpha 3)`$ are surface integrals over the surface of the particle $`i`$. For spherical particles, it can be omitted. In general, it is given by

$$C_i(\widehat{𝐫})=\sqrt{det𝐠_i(\widehat{𝐫})/det𝐠_{i,\text{sph}}(\widehat{𝐫})},$$
(14)
where $`𝐠_i(\widehat{𝐫})`$ is the local metric tensor of the particle $`i`$, and $`𝐠_{i,\text{sph}}(\widehat{𝐫})`$ that of a reference sphere of radius $`R_i(\widehat{𝐫})`$. (In polar coordinates, $`\sqrt{det𝐠_{i,\text{sph}}(\widehat{𝐫})}=R_i^2(\widehat{𝐫})\mathrm{sin}(\theta )`$). Finally, the factor $`\lambda `$ in Eqn. (11) ensures that the dimensional crossover for hard spheres is reproduced correctly and that the hard sphere system exhibits a solid-fluid transition in three dimensions .

$$\lambda =\frac{\left(1\xi (𝐫)^2\right)^3}{13\xi (𝐫)^2}\text{with}\xi =\left|\frac{𝐧_2(𝐫)}{n_2(𝐫)}\right|.$$
(15)
Alternatively, Eqn. (11) can be replaced by a more modern version due to Tarazona , which depends on an additional tensorial weighted density and does not contain the ad hoc factor $`\lambda `$.

## 3 Construction of the functional

As mentioned before, the Rosenfeld functional describes hard sphere fluids very successfully. On principle, it could also be applied to molecular fluids. However, it has the serious drawback that it does not reduce to the Onsager functional in the low density limit. Even worse, a closer inspection of Eqns. (8)–(15) reveals that it has no contribution at all which would favor parallel alignment of particles in a homogeneous fluid. Hence it cannot produce stable homogeneous anisotropic fluids.

The reason for this failure can be understood by looking at the relation between the Mayer-function and the weight functions in more detail. For hard spheres, the Mayer-function for a pair of particles at position $`𝐫_i`$ and $`𝐫_j`$ can be decomposed exactly as

$$\text{M}_{ij}=w_i^{(0)}w_j^{(3)}+w_i^{(3)}w_j^{(0)}+w_i^{(1)}w_j^{(2)}+w_i^{(2)}w_j^{(1)}𝐰_i^{(1)}𝐰_j^{(2)}𝐰_i^{(2)}𝐰_j^{(1)},$$
(16)
where $``$ denotes the convolution product:

$$w_i^{(\alpha )}w_j^{(\beta )}=𝑑𝐫w^{(\alpha )}(𝐫𝐫_i)w^{(\beta )}(𝐫𝐫_j).$$
(17)
This decomposition together with Eqns. (8)–(12) ensures that the functional reproduces the correct virial expansion at least up to second order. Unfortunately, a decomposition of the type above is no longer valid for anisotropic particles. Therefore, the Onsager limit is not recovered. Chamoux and Perera  have proposed ways to cure the virial expansion on the level of the direct correlation function, for the case of isotropic molecular fluids . A systematic way of dealing with the problem on the level of the functional itself would be to add a correction term on the right hand side of Eqn. (16), to deconvolute it (if possible), and then rederive a Rosenfeld type functional on the basis of this new decomposition. This can be done for rod-sphere interactions . An extension to the general case is currently under way . Unfortunately, it turns out that infinitely many additional terms are required in the decomposition in order to make Eqn. (16) exact, and the numerical treatment becomes difficult unless one resorts to approximations.

Here we propose a simpler, numerically perhaps more tractable Ansatz – a straightforward modification of Eqns. (9), (10) and (11). Instead of keeping $`n_1`$ and $`𝐧_1`$ as separate quantities, we suggest to replace the products $`n_1(𝐫)n_2(𝐫)`$ and $`𝐧_1(𝐫)𝐧_2(𝐫)`$ by new joint quantities $`n_{12}(𝐫)`$ and $`\stackrel{~}{n}_{12}(𝐫)`$ such that the functional reproduces correctly the Onsager limit. This is achieved as follows: we define the functions

$`g(\varpi _i,\varpi _j,𝐫_i𝐫_j)g_{ij}`$ $`:=`$ $`M_{ij}w_i^{(0)}w_j^{(3)}w_i^{(3)}w_j^{(0)}`$ (18)
$`h(\varpi _i,\varpi _j,𝐫_i𝐫_j)h_{ij}`$ $`:=`$ $`w_i^{(1)}w_j^{(2)}+w_i^{(2)}w_j^{(1)}𝐰_i^{(1)}𝐰_j^{(2)}𝐰_i^{(2)}𝐰_j^{(1)},`$ (19)
and

$$H^{\text{eff}}(\varpi _i,\varpi _j,𝐫_i𝐫_j)=g_{ij}/h_{ij}.$$
(20)
Then we replace $`n_1(𝐫)n_2(𝐫)`$ and $`𝐧_1(𝐫)𝐧_2(𝐫)`$ by

$`n_{12}(𝐫)`$ $`:=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑\varpi _i𝑑\varpi _j𝑑𝐫_i𝑑𝐫_j\rho (𝐫_i,\varpi _i)\rho (𝐫_j,\varpi _j)H^{\text{eff}}(\varpi _i,\varpi _j,𝐫_i𝐫_j)}`$
$`\times \left(w_i^{(1)}(𝐫_i𝐫)w_j^{(2)}(𝐫_j𝐫)+w_i^{(2)}(𝐫_i𝐫)w_j^{(1)}(𝐫_j𝐫)\right)`$
$`\stackrel{~}{n}_{12}(𝐫)`$ $`:=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑\varpi _i𝑑\varpi _j𝑑𝐫_i𝑑𝐫_j\rho (𝐫_i,\varpi _i)\rho (𝐫_j,\varpi _j)H^{\text{eff}}(\varpi _i,\varpi _j,𝐫_i𝐫_j)}`$
$`\times \left(𝐰_i^{(1)}(𝐫_i𝐫)𝐰_j^{(2)}(𝐫_j𝐫)+𝐰_i^{(2)}(𝐫_i𝐫)𝐰_j^{(1)}(𝐫_j𝐫)\right).`$
The effective function $`H^{\text{eff}}`$ is constructed such that the functional reduces to the Onsager limit (3) in the low density limit. Thus the mean curvature $`H(\widehat{𝐫})`$ in the weight functions $`w^{(1)}`$ and $`\text{w}^{(1)}`$ (Eqn. (13)) is effectively replaced by a function $`HH^{\text{eff}}`$ which depends on the orientations of a pair of particles, and their distance vector. In the hard sphere limit, $`H^{\text{eff}}`$ is constant, $`H^{\text{eff}}1`$, and we recover the original Rosenfeld functional. Note that $`n_{12}`$ and $`\stackrel{~}{n}_{12}`$ contain information on pairs of interacting particles. We thus give up the idea of formulating a functional which depends only on single, orientation independent weighted densities . An approach in the same spirit has been introduced by Schmidt for mixtures of spheres and needles .

The functional can be simplified by making the approximation that $`H^{\text{eff}}`$ only depends on the orientations $`\varpi _i`$ and $`\varpi _j`$. In that case we solely require that the second virial coefficient is exact in the *homogeneous* fluid for every distribution of orientations. In the homogeneous fluid the contributions to the free energy density having vectorial character vanish. Carrying out the integrals on the spatial variables, we obtain

$$\frac{\beta _R^{ex}}{N}=\mathrm{ln}\left(1\eta \right)+\frac{\rho \overline{R}S𝑑\varpi 𝑑\varpi ^{}f(\varpi )f(\varpi ^{})H^{\text{eff}}(\varpi ,\varpi ^{})}{1\eta }+\frac{1}{24\pi }\frac{\rho ^2S^3}{(1\eta )^2},$$
(23)
where $`\overline{R}`$, $`S`$ are the mean radius and the surface of the body, respectively, and $`\eta `$ is the packing fraction. The mean radius is defined as the integral of the mean curvature over the surface of the particle, $`\overline{R}=_S𝑑AH/4\pi `$. Eqn. (23) must reduce to the Onsager functional (4) in the low density limit $`\rho 0`$. With $`\eta =v_0\rho `$ this leads to the equation

$$H^{\text{eff}}(\varpi ,\varpi ^{})=\frac{\frac{1}{2}\nu (\varpi ,\varpi ^{})v_0}{\overline{R}S},$$
(24)
which replaces Eqn. (20) in this approximation.

In the case of homogeneous fluids, both (20) or (24) give the same bulk free energy as a function the orientation distribution $`f(\varpi )`$:

$$\frac{\beta _R^{ex}}{N}=\mathrm{ln}\left(1\eta \right)\frac{\eta }{1\eta }+\frac{\rho }{2}\frac{\varpi 𝑑\varpi ^{}f(\varpi )f(\varpi ^{})\nu (\varpi ,\varpi ^{})}{1\eta }+\frac{1}{24\pi }\frac{\rho ^2S^3}{(1\eta )^2}$$
(25)
Note that for isotropic fluids, $`f(\varpi )`$ is constant, and the free energy (25) is identical to that obtained with the original Rosenfeld functional.

## 4 Application to spherocylinders and ellipsoids

We have tested our approach by calculating the location of the nematic-isotropic phase transition for hard spherocylinders and uniaxial prolate ellipsoids. The results were compared with the corresponding phase diagrams obtained from the Onsager theory, the Parsons-Lee theory, and simulation data .

A spherocylinder consists of a cylinder of length $`L`$ and diameter $`D`$ capped by a hemisphere of diameter $`D`$ at both ends (see Figure 2).

The excluded volume of two spherocylinders which have the angle $`\theta `$ with respect to each other is

$$\nu (\theta )=\mathrm{\hspace{0.25em}2}DL^2\mathrm{sin}(\theta )+2\pi D^2L+4\pi D^3/3.$$
(26)
In the case of ellipsoids, the exact calculation of the excluded volume is quite involved. One possibility is to follow Ref. , another to exploit the Perram and Wertheim routine for the ellipsoids contact function . Here, we have adopted a scheme outlined in Ref. , because it is sensibly faster. We consider a pair of equal uniaxial prolate ellipsoids with semiaxes $`a`$, $`a`$ and $`b`$; one ellipsoid is mapped onto a sphere of radius $`a`$, and the other one is mapped onto a particular biaxial ellipsoid which has always one semiaxis equal to $`a`$. The remaining two semiaxes depend on the relative orientation of the two ellipsoids and can easily be evaluated. The excluded volume between the sphere and the biaxial ellipsoids can be calculated using an expression by Kihara for the orientationally averaged excluded volume between two convex hard particles . The covolume of the original pair of uniaxial prolate ellipsoids is then recovered by inverting the mapping. In the past, many theoretical studies  have considered systems of “hard Gaussian overlap” particles, i. e., ellipsoids where the contact distance is approximated by an expression due to Berne and Pechukas . However, it has been noted some years ago  and shown recently  that the comparison of these theoretical results with true ellipsoids behaviour is not appropriate.

Both in the cases of hard spherocylinders and hard ellipsoids, the properties of the fluids are fully determined by the density and the shape anisotropy parameters, $`x=L/D`$ and $`k=b/a`$, respectively.

In the homogeneous fluid, the orientation dependent part of the excess free energy has the same form in the Onsager model, the Parsons-Lee model, and our model:

$$\frac{\beta ^{ex}[f(\varpi )]}{N}=\frac{\lambda }{2}𝑑\varpi 𝑑\varpi ^{}f(\varpi )f(\varpi ^{})\nu (\varpi ,\varpi ^{}),$$
(27)
where $`\lambda =\rho `$ in the Onsager model, $`\lambda =\frac{\rho }{8}\frac{(43\eta )}{(1\eta )^2}`$ in the Parsons-Lee model, and $`\lambda =\frac{\rho }{2(1\eta )}`$ in our model. Once a model free energy functional is chosen, the next step is to determine the thermodynamically stable phase for a range of densities. We thus need to calculate the functions $`f(\varpi )`$ that extremize the free energy under the constraint $`𝑑\varpi f(\varpi )=1`$. For all three models under consideration, this amounts to finding the solution of the Onsager integral equation

$$\mathrm{ln}(\xi f(\varpi ))=\lambda 𝑑\varpi ^{}f(\varpi ^{})\nu (\varpi ,\varpi ^{}),$$
(28)
where the constant $`\xi `$ is determined from the normalisation condition. The equation always has an isotropic solution $`f(\varpi )=constant`$, and may have an anisotropic solution in addition at certain values of $`\lambda `$.

We have solved the Onsager integral equation for values of $`\lambda `$ in the range of $`\lambda ^{}=2DL^2\lambda [8:12]`$ for spherocylinders, and $`\lambda ^{}=2ab^2\lambda [8:25]`$ for ellipsoids, in steps of 0.01. To this end, we have applied an iterative numerical method, which is simple and reliably convergent . The integrals were calculated by the Gauss-Legendre quadrature, using 50 points per integral and checking that with 100 points we get the same results. The initial guess for the highest value of $`\lambda ^{}`$, $`\lambda ^{}=12`$ or $`\lambda ^{}=25`$, was

$$f(\theta ,\varphi )=\frac{1}{2\pi }\frac{\mathrm{exp}\left(P_2(\mathrm{cos}\theta )\right)}{𝑑u\mathrm{exp}\left(P_2(u)\right)},$$

where $`P_2(u)`$ is the second Legendre polynomial. At lower $`\lambda ^{}`$, the iteration was started with the anisotropic solution for the previous, next higher, value of $`\lambda ^{}`$. The convergence criterion was

$$\left[\underset{i=1}{\overset{n_p}{}}\left[f_{m+1}(\theta _i)f_m(\theta _i)\right]^2\right]^{\frac{1}{2}}<10^6$$
(29)
where *$`n_p`$* refers to the number of points used to perform integrals and *m* to the *m*th iteration. The dependence on $`\varphi `$ has been omitted because the distribution function does not actually depend on it. The resulting distribution functions for $`\lambda ^{}=11.01`$, 10.01, 9.01, 8.01 for spherocylinders are shown in Figure 3. Our orientational distribution functions compare well with those in reference . Due to the particular shape of the covolume between a pair of spherocylinders, the solutions depend on $`\lambda ^{}`$, but not on $`x`$. We found anisotropic solutions for values $`\lambda ^{}8.88`$, in agreement with Lasher .

In the case of ellipsoids, the solution for fixed $`\lambda ^{}`$ depend on the shape anisotropy parameter $`k`$. In Figure 4, solutions at different value of $`k`$ are shown for $`\lambda ^{}=15.01`$. The lowest value of $`\lambda ^{}`$ which still yields an anisotropic solution is a decreasing function of $`k`$.

Next we must identify the stable phases and the coexistence line. At coexistence, both the pressure and the chemical potential are equal:

$$P_N=P_i;\mu _N=\mu _i.$$
(30)
We have solved these two equations with Newton’s method, i.e., we found the density $`\rho _I^P`$ at which an isotropic phase has the pressure $`P_N`$ and the density $`\rho _I^\mu `$ at which an isotropic phase has the chemical potential equal to $`\mu _N`$. At the coexistence line, both densities must be equal. Hence we calculated the nematic $`\rho _N`$ at which $`|\rho _I^P\rho _I^\mu |`$ was minimal. Usually the value of the minimum was in the range $`|\rho _I^P\rho _I^\mu |[10^9;10^5]`$.

The results are shown as a function of the anisotropy in Figure 5 for spherocylinders and in Figure 6 for ellipsoids. They are compared with simulation data of Bolhuis and Frenkel , Frenkel and Mulder , and Camp *et al.* .

The figures demonstrate that our functional performs much better than the Onsager functional, and only slightly worse than the Parsons-Lee functional. As has already been demonstrated elsewhere for spherocylinders and ellipsoids , the predictions of the latter are almost exact. The slightly inferior performance of our functional is probably related to the fact that the Parsons-Lee approach is based on the Carnahan-Starling equation of state, whereas the Rosenfeld functional yields the Percus-Yevick equation of state for hard spheres, which is slightly inferior. On principle, one can modify the Rosenfeld functional such that it reproduces the Carnahan-Starling equation of state for fluids, along the lines of an approach suggested by Tarazona . Unfortunately, this is done at the expense of a less accurate description of crystals .

Our approach could possibly be improved by making the coefficient of $`\mathrm{\Phi }_3`$ in Eqn. (11) dependent on the orientation distribution $`f(\varpi )`$: formally, the Rosenfeld functional has the form of a third order $`Y`$-expansion . Mulder and Frenkel  and Tjipto-Margo and Evans  have generalized the latter to convex anisotropic bodies and applied it to hard ellipsoids, determining the third virial coefficient numerically. Their nematic-isotropic coexistence densities were lower than those observed in the simulations. In contrast, our functional tends to overestimate the coexistence densities. This suggests that introducing orientation dependent coefficients will shift the transition lines in the correct direction.

However, these modifications complicate the functional. We have shown that already our simple version describes uniform fluids reasonably well. Compared to the Parsons-Lee theory, it has the advantage of being based on a reference density functional which describes the local structure of hard spheres very accurately. It interpolates on the level of the direct correlation function, i. e. the local bulk structure, between the Percus-Yevick solution for hard spheres, which is very good, and the Onsager solution for infinitely elongated particles, which is exact. Therefore, we believe that it will be suited to describe inhomogeneous anisotropic fluids.

## 5 Summary and Outlook

We have introduced a new density functional for liquid crystals, which interpolates between the successful Rosenfeld functional for hard sphere fluids  and the Onsager functional for liquid crystals . As a first test, we have calculated the nematic-isotropic phase diagram for hard spherocylinders and hard ellipsoids and shown that the new functional produces reasonable results. In the next step, it shall be applied to calculate the local liquid structure in homogeneous fluids. For example, direct correlation functions in nematic fluids of ellipsoids have been determined recently from computer simulations . They can also be calculated from our density functional. This will be a second, much more sensitive test of the theory.

If the new functional is successful, it will allow to predict the structure and the properties of nematic liquid crystals in the vicinity of inhomogeneities. It will thus contribute to an improved microscopic understanding of surface phenomena such as surface anchoring on rough and structured surfaces, or of defect structures and defect interactions. It also provides a new approach to smectic ordering, and might lead to useful insight into the nature of the nematic-smectic transition .

## Acknowledgments

We thank Guido Germano for fruitful interactions. We are grateful to Yasha Rosenfeld for useful comments on the manuscript and for pointing out reference  to us. This work was funded by the Scuola Normale Superiore di Pisa and by the German Science Foundation.


format: text

text:
# Abundances of Molecular Species in Barnard 68

## 1 Introduction

Stars form out of cores of dense molecular gas. Such cores can be embedded deeply within molecular clouds, but numerous examples of isolated cores, e.g., Bok globules (Bok & Reilly 1947; Clemens & Barvainis 1988), have been also identified as dark patches of high visual extinction against rich stellar backgrounds. As isolated cores, Bok globules can be relatively free from the influences of nearby star-forming events, and so can provide very simple and pristine examples of the physical and chemical structures of cores prior to the formation of low-mass stars. In this paper, we examine the abundances of several species of molecules or molecular ions in the globule Barnard 68 (B68; Barnard 1919).

B68 (also known as CB82 or LDN57) is an Ophiuchus dark cloud, located at $``$125 pc (de Geus, de Zeeuw, & Lub 1989)111A distance range of 60-205 pc for the Ophiuchus dark clouds was formally estimated by de Geus et al. with 125 pc $`\pm `$ 25 pc defined as their center., within the Loop I superbubble. Figure 1 shows an $`R`$-band image of B68 obtained from the CDS/Aladin sky atlas, revealing its compact, roundish morphology. B68 has no evidence of star formation within its interior, e.g., from detections of outflow wings in lines of 12CO (Avery et al. 1987) or from an IRAS point source (Parker 1988)222The 24′′ $`\times `$ 5′′ (P.A. of 94) error ellipse of the IRAS point source 17194-2351, associated with B68 by Clemens & Barvainis, is centered $``$3 beyond the opaque edge of B68 seen in Figure 1.. Furthermore, continuum emission at $`\lambda `$ = 1.3 mm toward B68 has not been detected (Reipurth, Nyman, & Chini 1996; Launhardt & Henning 1997).

B68 is unique in that it is the only starless Bok globule to have its density structure well defined. Using the extinction of the background Galactic bulge K-giant population, Alves, Lada, & Lada (2001; ALL01) found an azimuthally-averaged radial column density profile for B68 that was well matched with that of an isothermal, self-gravitating Bonnor-Ebert sphere, specifically one with a near-critical center-to-edge density contrast of 16.5 (Bonnor 1956; Ebert 1955). The robust definition of density within B68 removes a major uncertainty towards determining its molecular abundances. Recently, Bergin et al. (2002) used single-dish maps of B68 in conjunction with models specifically incorporating the ALL01 Bonnor-Ebert density structure to ascertain radial variations in C18O and N2H+ abundance, with the lowest values occurring at the highest extinctions within B68. Hotzel et al. (2002) also found evidence for depletion of CO isotopes with extinction from their own single-dish maps and a re-analysis of the ALL01 data. Each group speculated these molecules sublimate from the gas phase onto dust grains at the high densities and low temperatures in the B68 interior. A similar interpretation was made by Tafalla et al. (2002) to explain the radial depletion of CS and C18O in the interiors of 5 other starless cores, although N2H+ was not depleted and NH3 was enhanced at these locations.

We explore further the abundances of several molecules or molecular ions within B68, incorporating the ALL01 Bonnor-Ebert density structure into our own models. New observations of B68 made with the BIMA millimeter array are used to obtain upper limits on abundances of 13CO and HCO+ near the extinction peak found by ALL01, i.e., at R.A./decl. (2000) = 17h22m38.6s, -234946.0′′. Upper limits to C18O and N2H+ abundances are also derived at this position, consistent with values derived by Bergin et al. In addition, single-dish detections of lines reported in the literature are used to derive abundance values for C18O, CS, N2H+, NH3, H2CO, and C3H2 toward positions in B68 offset from the extinction peak. We describe the BIMA observations and single-dish data in §2, and our method of determining abundances with a Monte Carlo radiative transfer code in §3. The derived results are discussed in §4. A concluding summary is found in §5.

## 2 Observations and Results

### 2.1 BIMA Observations

Figure 1 shows the position of the one pointing made toward B68 with the BIMA millimeter interferometer at Hat Creek, CA, i.e., at R.A./decl. (2000) = 17h22m38.6s, -234942.0′′, approximately 4′′ north of the extinction peak. Both 13CO 1–0 and C18O 1–0 were observed simultaneously during tracks on 2001 May 09 and May 20, and N2H+ 1–0, HCO+ 1–0, and HOC+ 1–0 were observed simultaneously on 2001 April 24 and May 14. Figure 1 also shows the minimum and maximum FWHM sizes of the BIMA primary beams at the frequencies of 13CO 1–0 and HCO+ 1–0. Each line was observed in its own correlator window of 6.25 MHz width, with each window set to have 512 channels of 12.2 kHz width (i.e., $``$0.03 km s-1 at 110 GHz.) The correlators were tuned so the central channel of each window corresponded to a $`V_{LSR}`$ of $`+`$3.4 km s-1, the velocity centroid of C18O 2–1 found by Wang et al. (1995).

Each track was only 4-5 hours in duration given the low maximum elevation of B68 from Hat Creek (i.e., $``$25). Visibility phases were calibrated by observing 1733-130 for 5 minutes approximately every 30 minutes over 800 MHz bandwidth in the lower or upper sideband. Visibility amplitudes were calibrated by similarly observing Uranus for 8 minutes at the end of each track. All data were reduced using standard routines within the MIRIAD software package (Sault, Teuben, & Wright 1995). All channels were cleaned to 2 $`\sigma `$ rms intensity levels.

Line emission was not detected at any of the frequencies observed, at any positions toward B68 within the $``$2 FWHMs of the BIMA primary beams. Binning together several channels or tapering the data with variously-sized Gaussians did not change this result. Table 1 summarizes these data, listing the lines, the source observatory, the synthesized beam FWHMs, and the angular offset of the BIMA pointing from the extinction peak. The 1 $`\sigma `$ rms sensitivities attained per channel ranged from 1.0 K for HOC+ 1–0 to 1.5 K for 13CO 1–0, and were measured by sampling all channels within a 40′′ $`\times `$ 40′′ box centered at the pointing center. Upper limits (at 3 $`\sigma `$) to the integrated intensities of these lines ranged from 0.21 Jy beam-1 km s-1 to 0.34 Jy beam-1 km s-1 (or 0.25 K km s-1 to 0.44 K km s-1) assuming a line width equal to 0.4 km s-1, the typical FWHM of optically-thin lines in B68 found by Wang et al.

### 2.2 Single-dish Observations

We surveyed the literature to find previous detections of molecular line emission toward B68 from which abundances could be measured. Fortunately, B68 was included in numerous surveys of dark globules for molecular line emission. Early detections of lines include those of 12CO 1–0 and 13CO 1–0 by Martin & Barrett (1978) and Leung, Kutner, & Mead (1982), 12CO 2–1 by Avery et al. (1987) and Clemens, Yun, & Heyer (1991), and 12CO 3–2 by Avery et al. Furthermore, NH3 (1,1) and (2,2) were detected by Martin & Barrett and Bourke et al. (1995). Although these data certainly associate line emission with B68, we exclude them from analysis because of their relatively low angular resolution, i.e., $`>`$60′′ FWHM.

More recent detections of line emission from various species toward B68 have been made with resolutions $`<`$60′′ FWHM, including those of C18O 2–1 and H2CO 312–211 by Wang et al., CS 2–1 by Launhardt et al. (1998), N2H+ 1–0 and C3H2 202–101 by Benson, Caselli, & Myers (1998), and NH3 (1,1) by Lemme et al. (1996). Other recent detections include those of C18O 1–0 and N2H+ 1–0 by Bergin et al., and 13CO 1–0, C18O 1–0 and 2–1 by Hotzel et al., but we exclude those data from our analysis because abundances have been specifically determined from them by the respective authors. An interesting non-detection in the literature is that of SO 10–01 by Codella & Muders (1997), but their map was centered on the nearby IRAS point source and did not overlap the region of extinction modeled by ALL01. The methods of observation and reduction relevant to these data are described in these references, and are not reproduced here.

In the higher-resolution single-dish studies, the positions where line characteristics were specifically reported vary widely. Figure 1 shows the positions and resolutions of observations from these studies against the $`R`$-band image of B68. Table 1 also lists the lines observed, the source observatories, the resolutions attained, and the angular offsets of the single-dish data from the extinction peak of the reported line characteristics. The positions listed in Table 1 for C18O 2–1, CS 2–1, and N2H+ 1–0 are those of peak line intensity from their respective maps, where line characteristics were specifically reported. (H2CO 312–211 and C3H2 202–101 were observed at only these positions by the respective authors.) Interestingly, these positions are not coincident with the extinction peak. The C18O and CS positions are approximately coincident with an arc of maximal C18O 1–0 integrated intensity centered at the extinction peak noted by Bergin et al. (see their Figure 1a), and the N2H+ position is approximately coincident with the “arc” of maximal N2H+ 1–0 integrated intensity about the extinction peak also noted by Bergin et al. (see their Figure 1b.)333An N2H+ 1–0 map toward B68 recently made by Caselli et al. (2002), confirms this position is that of peak intensity, although their data are not Nyquist sampled and are of lower resolution than those of Benson et al. and Bergin et al.

## 3 Monte-Carlo Radiative Transfer Models

Molecular line data can place strong limits on molecular abundances within B68 because its density structure has been so well constrained by ALL01. For example, abundances can be estimated by calculating with a Monte Carlo code the radiative transfer, molecular excitation, and line emission through a model of a dense core with a specific density profile, and varying abundances within the model until the output matches the observed data. For this purpose, we used the one-dimensional Monte Carlo code of Hogerheijde & van der Tak (2000) which can solve molecular excitation coupled with line and continuum radiative transfer in core models. This code is especially useful for application to B68. First, model fluxes can be found at positions arbitrarily offset from the central line-of-sight, to simulate observations of B68 at positions offset from the extinction peak. Second, the model fluxes can be sampled with the spatial frequency coverage of actual interferometer data to estimate how much emission may have been resolved out in those cases. Recently, van Zadelhoff et al. (2002) compared 8 molecular excitation/radiative transfer codes including that of Hogerheijde & van der Tak, and found that all codes give results to within 10% of one another, even at high opacity.

For models of B68, we assumed the thermal and density structure given by ALL01, i.e, an isothermal sphere with an outer radius of 12 500 AU and a “Bonnor-Ebert parameter” $`\xi _{\mathrm{max}}`$ of 6.9, corresponding to a center-to-edge density ratio of 16.5. Following ALL01, the temperature assumed in our models was 16 K, as derived by Bourke et al. from NH3 observations of B68. (An alternative temperature profile is considered in §4.1 below.) The $`\xi _{\mathrm{max}}`$ parameter identifies the particular solution from the family that solves the second-order differential equation characterizing a Bonnor-Ebert sphere. We determined the density in 30 linearly-spaced concentric shells by solving that differential equation with a fourth-order Runge-Kutta method (see Press et al. 1992.) The central density in the ALL01 model for B68 is $``$2.5 $`\times `$ 105 cm-3.

Molecular abundances were assumed to be constant with radius, although abundance gradients within B68 are quite possible for the molecular species considered here. For example, Bergin et al. found edge-to-center contrasts in the abundances of C18O and N2H+ of 100 and 2 respectively using their own maps of B68. With only single pointings for the single-dish data and non-detections for the interferometer data, we can derive abundance values or upper limits only at the positions listed in Table 1, averaged over the respective lines-of-sight and beam widths, and cannot probe for abundance gradients. Gradients should not produce dramatic differences in the abundances we derive, however, given the small sizes of the beam widths relative to the spatial extent of B68. For example, we note in §4.2 little difference between our abundances of C18O and N2H+ and those derived by Bergin et al. at positions observed in common, despite our different abundance profiles.

The velocity field of the gas was assumed to have a turbulent line width of 0.4 km s-1 (see Wang et al.; since the lines are optically thin at the abundances found, the integrated intensity is independent of the adopted turbulent line width.) We assume that each shell is stationary, although Bergin et al. suggest velocity variations within B68 may be significant. For reasons similar to those cited above for abundance gradients, we did not consider radial variations of line width.

In each shell, level populations were determined assuming standard collision rates for each molecule. The expected sky brightness distributions were then calculated using a $`256\times 256\times 100`$ cube with $`1^{\prime \prime }\times 1^{\prime \prime }\times 0.04`$ km s-1 elements. To compare models with interferometer data, “visibility” datasets were produced from the cubes using the MIRIAD task uvmodel and the antenna baselines from the original data. These latter datasets were inverted, cleaned and restored in the same manner as the actual data, and velocity-integrated line intensities calculated. By properly accounting for spatial filtering in this manner, we estimated BIMA would have recovered only $``$20% of the flux emitted from an ALL01 Bonnor-Ebert sphere in B68. To compare models with single-dish data, the original cubes were convolved with Gaussians with widths appropriately representing the resolution of the respective observations. The abundances were varied until the observational data were well matched.

## 4 Abundances

Table 2 lists the fractional abundances of the sampled molecular species at various positions toward B68, constrained with the data described in §2, and estimated by using the method described in §3 (assuming models with the Bonnor-Ebert density profile of ALL01.) The BIMA data yield the first upper limits to the abundances of HCO+ toward the extinction peak of B68, as well as upper limits to the abundances of 13CO, C18O and N2H+ at the same location. Since collision rates are not available for HOC+, we do not attempt to place limits on its abundance. In addition, single-dish data of B68 from the literature yield abundance values of H2CO, CS, NH3, and C3H2 toward B68 for the first time, but at various positions offset from the extinction peak. Finally, other single-dish data of B68 from the literature yield new abundance values of C18O and N2H+ at some of these latter locations.

### 4.1 Comparison with Other Clouds

How do the abundances derived for B68 here compare to those derived for other clouds? “Standard” cloud molecular abundances with which to base comparisons are difficult to define given wide variations in chemical environment or evolutionary epoch in the ISM. To provide some comparison, Table 2 also lists the fractional abundance estimates of the same species in clouds similar in character to B68, obtained directly from the literature but made without the detailed foreknowledge of cloud density structure now available for B68. The first set, for “translucent clouds,” are those listed by Turner (2000) for small, round Clemens & Barvainis clouds with edge-to-center visual extinctions of 2.0. These clouds were modeled with constant abundance, and the listed values are averages between the results of a hydrostatic equilibrium polytropic model and an $`n`$($`r`$) $``$ $`r^0`$ model444According to Turner (2002, private communication), these two models produce results that differ at most by a factor of 1.35 for species that are highly dependent on density.. The second set, denoted for “cold dense clouds,” are those compiled by Ohishi, Irvine, & Kaifu (1992) for a single position in TMC-1 or L134N, and were made using column density estimates from many authors, assuming N(H2) = 1022 cm-2 (see Turner, or Pratap et al. 1997, or Dickens et al. 2000 for alternative estimates.) The abundances from Turner and Ohishi et al. in Table 2 should not be regarded as universal. For example, Turner notes abundances can vary by an order of magnitude both between different translucent clouds and within larger, dense clouds (presumably determined using similar assumptions.)

Table 2 shows that every value of molecular abundance we derive for B68 is less than the lowest value derived for other clouds. From the BIMA data, we find the upper limits to the 13CO, HCO+, and N2H+ abundances in B68 are lower than values derived for other clouds by factors of 7, 14, and 4 respectively. (The upper limit to the C18O abundance is similar to the lowest value found for other clouds.) From the various single-dish data, we find the C18O, CS, and H2CO abundances differ least, by roughly an order of magnitude or less, with the B68 abundances lower than the lowest values derived for other clouds by factors of only 5, 3, and 16 respectively at the locations observed. However, the N2H+, NH3, and C3H2 abundances differ the most, by over an order of magnitude, with those of B68 lower by factors of 25, 29, and 170 respectively at the locations observed.

The assumption in our models of an isothermal Bonnor-Ebert sphere, at the same 16 K temperature derived by Bourke et al. and used by ALL01, has likely minimized the resulting abundance values we obtain. Other temperature profiles, such as those suggested for Bonnor-Ebert spheres by Zucconi, Walmsley, and Galli (2001; see also Evans et al. 2001) which have cooler temperatures at small radii, may be more appropriate. For example, Bergin et al. assumed for their models of B68 the Zucconi et al. radial temperature profile with a global reduction of 2 K, in addition to the ALL01 density profile. With this same temperature profile in our models, the values in Table 2 increase by factors of only 2–3. Therefore, using this radial temperature profile would bring the abundances of C18O, CS, and H2CO at the positions observed in B68 even more in line with the abundances found for other cloud types. However, the abundances of N2H+, NH3, and C3H2 at other positions would still remain lower than those of other clouds by an order of magnitude or more.

Only N2H+, NH3, and C3H2 are arguably depleted in B68, given the values shown in Table 2 and assuming a typical order of magnitude abundance variation between and within similar clouds. However, other species considered here still may be depleted in B68. Note that the degree to which the abundance of a particular species in B68 is low compared to those of other clouds appears related to the position across B68 that the respective data were obtained. The C18O, CS, and N2H+ data from the literature were those of maximum line brightness from maps of B68, and H2CO and C3H2 data were obtained only subsequently by the respective authors at those positions. (NH3 was observed only at the Clemens & Barvainis position of B68.) Figure 1 and Table 2 together reveal that the less-discrepant abundances are found at positions relatively far from the extinction peak but the more-discrepant abundances are found at positions closer to the extinction peak (except notably NH3.) This pattern suggests C3H2 may be also depleted by some process related to extinction, e.g., the sublimation of gas-phase molecules onto grains, as suggested by Bergin et al. for C18O and N2H+ and Hotzel et al. for 13CO and C18O. This same idea may also explain how CS and H2CO appear relatively undepleted in the outer, less-extincted radii of B68. Moreover, the lack of bright emission in CS near the extinction peak hint it may be also depleted at high extinction in B68. These speculations can be confirmed only after analyzing fully-sampled, high-resolution line maps of B68.

### 4.2 Comparison with Previous Abundance Determinations

How do the abundances we derive compare with those found previously for B68? Bergin et al. found the C18O abundance in B68 rises from very low values at $`A_V`$ $`<<`$ 1 to a peak of 1 $`\times `$ 10-7 at $`A_V`$ = 2, and decreases to 1 $`\times `$ 10-9 at $`A_V`$ $`>`$ 20, a contrast of $``$100. In addition, Bergin et al. found the N2H+ abundance rises from very low values at $`A_V`$ $`<<`$ 1 to a peak of 6 $`\times `$ 10-11 at $`A_V`$ = 3, and decreases to 3 $`\times `$ 10-11 at $`A_V`$ $`>`$ 20, a contrast of $``$2. Any evidence of N2H+ depletion is remarkable, given its oft-described utility as a non-depleting probe of dense core interiors (e.g., see Tafalla et al. 2002.)

Our interferometer data of C18O and N2H+ do not provide much additional support for the abundance model for B68 of Bergin et al. Our derived upper limits for C18O and N2H+ from the BIMA data are not particularly low, but they remain consistent with the still-lower values found by Bergin et al. throughout the core. Also, the non-detection of compact line emission with BIMA suggests the low abundances found by Bergin et al. are not due to the single-dish beam dilution of small clumps of relatively abundant material. A spatially-smoother model of gas density is indeed more appropriate for B68, such as the Bonnor-Ebert sphere suggested by the dust extinction maps.

The single-dish C18O and N2H+ data from the literature only pertain to one line-of-sight per transition toward B68, but provide strong support for the abundances derived by Bergin et al., if the same modified Zucconi et al. temperature profile is assumed. With this assumption, the C18O literature data, from a position 59$`\stackrel{}{\mathrm{.}}`$4 offset from the extinction peak, yield an apparent abundance of 6–9 $`\times `$ 10-8, quite consistent with the 3–9 $`\times `$ 10-8 abundance range expected along this line-of-sight from the model of Bergin et al. Also, the N2H+ literature data, from a position 16$`\stackrel{}{\mathrm{.}}`$3 offset from the extinction peak, yield an apparent value of 4–6 $`\times `$ 10-11, slightly larger than but still consistent with the 3–4 $`\times `$ 10-11 abundance range expected along this line-of-sight from the model of Bergin et al. Furthermore, the observed integrated intensities of these literature data at each respective location and our BIMA observations can be reproduced to within only a few percent with our own Monte Carlo models, assuming the ALL01 density profile, the modified Zucconi et al. temperature profile and the Bergin et al. radial abundance profiles. (Note that here we relaxed the earlier assumption of constant abundance with radius.)

## 5 Summary and Conclusions

Using our own interferometer line data and single-dish line data culled from the literature of the Bok globule B68 and a Monte-Carlo radiative transfer code, we derive abundances for C18O, CS, N2H+, NH3, H2CO, and C3H2 and upper limits for abundances of 13CO and HCO+. Foreknowledge of the Bonnor-Ebert density configuration of B68, as found from extinction mapping by ALL01, has removed a major uncertainty that affected previous estimations of abundances in other clouds. We find molecular abundances in B68 are lower than those estimated previously for clouds similar to B68. N2H+, NH3, and C3H2 have abundances lower than the lowest values determined for other clouds by 1.4, 1.5, and 2.2 orders of magnitude respectively, but abundances of the other species are only lower by one order of magnitude or less. Depletion of CS is suggested by the large offset between the position of maximum CS line brightness and the extinction peak. Furthermore, abundances of C18O and N2H+ derived using data from the literature at locations where line characteristics were reported are consistent with those derived by Bergin et al.

As shown by Bergin et al., fully-sampled maps from single-dish telescopes can more effectively determine abundances within B68. Further maps of B68 can be now easily obtained using the new multi-beam focal-plane arrays on current single-dish telescopes. Furthermore, the upcoming generation of millimeter interferometers will have the sensitivity to provide additional data for combined maps of even higher resolution. With these data, more accurate molecular abundances will be determined for B68. For example, it will be interesting to determine if abundances of CS, HCO+, NH3, H2CO, and C3H2 also vary with extinction, like 13CO, C18O, or N2H+. Similar abundance determinations across other isolated dense cores may also be possible, provided their density structures can be similarly well-defined.

###### Acknowledgements.

We thank Barry Turner, our referee, for insightful comments that improved this paper. In addition, we thank Leo Blitz for allowing this project to be pursued with the BIMA millimeter array. JD’s research in Berkeley was supported by the Radio Astronomy Laboratory. MRH’s research in Berkeley was supported by the Miller Institute for Basic Research in Science. We also thank Charles J. Lada, João Alves, Tracy Huard, and Jon Swift for their help. This research has made use of the SIMBAD database and Aladin, both operated by CDS, Strasbourg, France.


text:
# Exact ground state number fluctuations of trapped ideal and interacting fermions

## I Introduction

Recent experimental observation of quantum degeneracy in a dilute gas of trapped fermionic atoms at low temperatures by DeMarco and Jin  inspired our previous work on fermions . We considered an ideal gas of spinless fermions in an isolated harmonic trap in one- and two-dimensions. In this work we restrict our study to a system of fermi gas in one-dimension. In addition to reporting some new results we obtained for a system of ideal fermi gas, we now introduce inter-particle interaction. Since our system is isolated (microcanonical ensemble), the results given in this paper for both ideal and interacting fermi gases are exact.

Consider first a system containing $`N`$ non-interacting particles in a one-dimensional harmonic confinement, either bosons or (spinless) fermions. The canonical partition function of this system is given by

$$Z_N(x)=x^{E_N(0)}\underset{j=1}{\overset{N}{}}\frac{1}{(1x^j)}=x^{E_N(0)}\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{\Omega }(n,N)x^n,$$
(1)
where $`x=e^\beta `$, $`\beta =1/k_BT`$, $`E_N(0)`$ is the ground state energy of the system. $`E_N(0)=N/2`$ for bosons and $`E_N(0)=N^2/2`$ for fermions. For simplicity we have set $`\mathrm{}\omega =1`$, where $`\omega `$ is the oscillator frequency. Note that the expansion coefficient $`\mathrm{\Omega }(n,N)`$ is the same for bosons and fermions, and so is the canonical partition function, apart from the prefactor involving the ground state energy $`E_N(0)`$. The ground state of the system is defined at $`T=0`$. For fermions the ground state consists of a set of $`N`$ occupied energy levels, while for bosons it is the single lowest energy level. Given an excitation energy $`n`$, there are many ways in which the particles can share this energy and excite to the higher states, so that the particle number in the ground state fluctuates. This ground state number fluctuation for ideal particles may be calculated exactly using combinatorics. This was done by us in Ref.  for ideal fermions. The desired quantity is the microcanonical multiplicity $`\omega (n,N_{ex},N)`$, defined as the number of ways to distribute exactly $`N_{ex}`$ out of $`N`$ particles from the ground state to the excited states, given $`n`$ quanta of energy in a system of $`N`$ particles. Once $`\omega (n,N_{ex},N)`$ is found, the corresponding probability $`P(n,N_{ex},N)`$ and the exact ground state number fluctuation can be calculated straightforwardly:

$`P(n,N_{ex},N)`$ $`=`$ $`{\displaystyle \frac{\omega (n,N_{ex},N)}{_{N_{ex}=1}^N\omega (n,N_{ex},N)}},`$ (2)
$`N_{ex}`$ $`=`$ $`{\displaystyle \underset{N_{ex}=1}{\overset{N}{}}}P(n,N_{ex},N)N_{ex},`$ (3)
$`N_{ex}^2`$ $`=`$ $`{\displaystyle \underset{N_{ex}=1}{\overset{N}{}}}P(n,N_{ex},N)N_{ex}^2,`$ (4)
$`(\delta N_0)^2`$ $`=`$ $`N_{ex}^2N_{ex}^2`$ (5)
$`=`$ $`N_0^2N_0^2.`$ (6)
Note that by definition $`P(n,N_{ex},N)`$ is normalized to unity. Eqs. (2)-(6) apply to any system of particles, either bosons or fermions. Note, however, that the microcanonical multiplicity $`\omega (n,N_{ex},N)`$ is different for different types of particles. For fermions, apart from the fact that the ground state consists of more than one energy levels, the distribution of $`N_{ex}`$ particles must comply with the Pauli Principle. For clarity, therefore, we shall attach a bracketed superscript $`\omega ^{\{B\}}(n,N_{ex},N)`$ for bosons, and $`\omega ^{\{F\}}(n,N_{ex},N)`$ for fermions in what follows.

Although the combinatorial method may be used to determine the microcanonical multiplicities $`\omega ^B(n,N_{ex},N)`$ and $`\omega ^F(n,N_{ex},N)`$ for both ideal bose and fermi gases exactly, it is very time-consuming computationally. In the case of an ideal boson gas in a one-dimensional harmonic trap, the problem is greatly simplified due to its connection to number partitioning theory. This is because the energy spectrum of a one-dimensional harmonic trap is equally spaced, and the $`N_{ex}`$ bosons are excited from a single lowest energy level. Given $`n`$ excitation quanta imparted to the system, the number of ways to excite exactly $`N_{ex}`$ bosons to the excited states, $`\omega ^B(n,N_{ex},N)`$, is equivalent to the number of ways of partitioning an integer $`n`$ into exactly $`N_{ex}`$ parts. It is related to the expansion coefficient $`\mathrm{\Omega }(n,N)`$ of the canonical partition function $`Z_N(x)`$ via :

$`\omega ^B(n,N_{ex},N)`$ $`=`$ $`\mathrm{\Omega }(nN_{ex},N_{ex}),nN_{ex}`$ (7)
$`=`$ $`0,\text{ }otherwise.`$ (8)
Note that Eq. (8) implies that $`\omega (n,N_{ex},N)=\omega (n,N_{ex},N^{}),N^{}N_{ex}`$, i.e. independent of system size. In the language of number theory, the canonical multiplicity $`\mathrm{\Omega }(n,N)`$ is the number of ways of partitioning $`n`$ into $`1,2,3,\mathrm{},N`$ parts. Clearly,

$$\mathrm{\Omega }(n,N)=\underset{N_{ex}=1}{\overset{N}{}}\omega ^B(n,N_{ex},N).$$
(9)
Since the canonical multiplicity $`\mathrm{\Omega }(n,N)`$ is the same for bosons and fermions (Eq. (1)), the following must be true:

$$\mathrm{\Omega }(n,N)=\underset{N_{ex}=1}{\overset{N}{}}\omega ^F(n,N_{ex},N).$$
(10)
So even though the bosonic $`\omega ^B(n,N_{ex},N)`$ is very different from the fermionic $`\omega ^F(n,N_{ex},N)`$, their sum $`\mathrm{\Omega }(n,N)`$ is the same. For bosons, expression (8) implies that knowing the canonical multiplicity $`\mathrm{\Omega }(n,N)`$ is sufficient to study the system at a microcanonical level. It therefore allows an easy way to calculate the exact number fluctuation, instead of the more cumbersome combinatorics formulae. This is not the case for fermions, however, since Eq. (8) does not apply. This is because the (quantum mechanical) ground state of fermions consists of $`N`$ occupied energy levels, unlike bosons whose ground state consists of only a single lowest energy level. To see why this multi-level ground state of fermionic system prevents the application of number theory, consider a system containing $`N=2`$ particles. For simplicity take $`n=2`$ and consider exciting one particle, i.e., $`N_{ex}=1`$. In the language of number theory, this means partitioning the integer $`n=2`$ into one part. There is only one way of doing this. For bosons there is also only one way of taking $`1`$ out of $`N=2`$ particles from the ground state, and put it in the second level above, i.e., $`\omega ^B(2,1,2)=1`$. For fermions, however, either the particle which lies in the fermi level $`E_F`$ may get excited to the second level above, or the particle which lies just below $`E_F`$ gets excited to the first level above. Thus $`\omega ^F(2,1,2)=2`$, which is not equal to one as in the case of bosons. Note that $`n=2`$ may also be partitioned into two parts, and $`\mathrm{\Omega }(2,2)=2`$ since $`2=2,1+1`$. For bosons this means exciting both particles, i.e. $`\omega ^B(2,2,2)=1`$, whereas for fermions $`\omega ^F(2,2,2)=0`$ since it takes at least $`4`$ quanta to excite $`2`$ fermions. In either case the microcanonical multiplicities sum up to $`2`$ (Eqs. (9)and (10)). So the multi-level ground state of fermionic system compensates for the restriction in the number of ways the particles may be distributed imposed by the Pauli Principle. Fortunately, the mapping of the system of $`N`$ fermions in a one-dimensional harmonic trap to number partitioning theory is not impossible. In section II, we show that by looking at the excited energy levels and the multi-level ground state separately, we are able to connect each of these two spaces to a similar bosonic problem, which in turn pertains to number partitioning. We successfully derive a formula similar Eq. (8) which applies for fermions. Apart from mathematical interest, the formula allows us to calculate the exact ground state number fluctuation more efficiently than the combinatorics method.

In the conventional treatment of the ideal bose gas, the number fluctuation is calculated using the grand canonical ensemble (GCE). The result is given by

$$(\delta N)(\delta N)^2^{1/2}=\left[\underset{k}{\overset{\mathrm{}}{}}n_k(1+n_k)\right]^{1/2}.$$
(11)
where $`n_k`$ is the occupancy of level $`k`$. As the temperature approaches zero, the particles condense to the ground state. Relation (11) thus yields $`(\delta N_0)N_0N`$. This result is unphysical since as $`T0`$, there is no energy to excite particles and hence the number fluctuation should vanish. This is a failure of the GCE in describing the fluctuation of an ideal bose gas at low temperatures. Experimentally, the system has a fixed total number of particles $`N`$ confined in space by a trapping potential, so that it is more accurately described by a canonical ensemble (CE), or better still, a microcanonical ensemble (MCE). This has been shown by several authors; some used property (8) for the MCE, and the ground state number fluctuation has been found to vanish as $`T0`$ . Within the GCE framework, however, the divergence can be removed by introducing inter-particle interactions . Consider an ideal gas in a box, and recall that in GCE the particle number fluctuation is related to the thermal compressibility via

$$\chi __T=\frac{1}{V}\left(\frac{V}{P}\right)=\frac{V}{T}\left(\frac{\delta N}{N}\right)^2.$$
(12)
For $`T0`$, the pressure is independent of volume, so that $`\chi __T\mathrm{}`$ and hence $`\delta N\mathrm{}`$. With a repulsive interaction, however weak, there is a pressure due to interaction to ensure the fluctuation to be finite. In the case of fermions, there exists a Pauli pressure even if the gas is ideal so that the ground state fluctuation finite :

$$(\delta N_0)^2^{1/2}=\left[\underset{k}{}n_k(1n_k)\right]^{1/2},$$
(13)
the sum $`k`$ runs over the ground state levels occupied at $`T=0`$. As $`T0`$, $`n_k1`$ for $`k`$ below fermi level, so that the fluctuation goes to zero as expected. In this paper we shall use the inverse square two-body interaction in one-dimension (the Calogero-Sutherland model, abbreviated CSM) . This interacting model may be mapped onto the Haldane-Wu generalized exclusion statistics (also known as fractional exclusion statistics, FES) whose quasiparticles are non-interacting, characterized by a parameter $`g`$ . $`g=0,1`$ correspond to bosons, fermions respectively. Other values of $`g`$ represent various degrees of ’Pauli blocking’. The model is solvable, and the energy spectrum is exactly known. In  we calculated the fluctuation for bosons and fermions using this interaction in both GCE and CE. For an interacting bose gas we showed that the divergence of the number fluctuation at $`T=0`$ in the GCE is removed. The ground state number fluctuations of an interacting gas using CE in two- and three-dimensions have been previously calculated by a number of authors for different types of interaction . Our choice of the inverse-square two-body interaction model in one-dimension has a two-fold advantage. First, it serves as an interacting model which allows us to demonstrate that the divergence in the GCE is removed when an interaction is introduced. Second, owing to its mapping to non-interacting FES particles, the model enables us to show that the divergence is removed even for an ideal gas provided that the Pauli blocking is not zero, no matter how weak. To the best of our knowledge, the exact fluctuation of interacting particles in the MCE has never been calculated before. In section III of this paper, we report some results of the exact fluctuation of trapped particles interacting via CSM, or equivalently, of ideal FES particles.

## II Ideal Fermions

Consider $`N`$ spinless fermions in a one-dimensional harmonic confinement. The ground state is defined at $`T=0`$ where the $`N`$ particles occupy the lowest energy levels up to and including the fermi level $`E_F`$. Given $`n`$ quanta of energy, consider breaking up $`n`$ into two parts:

$$n=n_h+n_p,$$
(14)
where $`n_h`$ is the number of quanta it takes to bring $`N_{ex}`$ particles to the fermi level $`E_F`$ (which is equivalent to the distribution of $`N_{ex}`$ holes to the states below and including $`E_F`$), and $`n_p`$ is the number of quanta it takes to distribute these $`N_{ex}`$ particles in the excited states (above $`E_F`$). This effectively divides the fermionic energy levels into two sectors: the particle space above $`E_F`$ and the hole space below and including $`E_F`$. Let $`\omega _h^{\{F\}}(n_h,N_{ex},N)`$ be the number of ways to distribute $`N_{ex}`$ holes in the hole space, and $`\omega _p^{\{F\}}(n_p,N_{ex},N)`$ be the number of ways to distribute $`N_{ex}`$ particles in the particle space, both according to Pauli Principle, then given $`n`$ quanta and $`N_{ex}`$ particles,

$$\omega ^F(n,N_{ex},N)=\underset{\{n_h,n_p\}}{}\omega _h^{\{F\}}(n_h,N_{ex},N)\omega _p^{\{F\}}(n_p,N_{ex},N),$$
(15)
where the set $`\{n_h,n_p\}`$ satisfies Eq. (14) for a given $`n`$. The problem now pertains to finding $`\omega _h^{\{F\}}(n_h,N_{ex},N)`$ and $`\omega _p(^{\{F\}}n_p,N_{ex},N)`$. At first glance this seems to be more complicated than finding a single quantity $`\omega ^F(n,N_{ex},N)`$. However, recall that $`\omega ^F(n,N_{ex},N)`$ is the number of ways to distribute $`N_{ex}`$ particles above $`E_F`$ with respect to a set of $`N`$ ground state energy levels. By breaking up the fermionic energy levels into two parts we are now distributing $`N_{ex}`$ particles above $`E_F`$ and $`N_{ex}`$ holes below $`E_F`$, both with respect to a single energy level. As we shall see shortly, this allows us to connect to number partitioning theory as in the case of bosons. Let us now look at these two spaces separately.

### A Particle space

First consider the particle space. This space is unbounded starting from the fermi level $`E_F`$. We now consider a situation in which the $`N_{ex}`$ particles have already been taken out of the hole space, and are residing in the fermi level $`E_F`$ awaiting to be distributed in the particle space. Note that $`n_p`$ quanta are distributed among $`N_{ex}`$ fermions, with respect to only one energy level $`E_F`$. The problem is now similar to the bosonic one, except the distribution of particles must comply with the Pauli Principle. In terms of number partitioning theory, $`\omega _p^F(n_p,N_{ex},N)`$ is the number of ways of partitioning an integer $`n_p`$ into $`N_{ex}`$ different parts, with $`N_{ex}N`$. As in the case of bosons, it is related to the canonical multiplicity via :

$`\omega _p^{\{F\}}(n_p,N_{ex},N)`$ $`=`$ $`\mathrm{\Omega }(n_p_p,N_{ex}),`$ (16)
$`=`$ $`\mathrm{\Omega }(n_p{\displaystyle \frac{N_{ex}(N_{ex}+1)}{2}},N_{ex}).`$ (17)
Remarkably, this is the same form as Eq. (8) for bosons with the shifted energy $`N_{ex}`$ replaced by $`_p`$. This shifted energy, $`_p=N_{ex}(N_{ex}+1)/2`$, is in fact the minimum energy it takes to excite $`N_{ex}`$ particles from $`E_F`$,

$$n_p^{min}(N_{ex})=\frac{(N_{ex}+1)N_{ex}}{2}.$$
(18)
Note that the partition function of this space is no longer given by that of fermions. The canonical multiplicity in Eq. (17) should therefore be thought of as the bosonic multiplicity. In other words, the problem is now mapped onto a similar bosonic problem, with the restriction that the parts of an integer being partitioned are different. This notion is most helpful when we discuss the hole space.

### B Hole space

We now consider taking $`N_{ex}`$ particles out of the multi-level ground state and put them in the fermi level $`E_F`$ (or equivalently, creating $`N_{ex}`$ holes in the ground state). For a given $`n_h`$ quanta, we wish to find $`\omega _h^{\{F\}}(n_h,N_{ex},N)`$, the number of ways of doing this. Unlike the particle space, the dimension of the hole space is bounded, set by the number of particle $`N`$ of the system. Note that for a given number $`N_{ex}`$ of particles, the Hilbert space of available states for $`N_{ex}`$ holes is dependent on the value of $`N_{ex}`$ itself and is given by:

$$N_H=NN_{ex}.$$
(19)
We need to find the partition function of this space for each $`N_{ex}`$, and derive a formula similar to Eq. (17) for $`\omega _h^{\{F\}}(n_h,N_{ex},N)`$. This may be done by considering a new system containing $`N_{ex}`$ bosons, whose energy space is bounded and is given by $`N_H+1`$ including the ground state. The goal is to determine $`Z_{N_{ex}}^H(x)`$, and expand this in terms of the coefficient $`\mathrm{\Omega }_h(n_h,N_{ex})`$. The $`N_{ex}`$-hole partition function of this hypothetical system can be found using :

$$Z_{N_{ex}}^H(x)=\frac{1}{N_{ex}}\underset{j=1}{\overset{N_{ex}}{}}Z_1^H(jx)Z_{N_{ex}j}^H(\beta ),$$
(20)
where $`Z_1^H(x)`$ is the one-particle partition function of the system containing $`N_{ex}`$ bosons and is given by:

$$Z_1^H(x)=\underset{i=0}{\overset{N_H}{}}x^i.$$
(21)
Note that the one-particle partition function needs to be determined for a given $`N_{ex}`$. Using this, $`Z_{N_{ex}}^H(x)`$ may then be found from (20). Once $`Z_{N_{ex}}^H(x)`$ is found, we may expand it in powers of $`x`$:

$$Z_{N_{ex}}^H(x)=\underset{i}{}\mathrm{\Omega }_h(i,N_{ex})x^i.$$
(22)
We are now ready to determine a formula for $`\omega _h^{\{F\}}(n_h,N_{ex},N)`$. Because the hole space includes $`E_F`$, the minimum energy to create a hole (or dig a particle) to put in the fermi level is zero, since there already is a particle there; for two holes the minimum energy is one, for three holes it is three…etc. In general,

$$n_h^{min}(N_{ex})=\frac{(N_{ex}1)N_{ex}}{2}.$$
(23)
Similar to Eq. (17) with the energy shift $`_h`$ given by (23), the number of ways of creating $`N_{ex}`$ holes in the hole space $`\omega _h^F(n_h,N_{ex},N)`$ is given by:

$$\omega _h(n_h,N_{ex},N)=\mathrm{\Omega }_h(n_h\frac{N_{ex}(N_{ex}1)}{2},N_{ex}).$$
(24)
Using Eqs. (17) and (24), the number of ways of distributing $`N_{ex}`$ fermions to the excited states, Eq. (15), now reads:

$$\omega ^F(n,N_{ex},N)=\underset{\{n_h,n_p\}}{}\mathrm{\Omega }(n_p\frac{N_{ex}(N_{ex}+1)}{2},N_{ex})\mathrm{\Omega }_h(n_h\frac{N_{ex}(N_{ex}1)}{2},N_{ex}).$$
(25)
It is obvious that if there is no hole space, $`\omega _h(n_h,N_{ex},N)=\mathrm{\Omega }_h(n_h_h,N_{ex})=1`$, the sum over the set $`\{n_h,n_p\}`$ vanishes since $`n=n_p`$, the energy shift $`_p=n_p^{min}(N_{ex})=N_{ex}`$, and Eq. (15) reduces to Eq. (8) for bosons.

Using Eqs. (25) and (2)-(6), we calculated the ground state number fluctuation of fermions for $`N=100`$. The result is shown in Fig. 4. For comparison we also show the corresponding result in the CE.

## III Ideal FES particles

We now calculate the exact ground state fluctuation of particles interacting via CSM, or equivalently, non-interacting particles obeying FES . The Hamiltonian of the CSM is given by :

$$H=\underset{i=1}{\overset{N}{}}\left[\frac{\mathrm{}^2}{2m}\frac{^2}{x_i^2}+\frac{1}{2}m\omega ^2x_i^2\right]+\frac{\mathrm{}^2}{m}\underset{i<j=1}{\overset{N}{}}\frac{g(g1)}{(x_ix_j)^2}$$
(26)
with the dimensionless coupling parameter $`g0`$. The particles are confined in a harmonic well and the thermodynamic limit is obtained by taking $`\omega 0`$ as $`N\mathrm{}`$, with $`\omega N=constant`$. In the thermodynamic limit, the properties of the system are translationally invariant, and would be the same if the particles were on a line, or a circle, instead of a harmonic confinement. To make the problem well-defined quantum mechanically, we have to demand that the wave functions go to zero as $`|x_ix_j|^g`$ whenever two particles i and j approach each other. Since the particles cannot cross each other, we may choose the wave function to be either symmetric (bosonic) or antisymmetric (fermionic). For $`g=0`$ and $`1`$, the model describes free bosons and free fermions respectively.

In  we have calculated the fluctuation for FES particles in the GCE and CE. To do this in the MCE, we need to find $`\omega ^{\{g\}}(n,N_{ex},N)`$ for a given $`g`$. The energy spectrum of the CSM Hamiltonian (26) is exactly known, and may be expressed in terms of the quasiparticle energy spectrum $`ϵ_g`$ of non-interacting FES particles, which is given by :

$$ϵ_g=ϵ_k(1g)N_k,$$
(27)
with $`\mathrm{}\omega 1`$. $`N_k`$ is defined as the number of particles below energy level $`k`$, $`ϵ_k=(k\frac{1}{2}),k=1,2,\mathrm{}`$ is the harmonic oscillator energy spectrum. The canonical N-particle partition function of FES particles reads:

$`Z_N`$ $`=`$ $`x^{E_N^{\{g\}}(0)}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{(1x^j)}},`$ (28)
$`=`$ $`x^{E_N^{\{g\}}(0)}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{\Omega }(n,N)x^n.`$ (29)
Note that the canonical multiplicity $`\mathrm{\Omega }(n,N)`$ is still the same for FES particles as for bosons and fermions. The only effect of the interaction strength $`g`$ is to alter the overall ground state energy $`E_N^{\{g\}}(0)=gN(N1)/2+N/2`$ which reduces to that of fermions, bosons for $`g=1,0`$ respectively.

Using Eq. (27), one may determine the spectra of FES particles for any value of $`g`$. In the Fig. (4) we show an example of how to find the spectrum of semions ($`g=1/2`$) for some values of excitation quanta $`n`$ and number of particles $`N=5`$. Starting from the fermionic spectrum on the left, for each value of $`n`$, the spectrum of semions is drawn and the resulting particles which are in the excited states are determined. In the case of $`n=3`$, for instance, the only possibility for fermions is to excite one particle and there are 3 ways of doing this as shown. Therefore, $`\omega ^{\{F\}}(3,1,5)=3`$, and $`\omega ^{\{F\}}(3,2,5)=\omega ^{\{F\}}(3,3,5)=0`$ since $`3`$ quanta is too few to excite $`2`$ or more particles. For semions, however, $`\omega ^{\{1/2\}}(3,1,5)=1`$, $`\omega ^{\{1/2\}}(3,2,5)=2`$, and $`\omega ^{\{1/2\}}(3,3,5)=0`$. Note that $`\mathrm{\Omega }(3,5)=3`$ in both cases. Following the same procedure, we found the microcanonical multiplicities $`\omega ^{\{g\}}(n,N_{ex},N)`$ for $`g=3/4,1/2,1/4`$, $`n=1\mathrm{}16`$ quanta, and $`N=2\mathrm{}5`$ particles. The method gets more cumbersome, however, for larger values of $`n`$ and $`N`$. Ideally, one wishes to be able to determine $`\omega ^{\{g\}}(n,N_{ex},N)`$ from $`\omega ^{\{F\}}(n,N_{ex},N)`$ or $`\omega ^{\{B\}}(n,N_{ex},N)`$, or from $`\mathrm{\Omega }(n,N_{ex})`$ similar to Eqs. (8) and (25). This general formula for $`\omega ^{\{g\}}(n,N_{ex},N)`$ for any $`g`$, if it exists, is yet to be found. Here, we report the finding of the microcanonical multiplicities for only two values of $`g`$: $`g=\frac{N2}{N1}`$ (close to fermions), and $`g=\frac{1}{N1}`$ (close to bosons). For $`g=\frac{N2}{N1}`$, $`\omega ^{\{\frac{N2}{N1}\}}(n,N_{ex},N)`$ is found to be given by:

$$\omega ^{\{\frac{N2}{N1}\}}(n,N_{ex},N)=\omega ^{\{F\}}(n+N,N_{ex},N)\omega ^{\{F\}}(n+N,N_{ex},N1),$$
(30)
For $`g=\frac{1}{N1}`$:

$`\omega ^{\{\frac{1}{N1}\}}(n,N_{ex},N)\text{ }`$ $`=`$ $`\omega ^{\{B\}}(n,N_{ex},N),\text{ }N_{ex}N1,N`$ (31)
$`\omega ^{\{\frac{1}{N1}\}}(n,N1,N)`$ $`=`$ $`\omega ^{\{B\}}(n,N1,N)+\omega ^{\{B\}}(n1,N1,N),`$ (32)
$`\omega ^{\{\frac{1}{N1}\}}(n,N,N)\text{ }`$ $`=`$ $`\omega ^{\{B\}}(nN,N,N).`$ (33)
These multiplicities must, of course, satisfy

$$\mathrm{\Omega }(n,N)=\underset{N_{ex}=1}{\overset{N}{}}\omega ^{\{g\}}(n,N_{ex},N).$$
(34)
For $`N5`$ and $`n16`$, the values found using Eq. (30)-(33) were verified with those found using the direct combinatorial method as described above.

Although a general formula for the microcanonical multiplicities for any $`g`$ has not been found, an important point concerning them is observed. For a given number of particles $`N`$, consider a set of discrete values of $`g`$,

$$g=1,\frac{N2}{N1},\frac{N3}{N1},\mathrm{},\frac{1}{N1},0.$$
(35)
Since the levels of the FES particles are shifted by an amount which depends on the value of $`g`$ and the number of particles below (see Eq. (27)), for some values of $`g`$ a particle might lie very close to the last level of the ground state, $`E_F^{\{g\}}`$. However, the particles are considered excited if and only if they lie above $`E_F^{\{g\}}`$, no matter how close. This results in the multiplicities $`\omega ^{\{g\}}(n,N_{ex},N)`$, for $`\frac{Ni}{N1}>g\frac{N(i+1)}{N1},i=1,\mathrm{},N1`$ to be the same. Note the equal sign in $`g\frac{N(i+1)}{N1}`$. For instance,

$`\omega ^{\{g\}}(n,N_{ex},N)`$ $`=`$ $`\omega ^{\{B\}}(n,N_{ex},N),for{\displaystyle \frac{1}{N1}}>g0.`$ (36)
Note that for these values of $`g`$ given by Eq. (35), the ’Fermi level’ $`E_F^{\{g\}}=N(1g)(N1)`$ is integral. So for $`g=(Nk)/(N1)`$, where $`k`$ is an integer, the ’Fermi level’ lines up with the $`(k1)^{th}`$ level below $`E_F`$ of fermions. In box 1, $`n=0`$ of Fig. 4, for instance, where $`N=5`$ and $`g=2/4`$ ($`k=3`$), the ’Fermi level’ lines up with the second level below $`E_F`$. This explains why the microcanonical multiplicities for some range of $`g`$ such that the ’Fermi level’ lies between the $`k^{th}`$ and $`(k+1)^{th}`$ levels of fermions are the same. This is due the discrete nature of the energy levels.

## IV discussion

In our previous work (Ref. ) where the direct combinatorial method was used, the fermionic calculation of the ground state fluctuation was restricted to a low number of particles $`N`$ and quanta $`n`$ using a normal office computer (Pentium III, 500 cpu). For a relatively small number of particles (e.g., $`N=10`$), at higher excitation $`n`$ the combinatorics method is more time-consuming due to the rapid increase in the number of possibilities with $`n`$. The method described in section II translates the problem in combinatorics into the problem in calculating the partition functions of the hole space, the latter being simpler computationally. Although this method is still time-consuming and the calculation for larger $`N`$ ($`N100`$) is still not possible using our office computer, it is more effective for higher number of quanta and relatively small number of particles. For demonstration we display the ground state fluctuation of fermions as a function of energy quanta $`n`$ in Fig. 4 for $`N=30`$. We also show the corresponding curve in CE, which is the same as GCE except at very low temperature (see e.g., Ref. ) for comparison. As expected, both go to zero as $`T0`$, with the MCE fluctuation less than that given by CE for all $`n`$. Note that the two fluctuations are very different, even for very high excitations. At $`n=6000`$ which is $`200\times E_F`$, the CE curve still differs from the MCE one by about $`14\%`$. It was shown in Ref.  for $`N=15`$ that the ground state occupancy $`N_0=_k^{E_F}n_k`$ for the two ensembles are very similar. Clearly, the number fluctuation is more sensitive to the ensemble used. Therefore, for a relatively small particle number, while it may be adequate to use the CE (or GCE) to describe a thermodynamic quantity such as the ground state occupancy, it should be used with caution when calculating the number fluctuation and related quantities.

Using Eqs. (30), (33) and (2)-(6), we calculated the exact ground state number fluctuation for interacting particles. In Fig. 4a we display the fluctuation for $`N=5`$, $`g=3/4,1/4`$. For comparison we also show the fluctuation of free fermions and bosons. The values for $`g=3/4,1/4`$ obtained using the direct combinatorial method (see Fig. 4) for $`n16`$ are also shown. Note that the curves agree with these values exactly. Note also that the number fluctuation for free fermions and free bosons cross at a certain energy, with the fermionic one starting from smaller at small quanta, to larger at high quanta. This is because the number of possibilities of creating holes within the fermi sea and distributing particles above, which starts from low at small energy, increases more rapidly than for bosons whose ground state consists of only one level. Similar behaviours are observed for FES particles, whose $`g`$ values ($`g=3/4,1/4`$) represent partial Pauli blocking which are both less than that of fermions.

Eqs. (30) and (33) in principle may be applied for any $`N`$. However, since they involve the addition of two quantities which may be very large at large quanta and number of particles, there is a difficulty in obtaining their values accurately. Therefore, without lost of accuracy, we restricted the calculations to $`N=10`$. Fig. 4b shows the ground state number fluctuation for $`N=10`$, $`g=1,8/9,1/9,0`$. Note that the curves for $`g=8/9,1/9`$ are closer to those of fermions and bosons than for $`g=3/4,1/4`$, $`N=5`$, in Fig. 4a. As $`N`$ gets larger we expect the fluctuation curves for $`g=\frac{N2}{N1}`$ and $`\frac{1}{N1}`$ to come very close to those of free femions and free bosons. The formulae (30) and (33) are therefore useful only for systems with a small number of particles. Note also that in both graphs, the results for $`g=\frac{N2}{N1}`$ are closer to those of fermions than for $`g=\frac{1}{N1}`$ to bosons. This is can be understood from comparing Eqs. (30) and (33). Eq. (30) involves the difference of the fermionic microcanonical multiplicities of two system sizes, whose values may be very similar especially at low quanta. Eq. (33), however, involves the addition of two bosonic microcanonical multiplicities. This brings the microcanonical multiplicities of $`g=\frac{1}{N1}`$ further from those of bosons than $`g=\frac{N2}{N1}`$ from those of fermions.

Finally, In Figs. 4a,b we compare these exact fluctuation curves with the corresponding CE curves for $`g=\frac{N2}{N1},\frac{1}{N1}`$, $`N=5,10`$. From , the expression for the canonical ground state fluctuation is given by:

$$(\delta N_0)^2^{\{g\}}=g((\delta N_0)^2)^F+(1g)((\delta N_0)^2)^B$$
(37)
Clearly, the results between the two ensembles are very different especially for small $`n`$. In both graphs the microcanonical fluctuations are less than the canonical ones for all $`n`$. Note that the CE formula (37) allows one to calculate the fluctuation for a general $`g`$. The determination of a general formula for the microcanonical multiplicity, $`\omega ^{\{g\}}(n,N_{ex},N)`$, and hence the ground state number fluctuation for any $`g`$, if it exists, remains a challenge.

## V acknowledgments

MNT would like to acknowledge the Natural Sciences and Engineering Research Council (NSERC) of Canada for financial assistance. Special thanks to Dr. J. C. Waddington for suggesting an important point regarding ideal fermions, Dr. D. Sen for going through the exact combinatorial results of FES particles, and Dr’s. M. V. N. Murthy, R. K. Bhaduri for guidance and a careful reading of the manuscript.