text:
# Atom-Molecule Laser Fed by Stimulated Three-Body Recombination

On leave from the \]Institute for Space Sciences, Bucharest-Magurele 76900, Romania.

(Submitted May 14, 2003)

## Abstract

Using three-body recombination as the underlying process, we propose a method of coherently driving an atomic Bose-Einstein condensate (BEC) into a molecular BEC. Superradiant-like stimulation favors atom-to-molecule transitions when two atomic BECs collide at a resonant kinetic energy, the result being two molecular BEC clouds moving with well defined velocities. Potential applications include the construction of a molecule laser.

One of the major goals of the recent research in BEC physics is the creation of a molecular condensate zoller2002 . Many applications are envisioned for such a system, including coherent quantum chemistry moore2002 . One approach to forming a molecular BEC utilizes an atomic condensate along with a method of coherently driving pairs of atoms into (bound) molecular states. Two mechanisms have been proposed for this: (i) photoassociation, in which an external laser couples atomic states to molecular states juli98 , and (ii) control of the scattering length via external magnetic field ramps tuned near a Feshbach resonance timm99 . Another problem of interest to the BEC community is the realization of an atom laser hagley99 . A proposed technique kozuma99  for this is the application of a time-dependent light pulse resulting in Bragg scattering of the condensate cloud, in which a fraction of the condensate receives a large momentum kick and attains velocity relative to the parent cloud.

In this Letter we propose a new mechanism for the formation of a molecular BEC from an atomic BEC. This mechanism relies on a twofold bosonic stimulation of the three-body recombination process — both by the number of atoms present and by the number of molecules. We show that two atomic condensates with a well-defined (resonant) relative velocity will mutually stimulate transitions into molecular states, resulting in the production of two coherent molecular clouds, each with a well-defined momentum.

Consider first the elementary process involving only three isolated atoms. The atoms, initially in an unbound state, are converted through recombination into a bound molecular state with $`1/3`$ of the recoil energy, while a third atom carries away the remaining $`2/3`$. Conservation of momentum requires that the resulting atom and molecule have momenta of equal magnitude ($`\mathrm{}k_0`$), but in opposite directions. The value of $`k_0`$ is fixed by conservation of energy to be $`k_0=\sqrt{4mE_0/3\mathrm{}^2}`$, where $`E_0`$ is the molecular binding energy and $`m`$ is the atomic mass. This process is characterized by a transition matrix element $`T_{if}`$ between the initial and final three-body states. The spontaneous decay of an ultracold atomic sample of density $`n`$ through three-body recombination is usually characterized by a loss rate coefficient $`K_3`$, such that $`dn/dt=3K_3n^3/6`$. Here $`n^3/6`$ is the number of atomic triplets available for recombination in the considered volume and the factor of 3 denotes the atoms lost in each recombination. An additional factor of 1/6 multiplies the right-hand side if all three atoms share the same state, i.e., if we consider the decay of a pure condensate kagan85 . The connection between $`T_{if}`$ and $`K_3`$ is given by fed96

$`K_3={\displaystyle \frac{6/\mathrm{}}{(2\pi )^2}}{\displaystyle \frac{d^3𝐤}{dE}\left|T_{if}\right|^2\delta \left(E\frac{3\mathrm{}^2k_0^2}{4m}\right)𝑑E}={\displaystyle \frac{6k_0m}{\mathrm{}^3}}|T_{if}|^2.`$
If the recombination forms a weakly bound level with $`E_0\mathrm{}^2/(ma^2)`$ then $`k_0=2/(\sqrt{3}a)`$, where $`a`$ is the $`s`$-wave scattering length. In this case, working in the ultracold limit and using momentum normalized states, Ref. obtains $`T_{if}=\zeta \mathrm{}^2a^{5/2}/m`$ fed96 , where $`\zeta `$ is a numerical constant. More generally $`K_3`$ has been shown to scale with the fourth power of the scattering length in most cases, $`K_3=\beta (\mathrm{}/m)a^4`$ esry99 ; nielsen99  which for the weakly bound level gives $`\beta =(12\zeta ^2)/(\pi \sqrt{3})`$. Ref. fed96  finds $`\beta 23`$, while a more rigorous calculation, valid for a wide range of interatomic potentials esry99 , yields $`\beta 180`$ for $`a>0`$ and $`\beta 1014`$ for $`a<0`$. (We omit a $`\mathrm{sin}^2`$ modulation factor discussed in .)

We are interested in the process in which three-body recombination is stimulated by bosonic enhancement. We consider two atomic condensate clouds moving with relative velocity $`\mathrm{}k_0/m`$, and choose a reference frame in which the two clouds move with equal velocities towards each other (see Fig. 1) along the direction $`\widehat{k}`$. Let the wavevector of the atoms moving to the right (left) be denoted by $`𝐤_1`$ ($`𝐤_2`$) and the population of the mode by $`n_1`$ ($`n_2`$). Note that $`𝐤_1=𝐤_2=k_0/2\widehat{k}`$. Recombination of three atoms with momentum $`\mathrm{}𝐤_1`$ can lead to a final state in which the resulting atom has momentum $`\mathrm{}𝐤_2`$, in which case the resulting molecule has the momentum $`\mathrm{}𝐤_{\text{m1}}=2\mathrm{}k_0\widehat{k}`$. This transition will be strongly favored due to bosonic stimulation if the mode $`𝐤_2`$ is highly populated. We recall the fact that the magnitude of the momentum is fixed by energy and momentum conservation and only its direction has an arbitrary value. Therefore, it is possible that the atom-molecule pair is ejected along a direction different from $`\widehat{k}`$. However, this process does not benefit from bosonic stimulation. By symmetry, recombination of three atoms with momentum $`\mathrm{}𝐤_2`$ will be strongly enhanced if the resulting atom-molecule pair is ejected along the direction $`\widehat{k}`$. In this case the momentum of the molecule has the value $`\mathrm{}𝐤_{\text{m}2}=2\mathrm{}k_0\widehat{k}`$ while the resulting atom has momentum $`\mathrm{}𝐤_1`$. We denote the numbers of molecules with wavevectors $`𝐤_{\text{m}1}`$ and $`𝐤_{\text{m}2}`$ by $`n_{\text{m}1}`$ and $`n_{\text{m}2}`$, respectively. Finally, note that, for any three-body recombination the reverse process is also possible. For example, an atom with momentum $`\mathrm{}𝐤_2`$ that collides with a molecule with momentum $`\mathrm{}𝐤_{\text{m}1}`$ can induce it to break up, resulting in three atoms with momentum $`\mathrm{}𝐤_1`$. In summary, the four modes considered here are coupled by the transitions

$$\begin{array}{cc}\hfill 3\text{A}(𝐤_1)& \text{A}(𝐤_2)+\text{A}_2(𝐤_{\text{m}1}),\hfill \\ \hfill 3\text{A}(𝐤_2)& \text{A}(𝐤_1)+\text{A}_2(𝐤_{\text{m}2}).\hfill \end{array}$$
(1)
The many-body interaction Hamiltonian describing the reactions in Eqs. (1) is

$$\widehat{H}_i=\mathrm{}\omega \left(c_2^{}c_{\text{m}1}^{}c_1c_1c_1+c_1^{}c_{\text{m}2}^{}c_2c_2c_2\right)+\text{h.c.}$$
(2)
Here $`c_{1,2}`$ ($`c_{\text{m}1,\text{m}2}`$) are the annihilation operators for the two atomic (molecular) modes and the coupling frequency is given by $`\mathrm{}\omega =T_{if}/V^{3/2}`$, where $`V`$ is the quantization volume. For a weakly bound molecule $`\omega =[\zeta \mathrm{}/(ma^2)](a^3/V)^{3/2}`$. The two operators,

$$\widehat{N}=c_1^{}c_1+c_2^{}c_2+2c_{\text{m}1}^{}c_{\text{m}1}+2c_{\text{m}2}^{}c_{\text{m}2},$$

$$\frac{\widehat{P}}{(\mathrm{}k_0/2)}=c_1^{}c_1c_2^{}c_2+4c_{\text{m}1}^{}c_{\text{m}1}4c_{\text{m}2}^{}c_{\text{m}2},$$

commute with the Hamiltonian in Eq. (2) as a consequence of the conservation of the total number of atoms and of the total momentum.

We now derive the equations of motion for the creation and annihilation operators in the Heisenberg representation. Adopting the notation $`a_{1,2}(\tau )`$ and $`b_{1,2}(\tau )`$ for the atomic and molecular Heisenberg operators divided by $`\sqrt{N}`$ (in which fast oscillations with the frequency associated to the non-interacting energy were factored out), the equations of motion read,

$`i{\displaystyle \frac{d}{d\tau }}a_1`$ $`=3a_1^{}a_1^{}b_1a_2+b_2^{}a_2a_2a_2,`$ (3)
$`i{\displaystyle \frac{d}{d\tau }}a_2`$ $`=3a_2^{}a_2^{}b_2a_1+b_1^{}a_1a_1a_1,`$ (4)
$`i{\displaystyle \frac{d}{d\tau }}b_1`$ $`=a_2^{}a_1a_1a_1,`$ (5)
$`i{\displaystyle \frac{d}{d\tau }}b_2`$ $`=a_1^{}a_2a_2a_2.`$ (6)
We use a scaled time $`\tau =\mathrm{\Omega }t`$ with $`\mathrm{\Omega }=N^{3/2}\omega =n^{3/2}T_{if}`$ ($`n=N/V`$ is the total concentration of atoms). For a shallow molecular level $`\mathrm{\Omega }=[\zeta \mathrm{}/(ma^2)](na^3)^{3/2}`$. Equations (3)–(6) describe the dynamic population transfer between the atomic and molecular modes due to resonant three-body recombination. They are based on the assumption that all two-body elastic collisions (atom-atom, molecule-molecule and atom-molecule) and also collisional loss of the weakly bound molecules to deeper molecular states can be neglected. In addition, we have neglected the spatial and temporal changes of the slowly varying envelopes of the condensate clouds. The simplest approximation to the above equation consists of replacing the operators by complex numbers equal to their averages: $`\stackrel{~}{a}_{1,2}(t)=a_{1,2}(t)`$ and $`\stackrel{~}{b}_{1,2}(t)=b_{1,2}(t)`$, i.e., the classical field approximation. The resulting system of ordinary differential equations, with initial conditions, can be integrated numerically. We consider the experimentally relevant case of an initial state with zero molecules. We also assume that initially the $`N`$ atoms are evenly split between the two momentum eigenstates: $`N_{1,2}n_{1,2}(t=0)=N/2`$. Assuming real initial conditions, we set $`\stackrel{~}{b}_{1,2}(0)=0`$ and $`\stackrel{~}{a}_{1,2}(0)=1/\sqrt{2}`$. The approximate version of the system of Eqs. (3)–(6) can be solved to describe the population transfer from the atomic to the molecular modes in a process characterized by the time constant $`T_R=2\pi /\mathrm{\Omega }`$. A characteristic that is similar to the superradiance phenomenon is double amplification (stimulation): after an initial buildup period for the molecular population, rates of further transitions are proportional both to the number of existing molecules in the recombined momentum state and to the number of atoms in their relevant momentum state.

The question of the validity of using the classical field approximation can be studied in general by accounting for higher-order terms in a cluster expansion of the four-operator products in Eqs. (3)–(6) (see, e.g., Ref. kohler02 ). However it is sufficient for our purpose to compare the results of the approximate version of Eqs. (3)–(6) to the many-body Schrödinger equation written using a basis of Fock states. We designate these basis states by $`|n_1,n_2,n_{\text{m}1},n_{\text{m}2}`$ where $`n_{1,2}`$ ($`n_{\text{m}1,\text{m},2}`$) are the occupation numbers of the two atomic (molecular) modes. Due to the two conservation laws $`n_1=N_1+n_{\text{m}2}3n_{\text{m}1}`$ and $`n_2=N_2+n_{\text{m}1}3n_{\text{m}2}`$, we can restrict the basis set to only those accessible by three-body recombinations from our initial state and thus use only two occupation numbers, say $`n_{\text{m}1}`$ and $`n_{\text{m}2}`$, to parameterize the Fock states. With this notation, the many-body wavefunction is

$$|\mathrm{\Psi }(t)=\underset{n_{\text{m}1},n_{\text{m}2}}{}f_{n_{\text{m}1},n_{\text{m}2}}(t)|n_{\text{m}1},n_{\text{m}2}.$$
(7)
After eliminating the fast oscillations and introducing the time scaling $`\tau ^{}=\omega t`$, the Schrödinger equation becomes

$`i{\displaystyle \frac{d}{d\tau ^{}}}f_{n_{\text{m}1},n_{\text{m}2}}`$ (8)
$`=\sqrt{(n_1+1)(n_1+2)(n_1+3)n_2n_{\text{m}1}}f_{n_{\text{m}1}1,n_{\text{m}2}}`$
$`+\sqrt{n_1(n_11)(n_12)(n_2+1)(n_{\text{m}1}+1)}f_{n_{\text{m}1}+1,n_{\text{m}2}}`$
$`+\sqrt{(n_2+1)(n_2+2)(n_2+3)n_1n_{\text{m}2}}f_{n_{\text{m}1},n_{\text{m}2}1}`$
$`+\sqrt{n_2(n_21)(n_22)(n_1+1)(n_{\text{m}2}+1)}f_{n_{\text{m}1},n_{\text{m}2}+1}.`$
Solving these equations with the initial conditions $`f_{0,0}(0)=1`$ and $`f_{n_{\text{m}1},n_{\text{m}2}}(0)=0`$ for $`n_{\text{m}1}0`$ and $`n_{\text{m}2}0`$, allows us to calculate the average numbers of atoms populating any of the four considered modes. These averages are shown in Fig. 2 for different values of $`N_1=N_2=N/2`$. Comparison with the corresponding results obtained using Eqs. (3)–(6) shows agreement over the initial period of fast transitions of atoms into molecular states. As expected, the classical field approximation fails at shorter times for lower values of $`N`$, while for sufficiently large $`N`$ it can be used to accurately describe the initial transition burst converting atoms to molecules. Equation (8) allows us to describe the population transfer within the four-mode approximation without making any additional approximations. However, calculations using the Fock basis become prohibitive if more than a few thousand atoms are considered.

Experimental arrangements that use the process described above to produce coherent molecular clouds can be readily envisioned. The simplest experiment would begin with a time-dependent Bragg splitting of a single atomic cloud as in the experiment of Kozuma *et al.* kozuma99 , but with the particular choice of momentum kick that matches the three-body resonance condition. In addition, the value of the molecular binding energy can be tuned using an external magnetic field inouye98 ; Donley02 . In order to take full advantage of the stimulated transition into molecular states, (i.e., to transform a large fraction of the atomic clouds into molecular clouds) the collision time of the two atomic clouds should be comparable to the characteristic time of the three-body recombination process $`t_{\text{coll}}T_R`$. Here $`t_{\text{coll}}R_c/(\mathrm{}k_0/m)`$ and $`R_c`$ is the characteristic size of the condensate along the collision direction. Consider a sample of 85Rb with a concentration of $`10^{15}/\mathrm{cm}^3`$. Assuming that the scattering length is adjusted with the help of a Feshbach resonance to $`a=400`$ a.u. and that $`E_0=\mathrm{}^2/ma^2`$, then we estimate $`T_R=0.45`$ ms (using $`K_3`$ given by esry99 ). The resonant process described here is characterized by a finite width if additional processes that limit the lifetime of the involved states are considered. One such process is the collsionally induced decay of the molecular level to lower molecular states. Only a very rough estimation of the lifetime is available for molecules formed by atomic species commonly used in BEC experiments. Some estimations (see, e.g, timm99 ) cite values of at least 1 ms for this lifetime at molecular densities of $`10^{15}/\mathrm{cm}^3`$. In addition, recent experiments in 85Rb Donley02  suggest that this is a conservative estimation, and that the actual lifetime might be somehow longer.

Up to this point we have neglected elastic collisions between atoms and molecules. In fact, if the resonant value of the relative velocity is not enforced by the experimental setup, depletion of the atomic clouds by elastic atom-atom scattering is the main loss process taking place before the two condensate clouds separate band2000 . As discussed in Ref. band2000 , elastic scattering must be accounted for beyond the mean-field (i.e., Gross-Pitaevskii) approximation in order to observe the condensate depletion due to elastic scattering. We can include this process in our simple model as a loss term in the interaction Hamiltonian. Considering only atom-atom collisions involving one atom with wavevector $`𝐤_1`$ and one with $`𝐤_2`$, Eqs. (3) and (4) in the classical field approximation become

$$\begin{array}{cc}\hfill i\frac{d}{d\tau }\stackrel{~}{a}_1& =3\stackrel{~}{a}_1^{}\stackrel{~}{a}_1^{}\stackrel{~}{b}_1\stackrel{~}{a}_2+\stackrel{~}{b}_2^{}\stackrel{~}{a}_2\stackrel{~}{a}_2\stackrel{~}{a}_2i\gamma |\stackrel{~}{a}_2|^2\stackrel{~}{a}_1,\hfill \\ \hfill i\frac{d}{d\tau }\stackrel{~}{a}_2& =3\stackrel{~}{a}_2^{}\stackrel{~}{a}_2^{}\stackrel{~}{b}_2\stackrel{~}{a}_1+\stackrel{~}{b}_1^{}\stackrel{~}{a}_1\stackrel{~}{a}_1\stackrel{~}{a}_1i\gamma |\stackrel{~}{a}_1|^2\stackrel{~}{a}_2.\hfill \end{array}$$
(9)
Here the parameter $`\gamma `$ is the elastic scattering rate constant divided by 2$`\mathrm{}\mathrm{\Omega }`$, and is given by $`\gamma =(\mathrm{}k_0/m)\sigma n/(2\mathrm{}\mathrm{\Omega })`$ band2000 , where $`\sigma `$ is the elastic scattering cross section at the relevant two-body collision energy $`E_{\text{coll}}=\mathrm{}^2k_0^2/m`$. Eqs. (9) show that the two processes — resonant three-body recombination and condensate loss due to elastic scattering — compete with each other and that the value of the parameter $`\gamma `$ is crucial in determining which process dominates. Fig. 3 shows the solution of the Heisenberg equations that include the elastic loss terms for two extreme cases. The example with $`\gamma =0.1`$ illustrates the case of little elastic loss allowing almost $`70`$% of the atoms to be transformed into molecules if $`t_{\text{coll}}T_R`$. However, if $`\gamma =1.5`$, for the same collision time more than $`80`$% of the atoms are lost due to elastic collisions and only about 8% end up in the molecular clouds.

The simplest estimation of $`\gamma `$ can be attempted in the low energy limit ($`k_0a^1`$) where $`\sigma 8\pi a^2`$. This leads to $`\gamma =(8\pi /\sqrt{3}\zeta )(na^3)^{1/2}`$ and yields high values of $`\gamma `$ even for very large values of the diluteness parameter (i.e., $`na^3=1`$ corresponds to $`\gamma =1.5`$). However, note the importance of using an exact value of $`\sigma `$ corresponding to the correct $`E_{\text{coll}}`$. For example, in the case of a weakly bound molecular level \[$`E_0\mathrm{}^2/(ma^2)`$\], $`k_01/a`$, which is beyond the validity of the low energy approximation. Apparently, changing the value of the parameter $`na^3`$ by changing either $`n`$ or $`a`$ is the primary way of shifting the balance between the two processes. However, special conditions might enhance or suppress either process. One such situation may occur when $`a`$ approaches zero and therefore elastic scattering vanishes while three-body recombination still has a considerable probability. Calculations by Esry et al. esry99  predict a significant $`K_3`$ for $`a=0`$. Similarly, Stenger et al. stenger99  measured three-body recombination loss from Na condensates in the vicinity of Feshbach resonances, observing a large amount of loss even at small values of $`a`$. Consider the Feshbach resonance located near $`B=535.7`$G in 85Rb. Tuning the magnetic field such that $`a=0`$ ($`B=537.6`$ G) leads to $`E_0=75\mu `$K, and for a density of $`10^{16}/\mathrm{cm}^3`$ we obtain $`K_3=5.4\times 10^{32}`$ cm6/s (as estimated in Ref. esry99 ), giving $`T_R28`$ ms. For Na with $`a=0`$ ($`B=907`$ G), Ref. stenger99  measures $`K_310^{28}`$ cm6/s, giving $`T_R0.24`$ ms. Alternatively, one can choose the experimental conditions such that the elastic collisions are strongly suppressed exactly at the resonant value of the collision energy $`E_{\text{coll}}`$ due to the presence of a Ramsauer minimum (see, for example, Burke98 ). The exact resonant matching between the position of the Ramsauer minimum and that of the binding energy of the molecular state (i.e., $`E_{\mathrm{coll}}=E_0/3`$) can be achieved near some Feshbach resonances. For the above mentioned resonance in 85Rb, this matching occurs for $`B=547`$ G, when $`a=353`$ a.u., and $`E_0=1.05`$ mK. In this case, considering $`n=10^{15}/\mathrm{cm}^3`$, one finds $`T_R1.3`$ ms ($`K_3=9.2\times 10^{26}`$cm6/s cf. esry99 ). Since producing an atomic condensate at negative $`a`$ might be a problem, we assume here that it is produced at some positive $`a`$ and that the magnetic field is swept into resonance precisely when the two atomic clouds start colliding.

In conclusion, resonant three-body recombination that benefits from the superradiant-like bosonic stimulation could be a promising candidate for producing coherent samples of molecules from existing atomic condensates. The simplest experimental setup that takes advantage of this process involves a collision of two atomic condensates, which produces two molecular condensates in counterpropagating momentum eigenstates. This setup could be particularly promising for the construction of a molecule laser by colliding two atomic condensates.

This work was supported in part by the NSF.


format: text

text:
# The Opacity of Nearby Galaxies from Counts of Background Galaxies: II. Limits of the Synthetic Field MethodBased on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

## 1 Introduction

### 1.1 Some history

How much dust is there in galaxy disks, and where is it? Reliable answers to these questions have remained elusive, although they bear directly on our understanding of many astronomical problems. Dust has to be properly accounted for to infer luminosity, stellar composition, and mass distribution correctly. Dust could be veiling correlations between global properties of galaxies and affecting our view of their dynamical evolution. It could be hiding significant amounts of luminous baryonic matter, or could trace extremely cold gas, and hence have implications for the dark matter problem. Dust in disks and halos of foreground spirals could be distorting our perception of the high-redshift universe: it could be partly responsible for the apparent decrease of quasi stellar objects (QSOs) at high redshifts (e.g., Fall & Pei, 1993), and for the problems encountered in detecting primeval galaxies. And dust in the precursors of present day galaxies could completely change our most recent picture of the star formation history of the universe (Madau, et al., 1996; Meurer, et al., 1997). To give just one more example, the Tully-Fisher correlation as a measure of cosmological distance depends on inclination corrections which could be significantly sensitive to dust.

Much work has been devoted to the problem of dust in galaxies, but different authors have often come up with widely differing results, from disks that appear significantly opaque to those basically transparent, and everything in between (see Davies & Burstein, 1995 and references therein). The virtues and disadvantages of different tests of opacity of spiral disks have been discussed at length in González et al. (1998, Paper I), and will not be repeated here. In summary, only one method is sensitive to the total amount of dust in disks, independently of its temperature, and does not require that any assumptions be made, neither about the geometry of the background probes, nor about the distribution of the absorbing dust (e.g., uniform screen or a collection of opaque blobs): the comparison of counts of background field galaxies seen through the disk of the foreground galaxy to counts of galaxies in reference fields.

This technique was pioneered by Shapley (1951). Before Paper I, however, it had only been applied to the Magellanic Clouds (Shapley, 1951; Wesselink, 1961; Hodge, 1974; Gurwell & Hodge, 1990). Furthermore, Wesselink (1961) deemed his results as “inaccurate”, while Gurwell & Hodge (1990) qualified theirs as “tentative.” The conceptually very simple idea of galaxy counts is, unfortunately, very difficult indeed to implement in practice. The general distribution of disk stars in the foreground galaxy reduces the contrast of background galaxy images, and foreground star clusters and H ii regions add strongly to the confusion. A first attempt to measure this effect was made by MacGillivray (1975), who used a plot of galaxy counts vs. star counts in the Small Magellanic Cloud to subtract the effects of “masking by stars” from those of extinction. In Paper I, we developed the “Synthetic Field Method” (SFM) in order to quantify and calibrate the effects of crowding and confusion, and to determine the accuracy with which statements can be made about the total opacity and reddening of an average line of sight through a foreground galaxy. We then applied it to NGC 4536, a bright spiral galaxy in the direction of Virgo, at a distance of 16.2 Mpc (Saha et al., 1996), and to NGC 3446, a small Magellanic irregular at a distance of 11.9 Mpc. In the SFM, images of a suitable reference field are added directly to the images of the foreground galaxy under study. The extinction through that foreground galaxy is then determined by comparing the numbers of “real” background galaxies, not to the number of galaxies in the original reference field image, but to the number of galaxies in the reference field that can still be identified after its addition to the foreground galaxy image. This procedure was hinted at by Wesselink (1961), and by Gurwell & Hodge (1990); they superimposed negatives of fields with and without foreground galaxy as an experiment to estimate the reduction in galaxy counts owing to crowding and confusion, but unfortunately did not correct their results accordingly.

The advent of the Hubble Space Telescope (HST) has made the method more easily applicable, thanks to its superb resolution and sensitivity, as well as to the stability of the sky, the photometry, and the point spread function (PSF). Other factors that favor the application of the method nowadays are the availability of deep and medium deep HST images that can be used as reference fields (Williams et al. 1996; Griffiths et al. 1994; Ostrander et al. 1998 Ratnatunga, Griffiths & Ostrander, 1999a, b), and the access to automated surface photometry packages (for example, Bertin & Arnouts, 1996), which allow an accurate assessment of galaxy detection limits.

De Vaucouleurs (1995) suggested that, after the Magellanic Clouds, the best targets would be the next largest galaxies in angular size, i.e., M 31, M 33, M 51, M 81, M 101. Of course, the surface density of background field galaxies is roughly constant in any direction of the sky, but closer galaxies will cover a larger angular area, allowing for results with better spatial resolution. Indeed, M 31 appears as an excellent potential candidate to apply the SFM because it is so nearby, and so well studied. Not only it could offer a unique opportunity to study the distribution of dust on small scales, and in particular to investigate the differences between arm and interarm regions, and between different galactic radii. But also, large–scale, high–resolution maps of the 21–cm line of H i (Brinks & Shane, 1984), and of the 2.7–mm line of CO (Loinard et al., 1999; Nieten et al., 2000) have been obtained, as well as large–scale maps of the distribution of warm dust (Haas et al. 1998, and references therein). Once the total dust content is known, it would become possible to determine unambiguously whether there exists a significant fraction of very cold dust (at T $`<`$ 10 K) – as has been suggested in the past – that has so far escaped detection in the far–infrared. Also, assuming a canonical dust–to–gas ratio, the total gas content could be measured by means completely independent from the usual methods. In particular, the mass of molecular gas could be estimated without having to rely on the (highly uncertain – see Loinard & Allen 1998) CO to H2 “conversion factor”.

### 1.2 It is much harder than it looks

In a nutshell, the motivation for this paper is the fact that – as we shall see momentarily – applying the SFM to HST Wide Field Planetary Camera 2 (WFPC2; Holtzman, et al., 1995) data of M 31 or the LMC has not proven easy or very fruitful. Almost no background galaxies were visible through the disk of M 31, and very few could be seen through the Magellanic Clouds (see below). In the case of the LMC, we were particularly surprised because a very simple calculation (dividing 1010 stars by the area of the galaxy or $``$ 68 kpc2) shows that there is a star only every $``$ 4 linear Wide Field Camera (WFC) pixels; i.e., with 0$`\stackrel{}{\mathrm{.}}`$1 pixels, we should already be looking between the stars!111M 31 is about 14 times further and has more stars per unit area (see §4) than the LMC, so the WFC Andromeda data should not have “empty” pixels. But Virgo galaxies are more than 20 times further away than M 31 and hence have many more stars per pixel, and yet the SFM was successfully applied there. Moreover, since it is also very difficult to see any reference galaxies in the synthetic combined images, it is impossible to tell whether the paucity of “real” background galaxies is caused by extinction or by lack of contrast against the galaxy foreground. In order to illustrate the problem, figure 1 shows color cutouts of one of the WFPC2 subfields of the LMC we have used, without (left) and with (right) the HDF-N added in; finding galaxies, especially the simulated ones, behind NGC 4536 in Virgo, in data with similar surface brightness, had been comparatively much easier – albeit easy, it was not.

In summary, while the SFM provides useful results at the distance of Virgo, it seems to fail in nearby objects. The purpose of the present paper is to determine quantitatively the conditions in which the SFM can yield useful results, and whether we can expect (or when we can expect) to use it to determine accurately the total opacity of Local Group and nearby galaxies. To this end, we will first look into the characteristics, both of the NGC 4536 images where the SFM was applied successfully, and of those images where it failed. Subsequently, we will perform simulations at different levels. First, we will artificially reproduce the stellar fields of our selected images; we will constrain these simulations with published luminosity functions, as well as with the stellar counts, surface brightnesses, and pixel–to–pixel dispersions of the data. Next, we will simulate either changes of surface brightness at fixed resolution, or changes of resolution (to mimic improvements that can be expected with increasingly bigger, diffraction–limited telescopes). To each simulated image with different brightness/resolution, we will add the HDF-N (used as reference field), in order to estimate the number of background galaxies that we can expect to see in the absence of extinction.

## 2 The observational material

All the actual data used are archival HST WFPC2 images. Table 1 lists the coordinates of the targets, and the total exposure times and observation bands. With the exception of the HDF-N, the data were retrieved from the archive after pipeline calibration, and subsequently cleaned of “warm” pixels and cosmic rays, aligned, and coadded. A full description of these steps for the NGC 4536 data has been given in Paper I; there are two small differences for the additional data used here. First, cosmic rays in the M 31–Inner field and in the LMC data were cleaned by rejecting very high counts in pixel stacks (Voit, 1997), and not by anti-coincidence in CR-split exposures (Anderson, King, & Sosin, 1995). One of the 3 images in each of the two (F555W and F814W) M 31–Inner mosaics has a much shorter exposure time (300 and 200 s, respectively, in $`V`$ and $`I`$). Regarding the LMC, there were only two images available in each filter (F606W and F814W); we made a stack of the 4 images to identify the cosmic rays. The second difference in the reduction is that, for the Local Group data, no relative shift between individual images had to be found, since all the exposures for each field were taken, respectively, on the same epoch and with the same pointing.

### 2.1 HDF-N

We used the Version 2 Release WFC images of the Hubble Deep Field (HDF-N; Williams, et al., 1996). There are 3 contiguous WFC fields in this data set (and one PC field which we do not use here); we refer to them as WF2, WF3, and WF4. Each WFC field was recorded in the 4 photometric bands designated F300W, F450W, F606W, and F814W (cf. Voit, 1997); we have not used the images in the UV band F300W. Two processing steps, described at large in Paper I, were carried out on the HDF-N data set in order to render it compatible with the other archival WFPC2 images used here. Briefly, the HDF-N images were first re-binned to their original pixel size of 0$`\stackrel{}{\mathrm{.}}`$1; second, since there is no HDF-N image in the F555W passband, we created one by interpolating linearly between the existing F450W and F606W images, taking into account the differing responses of the three filters. These “interpolated” F555W image was not needed for the LMC simulations, given that the LMC data were taken through the F606W filter.

### 2.2 NGC 4536

The images of NGC 4536 were originally acquired with the WFPC2 for the calibration of type Ia supernovae in nearby galaxies as standard candles (Saha, 1996); they have already been published by Saha et al. (1996), who give a very detailed account of the observations and the pipeline calibration (i.e., bias- and dark- subtraction, and flat-fielding) of the data.

There are thirty–four 2000 s individual exposures of NGC 4536 in the F555W passband, taken at 17 discrete epochs between 1994 June 3 and 1994 August 9; in the F814W passband, ten 2000 s images were taken at 5 discrete epochs, over the same two-month period. We cleaned the pipeline-calibrated frames of “warm pixels” and cosmic-rays, aligned them, and combined all the data in each passband into a single deep image; the resulting mosaics have total exposure times of, respectively, $`6.8\times 10^4`$ seconds in the F555W passband, and $`2\times 10^4`$ seconds in the F814W passband.

### 2.3 M 31 and the LMC

Two fields were considered in M 31. The data on the “Inner” field (M 31–Inner) were obtained on 1997 August 14 for a program on stellar populations in the M 31 disk (Trauger, Holtzman, & Gallagher 1997). The total exposure times are 2400 s at F555W and 2500 s at F814W. The field is located at about 5 kpc from the center of M 31. It was chosen specifically to avoid spiral arms and to have low, uniform extinction; hence, it shows very little structure. The data on the “D478” field of M 31 (M 31–D478) were obtained in parallel mode on 1996 July 8, as part of the Medium Deep Survey (Griffiths, et al., 1994; Ostrander et al., 1998; Ratnatuga et al., 1999a, b). The total exposure times are 5600 s at F555W and 6400 s at F814W. We retrieved them because the field overlaps part of the dark dust cloud D478, which is located at 2.5 kpc from the center of the galaxy (Hodge, 1981), and has been extensively studied by Allen and co–workers (Allen & Lequeux 1993, 1994; Loinard & Allen, 1998).

The data on the LMC were obtained in parallel mode on 1996 January 30, also as part of the Medium Deep Survey. The field is located $``$ 1$`\stackrel{}{\mathrm{.}}`$3 from the center of the LMC, and these same observations have been amply discussed by Elson, Gilmore, & Santiago (1997).

## 3 Diagnostic of the problem

Before even trying to test for the presence of extinction in the M 31 and LMC fields, we added to each of the WFC subimages each of the 3 HDF-N WFC chips. These “baseline” simulations, where the 3 chips of the control field are added one at a time without attenuation, establish the capability of the method to discriminate between crowding and extinction (Paper I).

In the case of the M 31–D478 field, regardless of what WFC subimage we analyzed, we could barely see the one or two brightest galaxies of each of the HDF-N chips. We could also see at most one or two “real” background galaxies in the whole WFC field. In the M 31–Inner field, we could see an average of 7 HDF-N galaxies and 4 “real” galaxies per WFC chip. As for the LMC, the simulations yielded an average of 22 galaxies per chip, and we could count 10 “real” background galaxies per subimage, on average. Assuming that the error in the galaxy counts equals $`\pm `$3.5$`\times N^{1/2}`$ (see §5), these numbers translate roughly into $`A_V=1.2\pm 3.1`$ mag for the M 31-D478 field; $`A_V=0.3\pm 1.3`$ mag for the M 31-Inner images; and $`A_V=0.9\pm 1.1`$ mag for the LMC. Given the errors, these are not tremendously useful results! For comparison, in similar “baseline” simulations of NGC 4536 (Paper I), we found an average of 28 HDF-N galaxies in each of the 2 subimages that comprised the arm region (WFC2 + WFC3), and of 38 control galaxies in the interarm (WF4) region; the corresponding uncertainties on the opacity were 0.2–0.3 mag.

We hypothesized early that image “granularity” might have an important role in the success or failure of the SFM; hence, we searched for diagnostics of such image granularity. Average means or standard deviations calculated across entire images would not work, because “smooth” data with large-scale spatial gradients and “flat” data with large pixel-to-pixel differences would produce similar values (see Figure 2). Consequently, we considered instead the distributions of pixel means and rms variations in squares 100 pixels on the side. We compared the different distributions shown by the M 31, LMC, and NGC 4536 data; they are listed in Table 2. While the means reflect the surface brightnesses, the standard deviations constitute a diagnostic of the image granularity. For example, very broad distributions of mean values point to large-scale structure, due, for instance, to stellar associations and dust lanes.

The case of the M 31–D478 field is clearcut. It is on average more than two times brighter than NGC 4536 (Table 2). More significantly, it is much more “granular”, as shown by a mean rms 5 times larger, relative to the mean surface brightness. Large scale structure does not seem to contribute to the problem: the FWHM of the distribution of surface brightnesses in D478 is $``$ 14%, compared to the mean; this is of the same order as the spread in the interarm field of NGC 4536 (12%), the subimage with the least large-scale structure. The M 31–Inner field is indeed quite structureless, as attested by a distribution of surface brightnesses with an average FWHM of less than 1%! It is barely 10% brighter than NGC 4536 but, on the other hand, it is again 5 times more “granular.”

Finally, the LMC field is hardly brighter than NGC 4536; with the exception of the chip WF4, the structure shown by the data is also of the same order. However, the pixel-to-pixel variations are larger than 300% of the mean. Such a large value of the rms is, indeed, a sign that we are looking between the stars.

## 4 Simulations of stellar foreground

To quantify the effect of the granularity, we mimicked the various observations used with artificial stellar fields, and measured their effect on the detectability of background galaxies. These simulations were all produced in the $`V`$ passbands, i.e., as if they had been obtained through the HST filters F555W or F606W. In general, colors are crucial for the identification of background galaxies and hence for the success of the SMF. However, the present work is centered on simulations of stellar populations; population models, as well as magnitudes for different stellar types, are most readily available at $`V`$.

The simulations had several levels: first, we reproduced the LMC and M 31–Inner HST WFPC2 data, guided by published luminosity functions and/or the star counts of the data itself, and constrained by the mean and rms of the images. We chose to use the M 31–Inner field because of its ,“intermediate” surface brightness and because, as stated above (§2), it has very little large-scale structure and is dominated by image granularity. It has the additional advantage that it has roughly the same surface brightness as the NGC 4536 data (Table 1), so we could also use it as a basis to simulate a galaxy disk at the distance of Virgo (see below §6. In the case of the LMC, we really wanted to know why we could not detect more background galaxies in between the stars. At this stage, we always used a Tiny Tim (Krist & Hook, 2001) PSF.

The second step was to simulate the change in resolution offered, first by the smaller pixels of the Advanced Camera for Surveys (ACS), which has just been put on-board the HST, then by increasingly bigger, diffraction–limited telescopes. This was achieved by changing the number of stars per pixel. Here, we did two subclasses of simulations. In one case, we used Gaussian PSFs, and assumed that they were fully sampled. In the other, we used the Tiny Tim PSFs foreseen for the two cameras (Wide Field, ACS-WFC, and High Resolution, ACS-HR) of the ACS. Since the PSF for the ACS-HR should be Nyquist-sampled in the visible, we used its PSF to simulate the observations with bigger telescopes.

For completeness, we have also done a slightly different kind of simulation, which will have immediate applicability to predict the potential success of the SFM in existing archival images. We have mimicked changes in surface brightness with radius of the disk of M 31, at the resolution and with the PSF of the WFPC2 on board HST. Assuming a constant stellar population, this also boils down to changing the number of stars per pixel. Finally, for every object, and at each different brightness/resolution, we have added the HDF-N to the data, in order to assess the number of background galaxies that we can expect to see in the absence of extinction.

The background galaxies have been counted by eye. We have performed multiple checks since we started working on the subject and, as long as it is the same person who identifies the “real” and the control galaxies, this procedure yields robust results. Eventually, it would be desirable to use an automated photometry package especially adapted to this kind of study; however, automatic star/galaxy separation in crowded fields is quite difficult. For example, we have tested SExtractor (Bertin & Arnouts, 1996) on the F606W the WF3 field of the LMC (§2.3), both with and without the HDF-N added in (see Fig. 1). The program finds significantly more extended objects than we do galaxies by eye, but it is unclear which of the identified structures are actual galaxies, since SExtractor finds more (185 vs. 153) of these extended objects in the image without the HDF-N! In the case of our simulations, which do not include dust, HII regions, etc., star/galaxy separation should be in principle easier than in real life; but in any case we considered this effort as unnecessary and beyond the scope of the present work. Lastly, we do not investigate the application of the SFM to hypothetical images produced by a large IR telescope in space, like the New Generation Space Telescope (NGST) is anticipated to be. This would have necessitated simulations not only of the stellar foreground, but also of the background galaxies.

We present results for infinite exposure times (noiseless data), and also for “realistic” integration times, i.e., 1800 s or 7200 s with the HST. For telescopes with a hypothetical aperture D in meters, these times should be scaled by (2.4/D)2. The area of sky considered is always the same for all examples in the paper, and consequently the clustering error (§5) has always been assigned in the same way. 222We assume that the reduced field–of–view of bigger telescopes (which, in the absence of special optics/detectors, is also proportional to D-2) will be compensated by the smaller exposure time required to achieve the same depth, and that it will be possible to make mosaics with the same area in about the same total exposure time.

### 4.1 M 31–Inner

We used stellar population models, and the constraints offered by the data, to determine an adequate mixture of stars. The bulk of them were taken to belong to an old population with a main-sequence (MS) upper-cutoff type F5V and 0.05% by number of K giants (McClure & van den Bergh, 1968). To account for the fact that we were trying to simulate a disk, while McClure & van den Bergh (1968) studied the bulge of the galaxy, we then added 2 % by mass of a younger population (for this paper, percentage by mass is equivalent to percentage by number, because upper mass limits of different MS components are not so discrepant, and massive stars contribute very little to the total stellar mass). We have adopted this number for the contribution by young stars following our fit to the arm and interarm region populations of M 99 (González & Graham, 1996).333 There is also a hand–waving argument: if a galaxy has lived about 10 Gyr, and the spiral pattern goes around every 108 years, and the passage of the arm induces star formation (which is of course debatable and debated; e.g., Elmegreen, 1987), the youngest population would be of the order of 1 %. In these and all subsequent cases, main-sequences go down to M6 stars, or about the hydrogen-burning limit. At this point, the two free parameters left were the upper-cutoff stellar type of the young population, and the number of stars per pixel of the old population. The combination that best fit the mean and rms of the data was an upper cutoff type B5V, and 50 old stars pixel-1; the red giants and the young stars scaled accordingly. Stellar population synthesis models (Charlot & Bruzual, 1991; Bruzual & Charlot, 1993) with the described composition will have a mass-to-light (M/L) ratio of 1.8 which, assuming a Salpeter initial mass function (IMF), will yield about 100 stars per pixel. However, the factor of two discrepancy could very well be explained by a different percentage of young population, a different IMF, a different foreground Milky Way extinction towards M 31 (which is hard to measure precisely because of the dust and gas in the Andromeda galaxy itself), possible error in the distance to Andromeda, and by internal extinction in M 31 (which, of course, is what we are trying to measure). We also note here that including more luminous stars (B2V, say) in the young population makes it possible to match at the same time the surface brightness of the data and the number of stars per pixel required by the M/L ratio of the stellar population synthesis models; however, it is questionable that we see B2V stars in the data, and the inclusion of these brighter stars increased the standard deviation of the simulations much beyond that of the data.

Salpeter luminosity functions were used for the main-sequences, and a power-law of the magnitude with an exponent of 0.7 for the red giants (Iben, 1968). A uniform spatial density was used for all three components. The lower and upper absolute magnitudes for all components of the stellar population are given in Table 3. We also show the number of stars per pixel, in each component, for our simulation of the M 31–Inner field at the resolution of the WFPC2; the relative contributions of the three components remain fixed at all angular resolutions.

Figure 3 shows the V vs. V-I color–magnitude diagram (CMD) for the stars in the real M 31–Inner data. Stellar cores were fitted with DoPhot (Schechter, Mateo, & Saha, 1993), zero-points were found iteratively as prescribed by Holtzman, et al. (1995), and aperture corrections were applied (-0.62 at V, -0.77 at I). No reddening correction was performed. The box delineates the approximate area of the CMD populated by B stars for a distance modulus to M 31 of (M-m)o = 24.47 (Stanek & Garnavich, 1998); the rest of the stars in the diagram are either red giants or asymptotic giant branch (AGB) stars. Assuming a Galactic reddening towards M 31 of E(B-V) = 0.062 (Schlegel, Finkbeiner, & Davis, 1998); an internal reddening of E(B-V) = 0.05, or a third (Charlot & Fall, 2000) of that found by Massey et al. (1995) for OB associations in Andromeda; and R $`A_V`$ / E(B-V) = 3.1, B5 MS stars ($`M_V`$ = -1.1) will have $`m_V`$ = 23.7 in this diagram. (B2V stars, $`M_V`$ = -3, will have $`m_V`$ = 21.8.) The incompleteness in the diagram prevents us from using it to improve the estimate of the percentage of young population.

The positions and magnitudes for the artificial stars needed for the simulation were generated with the task starlist, within the software package IRAF444IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. (Tody, 1993). Computer memory limited the number of stars that could be generated in each run of starlist to about 5x105. This is sufficient to produce a square of 100 pixels on the side for the simulation of the M 31–Inner field with the resolution of the WFPC2.

Next, the task mkobjects, also within IRAF, was used to produce an actual image of the simulated stellar field out of the list of stars given by starlist. With Tiny Tim (Krist & Hook, 2001), we generated a PSF appropriate for a G5V star in the F555W passband, at the center of chip WF3; we used only one PSF for all stars, regardless of their stellar types or their positions in the image. We produced only one simulated 100x100 square to reproduce the M 31–Inner data; in order to create a frame with the angular size of the WFC subimages (80 ″on the side), we copied the artificial stellar field many times over. This means that the simulation with the resolution of the WFPC2 reproduces only the means, and not the spreads, of the distributions of surface brightnesses and standard deviations of the 100x100 squares into which the original data were divided (see §3); the simulation does not mimic the large scale structure of the data that manifests itself in such spreads. Table 3 lists the mean and pixel-to-pixel rms of the “mother” M 31–Inner simulation; a theoretical “sky background” produced by the WFPC2 Exposure Time Calculator (ETC; Biretta et al., 2001) is included for easy comparison with the data.

When changing resolution, two sub–types of simulations were implemented. In one case, we employed fully sampled Gaussians; in the other, we used Tiny Tim (Krist & Hook, 2001) to produce a PSF, again appropriate for a G5V star in the F555W passband, again the same for all stars. For the second sub–type, depending on the simulated resolution, we used either the model PSF for the center of the ACS-WFC, or that for the center of the ACS-HR camera. We produced once more only one square per angular resolution step. For simulations with Gaussian PSFs, we always produced 100x100 squares; for simulations with realistic PSFs, the sides of the squares scaled as the inverse of the pixel sizes. In all cases, we used the “tiling” technique to make frames 80 ″on the side. The final step of the simulations was to add the HDF-N chips and count how many galaxies could be seen, on average, in these simulated images the size of the WFC field of the WFPC2 (5.3 arcmin2).

### 4.2 LMC

As a basis for the LMC simulations, we used the luminosity function published by Jahreiss & Wielen (1997) for nearby stars, but corrected as in Holtzman et al. (1997) for the LMC; that is, numbers of stars with $`M_V`$ = 7 or fainter have been multiplied by 1.6. We assumed a distance modulus to the LMC of (M-m)o = 18.6 (Groenewegen & Oudmaijer, 2000), and an extinction correction of AV = 0.63 (Elson et al., 1997). The luminosity function $`\mathrm{\Phi }(M_V)`$ is presented in Table 4, in the same units used by Jahreiss & Wielen (1997), or stars within 20 pc in 1 mag bin.

In the case of these simulations, we did not try to match the characteristics of the original data, but those of WFPC2, chip 3, after editing out the saturated stars. These are labeled WF3-E in Table 2. There are three other differences with respect to the M 31–Inner simulations. First, we used a Tiny Tim PSF appropriate for an F6V star through the F606W filter for the simulations with a “realistic” PSFs. Second, the number of total stars in a square 80″ on the side is such ($`<`$ 2x104) that we had enough memory to generate a field in one pass of starlist, as opposed to having to cover it with many copies of the same “tile.” This led to the third difference, i.e., that we attempted to match with the “mother” simulation not only the mean surface brightness and the mean standard deviation of the data, but also the spread of surface brightnesses. The “mother” LMC simulation has 0.027 stars pixel-1 (or $``$ one star every 6 linear pixels, very close to our back of the envelope calculation in §1.2), $`\mu _V`$ = 21.68, FWHM = 0.29 mag (26% of the mean), and $`\sigma `$ = 4.79 mag (441% of the mean555This is 20 % too high compared to the edited data, but we could not improve it and keep the chosen luminosity function at the same time.).

### 4.3 Virgo

Here, we simulated four different situations based on the NGC 4536 data: an interarm region with the inclination to the line of sight of NGC 4536 (65°); an arm, also with such inclination, a face–on interarm region, and an arm in a face–on galaxy. For the interarm regions, we used the same stellar mixture as for the disk of M 31; for the face–on simulations, we “deprojected” the $`i`$=65° populations as if they were optically thin. Interestingly, given that the inclination of NGC 4536 is smaller than that of M 31 (78°), we needed the equivalent of 60 stars (instead of 50 or fewer) at the distance of Andromeda to match the characteristics of the NGC 4536 data. We assumed an apparent distance modulus to Virgo of (m-M)$`_{A_V}`$ = 31.23 (Saha et al., 1996); given the apparent distance modulus to M 31 of (m-M)$`_{A_V}`$ = 24.82 (see §4.1), we were trying to simulate a very similar stellar population located $``$ 19 times further away. Naturally, computer memory limited us to 5x5 pixel “tiles” for the “mother” simulation of the “inclined” Virgo galaxy; however, these were so small that, at each resolution step, we produced 400 of these “sub–tiles” to generate a 100x100 pixel “tile.” We used the same size of “tile” for the face–on cases, even though (assuming the same population) these had 2.5 times fewer stars per unit area \[cos(65°) $``$ 0.42\].

In order to mimic “arms,” we kept the old MS and the red giants in the same numbers as for the “interarm” simulation, but this time we added a young main sequence with an upper type B3V (10 M), and $``$ 2.4 % by mass of the old main sequence. For stars with masses between 9 M and 10 M, this is approximately equivalent to a density of 0.01 stars pc-2. Hunter (1995, and references therein) have found densities between 0.007 and 0.02 massive stars pc-2 for OB associations and selected regions of the Milky Way, the LMC, I Zw18 and M 31. Interestingly, we did not need to increase the numbers of the old population stars, as one would in principle have to do in the presence of a density wave; pushing up the fraction of old population augmented the mean surface brightness while decreasing the standard deviation under that shown by the data. As for M 31–Inner, the parameters and characteristics of the “mother” simulations are given in Table 3.

### 4.4 M 31 radial simulations

As we mentioned earlier, we have simulated the changes in surface brightness with radius of the disk of M 31, at the resolution and with the PSF of the WFPC2 on board HST. We have used the stellar mixture that worked for the M 31–Inner field, and changed the number of stars per pixel to reproduce the different surface brightnesses found between 10 and 100 galactocentric radius, at 10 steps. The surface brightness profile has been taken from Walterbos & Kennicutt (1988). The radial run of surface brightnesses numbers of old MS stars per pixel are given in Table 3.

## 5 Results

Figures 47 show the results of our simulations; figures 46 present information in basically the same fashion. Background galaxy counts observed are plotted vs. number of stars per pixel. The field of view of 5.3 arcmin2 is that of the combined 3 WFC subimages of the WFPC2. The left panels display the number–counts obtained with “realistic” Tiny Tim PSFs, while the right panels show results for Gaussian, fully–sampled PSFs (except for the leftmost square points, which show once more results for WFPC2 aboard HST). Empty symbols are for “noiseless” (infinite integration time) data, whereas full symbols are for 7200 s (figures 4 and 6, for M 31–Inner and Virgo, respectively) or 1800 s (Figure 5, LMC) exposures with HST; these times should be scaled by (2.4/D)2 for telescopes with a different aperture D in meters. It is indicated when results correspond to diffraction limited telescopes. For the “noisy” data, we have used the gain and read–out noise of the WFPC2 camera, except for the ACS simulations in the left panels, where values anticipated for the ACS were applied (Pavlovsky et al., 2001); the nominal, foreseen gain has been multiplied by 3, given that the ACS is expected to be 3 times faster than the WFPC2. In all cases, sky backgrounds were taken from the WFPC2 ETC (Biretta et al., 2001). For reference, in the absence of any foreground object and using –only in this case– the photometry package SExtractor (Bertin & Arnouts, 1996), we have identified $``$ 1700 galaxies, with an isophotal $`V`$ detection limit of 27 magAB arcsec-2, in the 3 WFC chips of the virtually noiseless HDF-N.

In the past (Paper I; González, Fruchter, & Dirsch, 1999), we have derived total (Poisson plus clustering) errors for the galaxy counts, as a function of total magnitude completeness limit, from values of the two–point correlation function, $`\omega (\theta )`$, published in the literature (Roche et al. 1993; Brainerd, Smail, & Mould, 1995). These theoretical total errors amount to $``$ 2$`\times \sigma _{Poisson}`$; $`\sigma _{Poisson}=N^{1/2}`$,where $`N`$ stands for the galaxy counts. However, in González et al. (1999) we found empirical errors reaching close to twice this value. Interestingly, Roche et al. (1993, columns 3 and 4 of their Table 2) obtain similar empirical errors, albeit with significantly larger fields ($``$ 23.5 arcmin2). Their predicted error ranges between 1.5–2$`\times \sigma _{Poisson}`$, while their empirical error varies between 2.5–3.5$`\times \sigma _{Poisson}`$. Tracking down the reason for this discrepancy is beyond the goals of the present work but, in the meantime, we will adopt $`\pm `$3.5$`\times \sigma _{Poisson}`$ error bars. This value should be conservative enough, as well as convenient, since in all the simulations we have identified background galaxies by eye (§3) and hence only could have an approximate estimate of the detection completeness limits.

Simulations for the M 31-Inner field with realistic PSFs (Figure 4a) and exposure times show that, while the ACS will multiply the number of background galaxies that can be detected by two compared to WFPC2, then there will be no significant further improvement until we can count with a diffraction–limited 10–m optical telescope. The galaxy counts in the simulations with the Gaussian PSFs (Figure 4b), however, improve steadily with the size of the telescope (although here we invite the reader to notice the jump in the counts around 1 star per pixel); already a diffraction limited 3–m telescope would be almost as good as the 10–m telescope in the left panel, or 5 times better than WFPC2. The 4.5–m telescope with a Gaussian PSF would be 10 times better than today’s WFPC2. With our extremely conservative error bars, we would expect the uncertainty on the extinction determination to decrease from $`\pm `$ 1.2 mag with WFPC2 to $`\pm `$ 0.9, $`\pm `$ 0.5, and $`\pm `$ 0.4 mag, respectively, with 2, 5, and 10 times more galaxies.

The LMC represents the most dramatic case. Figure 5a shows basically constant number–counts; the ACS will already be close to the meager factor of two improvement we expect from a diffraction–limited in the optical, 10–m telescope. Figure 5b again shows a steady increase in the number–counts with resolution, with a factor of 8 more counts already expected for a 3–m diffraction limited telescope! It also illustrates the situation where one benefits the most from long exposure times. The values above translate into uncertainties in the extinction of $`\pm `$ 0.8 mag now, $`\pm `$ 0.5 mag with the ACS, and $`\pm `$ 0.3 mag with 8 times more background galaxies.

Figure 6 displays results for a Virgo disk. The four top panels (a–d) summarize simulations for a galaxy with an inclination to the line of sight of 65°; the bottom panels (e–h) correspond to a face–on disk. Finally, panels (a–b) and (e–f) are for “interarm” regions, while panels (c–d) and (g–h) represent “arms.” The most interesting lesson from this figure is that, when using “realistic” exposure times, all panels are basically flat, with improvements that are, at most, of a factor of two. Also, for each different case (from top to bottom: interarm inclined, arm inclined, interarm face–on, arm face–on), the results are the same regardless of the PSF used. Thus, in the case of Virgo we are already close to the best results we will be able to get with the present day HST. The extinction uncertainties range in most cases between $`\pm `$ 0.5 and $`\pm `$ 0.3 mag, in perfect agreement with Paper I.

Finally, Figure 7 shows background galaxy counts that should be observable at different radii through the disk of M 31 with WFPC2 data, assuming everywhere the same population mixture used to simulate the M 31—Inner field. Numbers of stars per pixel and $`\mu _V`$ from Walterbos & Kennicutt (1988) are also shown. All the simulations in this series employ a “realistic” Tiny Tim PSF; “noisy” ones (empty symbols) have exposure times of 1800 s. With existing data, it will be hard to reduce the extinction uncertainty below $`\pm `$ 0.5 mag, even with long exposure times and for the largest radii investigated.

Table 5 lists, for each simulation, the pixel size, the linear size of each “tile,” the number of stars per pixel, the PSF used, the background added, exposure time, gain, read–out noise, and galaxy counts in a field of 5.3 arcmin2. For comparison, the number–counts in a simulation where the HDF has been added to the real data are also shown.

## 6 Discussion

As hypothesized initially, the granularity of the foreground galaxy (as measured by the number of stars per pixel) does appear to play a key role in the detectability of background galaxies. However, the comparison between diffraction–limited telescopes with realistic and Gaussian PSFs (Figures 46) shows that the “quality” of the PSF often plays a significant role as well. This points to two distinct sources of granularity in the images. A “primary” granularity due to the stellar distribution itself – and best measured in the simulated images with Gaussian PSFs – and a “secondary” granularity produced by the wings and spikes of the PSF. At the distance of Virgo, the number of stars per pixel is always very large or large (between 50 and 104, depending on resolution), and each individual star is very faint, so the stellar content of the foreground galaxy appears as a smooth continuum with very little sub–structure (i.e. almost without primary granularity). Changing the quality of the PSF (from a perfect Gaussian to a realistic case) has almost no influence, because the additional “secondary” granularity introduced by an imperfect PSF has essentially as little sub–structure as the stellar field itself.

The situation in the LMC is almost exactly opposite to that of Virgo. In the LMC, there is only one star every several tens or hundreds of pixels (again, depending on resolution); in addition, owing to the proximity of the LMC, this star is often bright. Although this implies a “primary” granularity significantly larger than in Virgo, one would a priori expect to be able to see “between” the stars. Only a few pixels out of several tens or hundreds would image (and, therefore, be affected by) stars proper in the LMC, and most pixels would be expected to be free of such emission. That so few background galaxies can be identified in the LMC implies that, there, the secondary granularity dominates completely over the primary one. A quite natural explanation to this phenomenon is that a realistic (non-Gaussian) PSF will spread a significant fraction of the stellar light over many pixels, in effect “propagating” the granularity to pixels that would otherwise be emission–free. For a 10-m diffraction–limited telescope, up to 90% of the background galaxies can be lost to this secondary granularity (Figure 5). Figure 8 demonstrates this dramatic effect for the case of the ACS-HR camera (i.e., a diffraction limited 2.4-m telescope). This should be a good general example because the PSF of this instrument is well sampled (which is also the reason why we have used the ACS-HR PSF to simulate images obtained with larger telescopes). For reference, the top panel of this figure shows WF3 subfield of the WFPC2 F606W data of the LMC. The middle left panel presents our simulation of the same field with the resolution of the ACS-HR and the Tiny Tim PSF for said camera and wavelength (we remind the reader that we have omitted saturated stars from the simulation). The middle right panel shows our simulated field with a realistic PSF plus the HDF-N; like in the case of the right panel of figure 1, reference galaxies can only barely be seen. The bottom panels display, respectively, the simulation with a Gaussian PSF without (left) and with (right) the HDF-N added in; tens of background reference galaxies can easily be counted in the right bottom panel. Figure 9 shows linear (a and c) and logarithmic (b and d) surface plots of the ACS-HR (top) and Gaussian (bottom) PSFs. As can be seen in the linear plots, the cores are similarly sampled, but the core of the Gaussian out to 3$`\sigma `$ receives 97% of the total light, whereas only 58% enters in the same area of the ACS-HR PSF. The wings of the PSFs are better appreciated in the logarithmic plots; the light in the wings and spikes of the realistic PSF reduces the signal–to–noise ratio of background galaxies and explains the enormous disparity in the efficiency with which they are identified in both cases. A bad subtraction of the realistic wings can only introduce more noise and further reduce the detectability of background objects.

The case of M 31 is intermediate between Virgo and the LMC. There, changing the resolution of the telescope helps little, until the number of stars per pixel falls below about 1 where a jump can clearly be seen on Figure 4. This jump corresponds to a change of regime: with less than 1 star per pixel, one encounters a situation similar to that of the LMC, where the secondary granularity becomes strongly dominant. Indeed, the difference between Gaussian and realistic PSFs becomes important only in that regime. With a few bright stars per pixel, the primary granularity is the dominant factor. This is also quite natural, because this is precisely the case where all pixels image stars, and at the same time large pixel-to-pixel variations (primary granularity) are expected. This situation, where each pixel contains more than 1 star, and where stars cannot be identified independently or their positions within a pixel determined accurately, essentially corresponds to the regime of surface brightness fluctuations (SBFs; Tonry & Schneider, 1988).

The comparison between M 31 and Virgo shows that SBFs largely prevent the detection of background objects only in Local Group galaxies, where each pixel contains a small number of bright stars. At the same time, the simulations with realistic PSFs show that there is always at least a slight improvement in the detectability of background galaxies with better resolution (which appears paradoxical, since fewer stars per pixel seem to imply more granularity). This improvement, moreover, is significantly smaller in Virgo than in the LMC. All these facts can be understood in a statistical way. For simplicity, we will assume that all stars are identical, and have an intrinsic flux $`f_{}`$ (in photons sec-1 received at the detector). The total number of photons $`\gamma `$ received in each pixel of the CCD in an integration time $`t`$, assuming that there are $`n`$ stars per pixel is:

$$\gamma =ntf_{}$$
(1)
The integration time $`t`$ can be considered a constant here, but both $`f_{}`$ and $`n`$ contribute to the noise. Therefore, neglecting the sky background and the read-out noise (RON) of the detector, the total noise is:

$$\sigma =\sqrt{\left(\sigma _nf_{}t\right)^2+\left(n\sigma _f_{}t\right)^2}$$
(2)
Since the uncertainty on $`f_{}`$ is $`\sqrt{f_{}}`$, and that of $`n`$ is $`\sqrt{n}`$, we obtain the following expression for the total noise:

$$\sigma =\sqrt{nf_{}^2t^2+n^2f_{}t^2}=t\sqrt{nf_{}^2+n^2f_{}}$$
(3)
A background galaxy will be detected with a signal-to-noise ratio equal to the quotient of the number of photons received from that background object ($`f_{bg}t`$) divided by the noise computed above (the Poisson noise due to the background galaxy can be neglected, since the surface brightness of such objects is always so much fainter than the foreground stars):

$$\left(\frac{S}{N}\right)_o=\frac{f_{bg}t}{t\sqrt{nf_{}^2+n^2f_{}}}=\frac{f_{bg}}{\sqrt{nf_{}^2+n^2f_{}}}$$
(4)
When the pixels are made $`L`$ times smaller (on the side), the flux received from the background galaxy ($`f_{bg}`$) and the number of stars per pixel $`n`$ both diminish by $`L^2`$ (the background galaxy is assumed to be resolved – if it was not, it could not be identified as such in our simulations), but the flux of each star stays constant (given that stars are unresolved and always fall completely in one pixel). The signal-to-noise ratio per pixel, therefore, becomes:

$$\frac{S}{N}=\frac{f_{bg}}{\sqrt{L^2nf_{}^2+n^2f_{}}}$$
(5)
But the pixels are now $`L`$ times smaller on the side, so the background galaxy now occupies $`L^2`$ more pixels; and since the noise adds in quadrature, the background galaxy is now detected with a total signal–to–noise ratio:

$$\left(\frac{S}{N}\right)_L=\frac{Lf_{bg}}{\sqrt{L^2nf_{}^2+n^2f_{}}}$$
(6)
Finally, if the same foreground stellar population is pushed $`d`$ times farther away, the number of stars per pixel $`n`$ is multiplied by $`d^2`$, while the flux of each star $`f_{}`$ is divided by $`d^2`$. The flux of the background galaxy obviously remains unchanged, and the signal-to-noise ratio becomes:

$$\left(\frac{S}{N}\right)_d=\frac{Lf_{bg}}{\sqrt{\frac{L^2nf_{}^2}{d^2}+n^2f_{}d^2}}$$
(7)
The behavior of this expression as a function of pixel size (set by $`L`$) for the same stellar population at three different distances $`d`$ is shown on the left panel of Figure 10. The value of the various parameters have been chosen to reproduce (although with arbitrary units) the conditions found in the LMC, M 31 and Virgo. The exact normalization factors will depend on the way signal-to-noise and background galaxy surface brightness detection limits are related, an issue that we do not investigate here. Figure 10b shows again the number of galaxies that are detected in the LMC, M 31, and Virgo as a function of the pixel size (taken from Figures 46). Those numbers are the results of the simulations with infinite integration time, and Gaussian PSFs – the best comparison data, since integration time, sky background, RON, and PSF “imperfections” are not accounted for in Equations 1–7. Also, an arbitrary constant has been added to the logarithm of the number of background galaxies seen through each object, to facilitate the comparison of the data with the model. The qualitative agreement between the left and right panels of Figure 10 is remarkable, especially given the extremely simple model we have adopted for the noise; in addition to illustrating the origin of the detection efficiency of background galaxies, the concordance between model and observations validates the approach used to count the galaxies. The number of background galaxies detected behind the LMC increases rapidly as the pixels are made smaller, until it reaches a plateau. This plateau corresponds to the situation where the signal–to–noise ratio of the galaxy detection is dominated by “the space between the stars,” or the smooth sky background in real data; the number of background objects seen is not improved with better resolution, since all the galaxies detectable with the chosen integration time are already identified. In M 31, the detectability also increases rapidly (although less rapidly than in the LMC), but never reaches a plateau (note that it would, eventually, if the pixels were made even smaller). In Virgo, the improvement in detectability as pixels are made smaller still occurs, but is much less pronounced than in the LMC or M 31, precisely because the foreground is approximately smooth. The two main features of the detectability of background galaxies in the simulations (increase followed by leveling–off in the case of the LMC; and smaller slope of the improvement for farther objects) are, thus, reproduced by Equation 7. Our approach of considering separately the contributions of $`n`$, the number of stars, and $`f_{}`$, the flux of each star, to the noise, also allows us to explain away the paradox stated at the beginning of the section. For an individual foreground object, the increasing number of stars per pixel as resolution worsens translates into smaller pixel–to–pixel fluctuations, which ought to lead to an increase in the counts of background galaxies. But of course (even though the surface brightnesses of both foreground and background galaxies stay constant), the total brightness of each pixel is also increasing, which means an increase in photon noise per pixel and a reduction in the background galaxy contrast, leading to a reduction in the counts of identified background galaxies. Of the two competing effects, the simulations show that the latter “wins.”

In the case of Virgo, we are, and will remain for the foreseeable future, within the regime of faint SBFs. Better resolution will yield only small improvement, and results will be independent of the quality of the PSF (Fig. 6 and Fig. 10). Conversely, in the LMC, we are already outside of the (bright) SBF regime with $`HST`$ WFPC2 images; in theory, improvements in resolution should translate into enormous gains in galaxy detection (Fig. 5b and Fig. 10), but Fig. 5a shows that these will not be obtained with a realistic PSF. The successful application of the SFM in this case will require an adequate subtraction of the foreground stars with their diffraction pattern, a notoriously difficult task. In the case of WFPC2 data, the difficulty arises from the undersampling, chromaticity, and spatial and time variability of the PSF. Model PSFs, such as those produced by Tiny Tim, are clearly not sufficiently accurate, especially in the wings. Better, possibly iterative, procedures will have to be developed (Anderson & King, 2000, and references therein). An attempt to tackle WFPC2 images of the LMC in this way will be presented in a subsequent paper (González et al., 2003, in prep.). It should be emphasized, however, that much work will probably be needed in order to fully master optical/near-IR PSF subtraction procedures, both for WFPC2 and for future telescopes and instruments. For ACS data, for instance, aside from the undersampling of the PSF by the WFC, the biggest problem we can anticipate is the huge geometric distortion of the field–of–view.

As for M 31, Fig. 4 and Fig. 10 show that small improvements can still be obtained with better resolution, regardless of PSF quality. Furthermore, images obtained with a 10-meter diffraction-limited optical telescope (which could be available in the not–so–distant future) will have less than 1 star per pixel. Except for the apparent brightness of the stars, such images of M 31 will be similar to the existing $`HST`$ images of the LMC, where the detectability of background galaxies is limited by the secondary granularity produced by structure in the PSF. Accurate removal of individual stars will be required at that point to improve the results of the SFM. For existing WFPC2 data, one could propose to remove at least the brightest stars. It would be very difficult, if not impossible, to do a good job, since this would require the accurate determination of both the position of the stars and the shape of the PSF in a very crowded field, with many more than 1 star per pixel. To investigate whether it would be worthwhile to pursue this idea, we generated an artificial M 31–Inner image with a list of stars from which we deleted the red giants and young MS stars brighter than magnitude 25. The identified background galaxies in a simulated 5.3 arcmin2 field went up from 18$`\pm `$15 (Table 6) to 26$`\pm `$18; the uncertainty in the extinction would go down from 1.3 to 1.1 mag. One could not remove fainter stars, since incompleteness sets in at 25 mag (Fig. 3).

Two final comments should be made here. First, integration times are never an issue for realistic PSFs. The results of the simulations for “short” (0.5–2 hours) integration times are always nearly as good as those with infinite integrations. Second, a note of caution. We have seen here that, using HST WFPC2 data, background galaxies can be detected through the disks of Virgo galaxies, but not through Local Group galaxies. It could, therefore, be proposed to simply smooth the data on Local Group galaxies to match the resolution achieved by HST in Virgo. This – of course – would be foolish, because it would also smooth the background galaxies. What one needs is a smooth foreground stellar field, but sharply defined background objects.

## 7 Conclusions

The analysis of the simulations presented in this paper shows that, for the time being, Virgo seems to offer the best prospects for the application of the SFM. Reasonably accurate statements can be made about the opacity of the disks at that distance (Paper I), and long integration times are not required. Arcsecond resolution H i VLA observations could be obtained for a significant number of the spirals in Virgo, and the Atacama Large Millimeter Array (ALMA) will be able to provide high–sensitivity CO maps with the same resolution. Obtaining the opacity of a large sample of spiral in Virgo using the SFM could, therefore, be used to statistically investigate the properties of the ISM there. Objects more distant than Virgo would tend to present even smoother foreground stellar fields, but the opacity there would be determined with less accuracy because their smaller angular size would reduce the number of background galaxies that could be identified.

For Local Group galaxies like M 31 and the LMC, the SFM will hardly be usable, even with telescopes significantly better than those currently at hand, unless the PSF is better understood. An adequate subtraction of individual point sources from WFPC2 data is an extremely difficult task, owing both to the undersampling (Anderson & King, 2000), and to the variability of the PSF with time, position in the chip, and spectral energy distribution of the imaged objects; the variability affects especially our present ability to model the wings (Krist & Hook, 2001), i.e., the source of the secondary granularity that hampers the detection of the background galaxies. Since the quality and stability of the PSFs of future instruments are unlikely to be much better than those of the WFPC2 PSF, a detailed understanding of optical/near-IR PSFs will be required to improve the results of the SFM in any Local Group galaxy.

R.A.G. thanks M. Guélin and the Institut de Radio-Astronomie Millimétrique (IRAM), Grenoble, for their hospitality during the preparation of this work; we are grateful to B. Turner and T. Dame for locating and faxing the Shapley 1951 reference within minutes of our request. R.A.G. was partially funded by the following grants: STScI AR-08360.01-A, NSF INT-9902653, and Mexican CONACyT I-35675E and 593-2000. L.L. acknowledges the support of CONACyT grant 592-2000.


text:
# Liquid-Gas Coexistence and Critical Behavior in Boxed Pseudo-Fermi Matter

## I Introduction

The possibility that a liquid-gas phase transition in finite nuclei may manifest itself via copious production of intermediate-mass fragments (IMF) in energetic heavy-ion reactions has driven both, theoretical and experimental studies of nuclear multifragmentation over more than a decade. A prominent role in theoretical considerations is played by the concept of a freezeout configuration.Gross (1997); Bondorf et al. (1995) This concept implies the existence of a definite volume, within which the system reaches a state close to thermal equilibrium. While the existence of an effective freezout volume may be debatable, the concept itself is useful for understanding the possible behavior of nuclear matter under various conditions. With such a didactic strategy in mind, and as an extension of earlier studies modelling the behavior of finite Fermi systems, Tõke and Schröder (2002, 1999); Tõke et al.  the present study considers pseudo-Fermi matter confined to a box of definite volume. It evaluates the isothermal behavior of such matter and its dependence on box volume and temperature. The utmost simplicity of the formalism allows one to gain insight into physical phenomena that may be obscured in more rigorous approaches.

## II Theoretical Formalism

The present study considers a scenario of nuclear pseudo-Fermi matter of mass number $`A`$, locked in a spherical box of a volume $`V`$ and kept at constant temperature $`T`$. In the proposed formalism, thermostatic properties of the matter are modelled by two equations, the isochoric caloric equation of state, and the zero-temperature equation of state. The equilibrium state of the system is then found as a function of $`V`$ and $`T`$, based on the requirement that the free energy of the system be minimal.

First, we consider homogeneous systems with constant density throughout their volume, i.e., systems in states where only one phase is present. The isochoric caloric equation of state for such a system is taken in a simple form adequate for low-temperature Fermi gases

$$E_{therm}^{}=aT^2,$$
(1)
where $`E_{therm}^{}`$ is the thermal energy and $`a`$ is the level density parameter. The latter parameter is assumed to depend on matter density as

$$a=a_o(\frac{\rho }{\rho _o})^{\frac{2}{3}},$$
(2)
in accordance with the low-temperature Fermi gas model.

It is worthwhile keeping in mind, that for Fermi gases, Eq. 1 is a good approximation only for temperatures that are small compared to the Fermi energy, i.e., for $`T<<E_F`$. Notably, for diluted or very hot Fermi systems, with the matter density approaching zero, $`\rho >0`$, or for high temperatures, $`T>\mathrm{}`$, the caloric equation of state approaches asymptotically that of a classical gas:

$$E_{therm}^{}=\frac{3}{2}AT$$
(3)
For the sake of simplicity and without loss of generality, the present study uses Eq. 1 over the full range of matter densities and temperatures considered. The term “pseudo-Fermi” matter is used to distinguish the matter considered here from true Fermi matter.

The second defining equation, the zero-temperature equation of state, expresses the compressional (potential) energy of the system as a function of the system volume or as a function of matter density. The present study adopts the harmonic approximation in the form used in the Expanding Emitting Source Model (EESM) Friedman (1988), whereby the in-medium nucleonic (potential) energy changes quadratically with the relative deviation of the actual matter density from the ground-state density, by an amount

$$ϵ_{compr}=ϵ_B(1\frac{\rho }{\rho _o})^2.$$
(4)
In Eq. 4, $`ϵ_{compr}`$ and $`ϵ_B`$ are compressional and ground-state binding energies per nucleon, respectively, and $`\rho `$ and $`\rho _o`$ are the actual and the ground-state matter densities, respectively. Equation 4 implies an effective ground-state incompressibility constant of $`K_o`$=-18$`ϵ_B`$. Assuming $`ϵ_B`$=-8 MeV, the effective incompressibility, including the effects of surface tension, equals $`K_o=144`$ MeV. Note, that within the harmonic approximation, for infinite nuclear matter characterized by $`ϵ_B16`$ MeV, the incompressibility constant is $`K_o288`$ MeV. The latter value places the present harmonic approximation “neutrally” between the currently considered limits of “soft” and “hard” equations of state for nuclear matter.

While the two defining equations 1 and 4 may be considered rather crude approximations, they do contain the essential physics responsible for first-order phase transitions and critical phenomena. Given these two equation, one can write expressions for all thermodynamic quantities characterizing the system, including the Helmholtz free energy $`F`$. The state of the system can then be found by minimizing the free energy, for any box volume $`V`$ and temperature $`T`$.

Based on Eq. 1, one can write for the entropy $`S`$ of a homogeneous system

$$S=_0^E\frac{1}{T}𝑑ϵ=2\sqrt{aE^{}}.$$
(5)
The free energy, $`F`$, for a homogeneous, single-phase system is given by

$$F=E_{total}^{}ST=E_{compr}^{}+E_{therm}^{}2aT^2=E_{compr}^{}aT^2.$$
(6)
In Eq. 6, for sake of simplicity, the free energy is given relative to a constant ground state energy of the system. Further, using Eqs. 2 and 4, Eq. 6 can be rewritten in a form revealing explicitly the important dependence on the matter density $`\rho `$, i.e.,

$$F=E_B(1\frac{\rho }{\rho _o})^2a_o(\frac{\rho }{\rho _o})^{\frac{2}{3}}T^2.$$
(7)
Eq. 7 allows one to write equations for the system pressure $`p`$ and the chemical potential $`\mu `$ as functions of volume (matter density) and temperature.

The pressure $`p`$ can be expressed generally as the negative partial derivative of the free energy $`F`$ with respect to volume $`V`$, at fixed number of nucleons $`A`$ and fixed temperature $`T`$, i.e.,

$$p=(\frac{F}{V})_{A,T}=\frac{1}{A}\rho ^2(\frac{F}{\rho })_{A,T}.$$
(8)
Thus, for the case of pseudo-Fermi matter considered here, one obtains based on Eqs. 7 and 8

$$p=2ϵ_B\rho _o(1\frac{\rho }{\rho _o})(\frac{\rho }{\rho _o})^2+\frac{2}{3}\alpha _o\rho _o(\frac{\rho }{\rho _o})^{\frac{1}{3}}T^2,$$
(9)
where $`ϵ_B`$ and $`\alpha _o`$ are the binding energy per nucleon and the level-density parameter per nucleon, respectively.

The chemical potential $`\mu `$ can be expressed generally as a partial derivative of free energy $`F`$ with respect to the number of nucleons $`A`$, taken at fixed volume $`V`$ and fixed temperature $`T`$, i.e.,

$$\mu =(\frac{F}{A})_{V,T}.$$
(10)
Using Eqs. 7 and 10, and noting further that $`\rho =A/V`$, $`E_B=Aϵ_B`$ and $`a_o=A\alpha _o`$, one obtains for the chemical potential

$$\mu =ϵ_B[14\frac{\rho }{\rho _o}+3(\frac{\rho }{\rho _o})^2]\frac{1}{3}\alpha _o(\frac{\rho }{\rho _o})^{\frac{2}{3}}T^2.$$
(11)
It is worth noting in Eq. 11 that, as a result of the requirement for $`V`$ to be constant, the magnitude of the chemical potential differs significantly from the value of the free energy per nucleon. For example, the contribution of thermal excitation to the chemical potential is only one-third of what constitutes the thermal part of the free energy per nucleon.

## III Liquid-Gas Coexistence

When confined to a thermostatic box, matter will eventually fill all of the available volume $`V`$ so as to minimize the free energy of the system. For non-interacting matter, the minimum free energy always corresponds to a uniform density distribution. This is generally not true for an interacting system. In particular, for interacting pseudo-Fermi matter at the minimum free energy, low-density and high-density phases coexist.

For a two-phase system of $`A`$ nucleons confined to volume $`V`$ at temperature $`T`$, the free energy can be expressed as a function of two variables, e.g., in terms of the volume of the gaseous phase, $`V_{gas}`$, and the number of nucleons contained in this phase, $`A_{gas}`$.

$$F=F_{gas}+F_{liquid},$$
(12)
Inserting for $`F_{gas}`$ and $`F_{liquid}`$ the expressions given by Eq. 7 for the corresponding numbers of nucleons $`A_{gas}`$ and $`A_{liquid}`$, respectively, one obtains

$$F_{gas}=A_{gas}ϵ_B(1\frac{A_{gas}}{V_{gas}\rho _o})^2A_{gas}\alpha _o(\frac{A_{gas}}{V_{gas}\rho _o})^{\frac{2}{3}}T^2$$
(13)
and

$$F_{liquid}=(AA_{gas})ϵ_B[1\frac{AA_{gas}}{(VV_{gas})\rho _o}]^2(AA_{gas})\alpha _o[\frac{AA_{gas}}{(VV_{gas})\rho _o}]^{\frac{2}{3}}T^2.$$
(14)
The condition for the minimum free energy can be expressed in form of two equations, reflecting requirements of dynamical and chemical equilibrium, respectively

$$(\frac{F}{V_{gas}})_{T,A_{gas}}=(\frac{F_{gas}}{V_{gas}})_{T,A_{gas}}(\frac{F_{liquid}}{V_{liquid}})_{T,A_{liquid}}=p_{gas}p_{liquid}=0$$
(15)
and

$$(\frac{F}{A_{gas}})_{T,V_{gas}}=(\frac{F_{gas}}{A_{gas}})_{T,V_{gas}}(\frac{F_{liquid}}{A_{liquid}})_{T,V_{liquid}}=\mu _{gas}\mu _{liquid}=0.$$
(16)
Results of a numerical minimization of the free energy of a two-phase system are shown in Figs. 1-4. Fig. 1 presents system isotherms, as predicted by the present formalism, in a representation of system pressure $`p`$ versus system volume $`V`$. Note that Eqs. 13 and 14 represent single-phase states as special cases. A pure liquid/gas state is thus among the possible outcomes of the calculations with $`A_{gas/liquid}=0`$. As seen in Fig. 1, the isotherms feature segments representing both, pure liquid and pure gas and, notably, the liquid-gas coexistence “plateaus”. The presence of coexistence plateaus does not come as a surprise as the harmonic interaction term of Eq. 4 has the salient characteristics of a Van der Waals interaction.

It is worth noting that the coexistence plateaus in this calculation result naturally from the actual minimization of the free energy and are not obtained via the well-known phenomenological Maxwell construct. Isotherms for hypothetical single-phase states are shown as dotted lines. These latter isotherms feature domains of spinodal instability characterized by negative compressibility. The dashed line in Fig. 1 illustrates the boundary of the liquid-gas coexistence domain. The “summit” point of this curve represents the critical point for the system and corresponds to a critical temperature of $`T_c10.0`$ MeV. At temperatures higher than $`T_c`$, the system can reside only in a single-phase, vapor state.

A different representation of the liquid-gas coexistence line is illustrated in Fig. 2. In this case, the temperature $`T`$ is plotted versus the matter density, for points along the boundary of the liquid-gas coexistence line (dashed line in Fig. 1). In the domain below this curve, the system is in a two-phase state. In this coexistence domain, the densities of gaseous and liquid phases at a given temperature are depicted by the intersection points of the coexistence curve with the horizontal line $`T`$=const. On the left shoulder of the curve and in the domain further left to it, the system is in a pure gaseous state, while on the right shoulder and in the domain further right to it, the system is in a pure liquid state. The difference between liquid and gaseous phases vanishes for temperatures above $`T=T_c=10.0`$ MeV.

The search for the minimum free energy for a given temperature $`T`$ and given volume $`V`$ yields, along with the volumes $`V_{gas/liquid}`$ and densities $`\rho _{gas/liquid}`$ of gaseous and liquid phases, the total excitation energy $`E_{total}^{}`$ of the boxed matter. This allows one to construct the isochoric caloric curves for the modelled system at different volumes of the confining box. A set of such caloric curves is illustrated in Fig. 3. As expected, these curves do not feature plateaus reported in various experiments,Natowitz et al. (2002) but rather inconspicuous kinks at the locations on the boundary of the coexistence domain.

## IV Critical Behavior

One of the salient features of critical behavior in Van der Waals systems is the presence of a singularity at the critical point. Here, when the system temperature T approaches the critical temperature $`T_c`$ from below, the difference between the densities of coexisting liquid and gaseous phases vanishes according to a power law

$$\rho _{liquid}\rho _{gas}=C(T_cT)^\beta ,$$
(17)
In Eq. 17, C is a constant and $`\beta `$ is the critical exponent. The magnitude of the critical exponent can be extracted conveniently by fitting a straight line to a double-logarithmic plot of $`(\rho _{liquid}\rho _{gas})`$ versus $`(T_cT)`$. Such a plot obtained in present calculations is shown in Fig. 4. The plot features, indeed, an approximately 2-MeV wide linear domain extending from the critical temperature $`T_c`$ down to lower temperatures. A linear fit to this domain allows one to extract the “coordinates” of the critical point of $`T_c`$= 10.0 MeV, and $`\rho _c/\rho _o`$=0.42. Furthermore, it yields a value of $`\beta `$=0.51 for the critical exponent.

Given the very schematic nature of the present formalism, the predicted characteristics of the critical point are well within a reasonable domain that can be inferred from more sophisticated, but less transparent calculations.

## V Discussion

A simple formalism has been presented that allows one to model the behavior of interacting pseudo-Fermi matter under the conditions of controlled volume and temperature. The formalism is shown to capture the essential physics of a first-order liquid-gas phase transition. The calculation demonstrates the characteristics of nuclear liquid-gas coexistence in a certain domain of system parameters and the presence of a critical point. While such characteristics are well expected, based on the the similarity of the utilized form of interaction to that of the Van der Waals interaction, the present formalism offers many didactical benefits. For example, due to its simplicity it offers a clear insight into physical phenomena it purports to describe and a relatively strict test-bench for possible qualitative or “hand-waving” arguments. The results obtained, while almost trivial, may alert one to possible challenges found by more complete models. As an example, it is evident from Fig. 1 that models of nuclear multifragmentation such as SMM Bondorf et al. (1995) and MMMC Gross (1997) tend to overpredict fragment multiplicities significantly. This is so, because these models assume that the liquid phase is always at ground-state density $`\rho _o`$. This forces the pressure of the surrounding gas within the freezeout volume to be higher than would be necessary to counter the lower pressure of the expanded liquid at equilibrium density. A higher gas pressure results in an overestimate of particle and fragment multiplicities.

The formalism presented here leaves ample room for further refinements such as a more strict modelling of the diluted Fermi matter, or the inclusion of isotopic effects. At any rate, it offers a convenient didactic tool to achieve a better understanding of nuclear thermodynamics.

###### Acknowledgements.

This work was supported by the U.S. Department of Energy grant No. DE-FG02-88ER40414.