text:
# Gravitational instability of finite isothermal spheres

## 1 Introduction

The thermodynamics of self-gravitating systems is a fascinating subject. It started with Antonov (1962)’s discovery that, when a self-gravitating system is confined within a box of radius $`R`$, no maximum entropy state can exist below a certain critical energy $`E=0.335GM^2/R`$. This intriguing result was further discussed by Lynden-Bell & Wood (1968) who conjectured that for $`E<0.335GM^2/R`$ the system would collapse and overheat. This is called “gravothermal catastrophe” or “Antonov instability”. Lynden-Bell & Wood have related this phenomenon to the very particular property of self-gravitating systems to possess negative specific heats. The gravothermal catastrophe picture has been confirmed by sophisticated numerical simulations (Larson 1970, Cohn 1980, Lynden-Bell & Eggleton 1980) and is expected to play a crucial role in the evolution of globular clusters. It is found that the collapse proceeds self-similarly (with power law behaviors) and that the central density becomes infinite in a finite time. This instability has been known as “core collapse” and many globular clusters have probably experienced core collapse (Binney & Tremaine 1987). In the case of dense clusters of compact stars (neutron stars or stellar mass black holes), the gravothermal catastrophe can lead to the formation of supermassive black holes of the right size to explain quasars and AGNs (Shapiro & Teukolsky 1995). Statistical mechanics is also relevant for collisionless stellar systems (e.g., elliptical galaxies, dark matter,…) undergoing a violent relaxation (Lynden-Bell 1967, Chavanis et al. 1996, Chavanis & Sommeria 1998, Chavanis 1998a,2001a). In particular, the inner regions of elliptical galaxies are close to isothermal and this is an important ingredient to understand de Vaucouleurs’ $`R^{1/4}`$ law (Hjorth & Madsen 1993).

On a theoretical point of view, the stability of isothermal spheres has been first investigated by Katz (1978) with a very powerful method extending Poincaré’s theory of linear series of equilibrium. He found that instability sets in precisely at the point of minimum energy. This stability analysis was reconsidered by Padmanabhan (1989) who studied the sign of the second variations of entropy and reduced the problem of stability to an eigenvalue equation. This leads to the same stability limit as Katz but the method of Padmanabhan provides in addition the form of the perturbation that induces instability at the critical point. It is found that this perturbation presents a “core-halo” structure.

The analysis of Padmanabhan (1989) was performed in the microcanonical ensemble in which the energy is fixed. The microcanonical ensemble is probably the most relevant for studying stellar systems like elliptical galaxies or globular clusters (Binney & Tremaine 1987). Indeed, apart from a slow evaporation, these systems can be assumed isolated so the evolution conserves energy $`E`$ and mass $`M`$. In addition, from the viewpoint of statistical mechanics, only the microcanonical ensemble is rigorously justified for non extensive systems, as discussed in the review of Padmanabhan (1990). However, it is always possible to define formally a canonical ensemble or a grand canonical ensemble in which the temperature is fixed instead of the energy. As suggested by de Vega et al. (1996a, 1996b) these ensembles may be suitable for describing the cold interstellar medium where the temperature is imposed by the cosmic background radiation at $`T3K`$ in the outer parts of galaxies, devoid of any star and heating sources (Pfenninger et al. 1994, Pfenninger & Combes 1994). In particular, by working out the statistical mechanics of the self-gravitating gas, de Vega et al. (1996a, 1996b) have shown that self-gravity can provide a dynamical mechanism to produce the fractal structure of the interstellar medium. They used the same approach to explain the fractal structure of the universe (de Vega et al. 1998), assuming that galaxies have reached a quasi-thermodynamical equilibrium like in the early work of Saslaw & Hamilton (1984).

For self-gravitating systems, it is well known that the thermodynamical ensembles do not coincide in the whole range of parameters (Padmanabhan 1990). Using toy models, Lynden-Bell & Lynden-Bell (1977) and Padmanabhan (1990) demonstrated that the region of negative specific heats allowed in the microcanonical ensemble is replaced by a phase transition in the canonical ensemble. This phase transition separates a dilute “gaseous” phase from a dense “collapsed” phase. Since these toy models are not very realistic, the self-gravitating gas was also studied in a meanfield approximation. In this viewpoint, an isothermal sphere is stable if and only if it is a local maximum of an appropriate thermodynamical potential (the entropy in the microcanonical ensemble and the free energy in the canonical ensemble). As expected from physical grounds, phase transition occurs when the gaseous sphere ceases to be a local maximum of this potential and becomes a saddle point. In this paper, we investigate the stability of isothermal gaseous spheres in the canonical ensemble by studying the sign of the second variations of the free energy. Our analysis is a direct extension of Padmanabhan (1989)’s approach in the microcanonical ensemble. The two studies therefore provide a unified description of the stability of isothermal spheres in the meanfield approximation in terms of thermodynamical potentials.

For a long time, the thermodynamics of self-gravitating systems was considered exclusively in a meanfield approach (or with toy models). However, recently, de Vega & Sanchez (2001a,b) have developed a rigorous statistical mechanics of self-gravitating systems by using field theoretical methods. In particular, they evidenced the existence of a thermodynamic limit in which the number of particles $`N`$ and the volume $`R^3`$ go to infinity keeping $`N/R`$ fixed. This very unusual thermodynamic limit proves to be appropriate to non extensive systems. By using Monte Carlo simulations and analytical calculations, they showed that the meanfield approximation correctly describes the thermodynamic limit except near the critical points where a phase transition occurs. They also derived the equation of state of a self-gravitating gas instead of assuming it, as is done in the meanfield treatments. Therefore, their work fully justifies the studies of previous authors and specifies their range of validity.

This paper is organized as follows. In section 2, we introduce a meanfield description of the system in the canonical ensemble and show that critical points of free energy $`J`$ at fixed temperature $`T`$ and mass $`M`$ correspond to isothermal spheres like those studied in the context of stellar structure (Chandraskhar 1942). We show that there is no global maximum of free energy. There is not even a local maximum of free energy in an unbounded domain: unbounded isothermal spheres have an infinite mass! We restrict therefore our analysis to the case of self-gravitating systems confined within a spherical box of radius $`R`$. In this case, there exists local maxima of $`J`$ (metastable states) if the normalized temperature $`\eta =\frac{\beta GM}{R}`$ is less than $`2.52`$ and the density contrast $`=\rho (0)/\rho (R)<32.1`$. Critical points of free energy with density contrast $`>32.1`$ are unstable saddle points. For $`\eta >2.52`$, there are not even critical points of free energy: in that case, the system is expected to undergo a phase transition and collapse. This “isothermal collapse” is the counterpart of the “gravothermal catastrophe” in the microcanonical ensemble (see Figs. 1-2).

In section 3, we study the sign of the second variations of free energy by using the methods of Padmanabhan (1989) introduced in the microcanonical ensemble. We show analytically that instability sets in precisely at the point of minimum temperature in agreement with the theorem of Katz (1978). The perturbation that induces instability at this point is calculated explicitly; it has not a “core-halo” structure contrary to what happens in the microcanonical ensemble.

In section 4, we study Jeans type gravitational instability of isothermal gaseous spheres described by Navier-Stokes equations. The introduction of a container removes the problems associated with an infinite homogeneous medium and avoids the Jeans swindle. We show analytically the equivalence between dynamical stability and thermodynamical stability and the fact that the stability of isothermal gas spheres does not depend on the viscosity. This confirms the findings of Semelin et al. (2001) who used numerical methods or approximations. We also give a simpler derivation of the geometric hierarchy of scales inducing instability discovered by these authors using sophisticated renormalization group technics (Semelin et al. 1999). This provides a more illuminating interpretation of their results.

In section 5, we make speculations about the fragmentation and the fractal structure of an isothermal self-gravitating gas. We distinguish between the Jeans length $`L_J`$ defined with the mean density and the King’s length $`L_K`$ defined with the central density. If we fix the Jeans length, instability occurs for $`R>\sqrt{\frac{2.52\mathrm{}}{3}}L_J`$ and is marked by the absence of critical point of free energy (i.e., hydrostatic equilibrium) above this threshhold. In that case the system is expected to collapse without fragmenting. If we fix the King’s length (or core size), a first instability occurs for $`R=\frac{8.99\mathrm{}}{3}L_K`$. Above this threshold critical points of free energy still exist but they are unstable saddle points. Secondary instabilities occur at larger box radii that asymptotically follow a geometric progression $`R_n[10.74\mathrm{}]^nL_K`$. The density profiles that trigger these high order modes of instability are calculated analytically. They present more and more oscillations whose zeros also follow a geometric progression $`r_n[10.74\mathrm{}]^nR`$. The profile that destabilizes the singular isothermal sphere has an infinite number of nodes! Each oscillation can be interpreted as a ‘germ’ in the langage of phase transition and the above picture suggests that the system will fragment into a series of ‘clumps’. It is expected that these ‘clumps’ will evolve by achieving higher and higher density contrasts, and finally fragment in turn into substructures. This yields a hierarchy of structures fitting one into each other in a self-similar way. This picture is given further support by the fact that both the domain sizes inducing instability and the zeros of the perturbation profile in each domain follow a geometric progression with the same ratio. This double-geometric progression may explain in a natural way the fractal structure of a self-gravitating gas like the interstellar medium and the large scale structures of the universe. This gives further support to the interpretation of de Vega et al. (1996a, 1996b, 1998) who first realized the importance played by self-gravity in building a fractal distribution of matter.

## 2 Thermodynamical stability of self-gravitating systems in the canonical ensemble

### 2.1 The free energy

Consider a system of $`N`$ particles, each of mass $`m`$, interacting via Newtonian gravity. The particles can be galaxies, stars, atoms, etc… We assume that the system is non rotating and non expanding. Let $`f(𝐫,𝐯,t)`$ denote the distribution function of the system, i.e. $`f(𝐫,𝐯,t)d^3𝐫d^3𝐯`$ gives the mass of particles whose position and velocity are in the cell $`(𝐫,𝐯;𝐫+d^3𝐫,𝐯+d^3𝐯)`$ at time $`t`$. The integral of $`f`$ over the velocity determines the spatial density

(1)
$$\rho =fd^3𝐯.$$

On the other hand, in the meanfield approximation, the total mass and the total energy of the system can be expressed as

(2)
$$M=Nm=\rho d^3𝐫,$$

(3)
$$E=\frac{1}{2}fv^2d^3𝐫d^3𝐯+\frac{1}{2}\rho \mathrm{\Phi }d^3𝐫=K+W,$$

where $`K`$ is the kinetic energy and $`W`$ the potential energy. The gravitational potential $`\mathrm{\Phi }`$ is related to the star density by the Newton-Poisson equation

(4)
$$\mathrm{\Delta }\mathrm{\Phi }=4\pi G\rho .$$

Finally, the Boltzmann entropy is given by the standard formula

(5)
$$S=k\frac{f}{m}\mathrm{ln}\frac{f}{m}d^3𝐫d^3𝐯,$$

which can be obtained by counting the number of microstates corresponding to a given macrostate and taking the logarithm of this number (Ogorodnikov 1965).

We shall work in the canonical ensemble in which the temperature is fixed, allowing the energy to fluctuate. In that case, the relevant thermodynamical potential is the Massieu function related to the Helmholtz free energy $`F=ETS`$ by $`J=\frac{1}{T}F`$. Hence

(6)
$$J=S\frac{1}{T}E.$$

At equilibrium, the system is expected to be in a state that maximizes the Massieu function (6) for a fixed total mass $`M`$ (in the following we shall call $`J`$ the free energy, although this is only the free energy up to a negative proportionality factor).

### 2.2 The isothermal gaseous spheres

Following Padmanabhan (1989)’s procedure, we start to maximize the free energy $`J`$ for a given density field $`\rho (𝐫)`$. This yields the Maxwell-Boltzmann distribution

(7)
$$f=\left(\frac{m}{2\pi kT}\right)^{3/2}\rho (𝐫)e^{\frac{mv^2}{2kT}},$$

which is a global maximum of $`J`$ with the previous constraint. Substituting this optimal distribution function in Eqs. (3)-(5), we can express the energy and the entropy in terms of the spatial density in the form

(8)
$$E=\frac{3}{2}NkT+\frac{1}{2}\rho \mathrm{\Phi }d^3𝐫,$$

(9)
$$\frac{S}{k}=\frac{3N}{2}+\frac{3N}{2}\mathrm{ln}\left(\frac{2\pi kT}{m}\right)\frac{\rho }{m}\mathrm{ln}\frac{\rho }{m}d^3𝐫.$$

We can now determine the variations of $`J`$ around a given density profile $`\rho (𝐫)`$. To second order in the expansion we get

$`\delta J={\displaystyle \frac{k}{m}}{\displaystyle \delta \rho \left(\mathrm{ln}\frac{\rho }{m}+1\right)d^3𝐫}{\displaystyle \frac{k}{m}}{\displaystyle \frac{(\delta \rho )^2}{2\rho }d^3𝐫}`$
(10) $`{\displaystyle \frac{1}{T}}{\displaystyle \mathrm{\Phi }\delta \rho d^3𝐫}{\displaystyle \frac{1}{2T}}{\displaystyle \delta \rho \delta \mathrm{\Phi }d^3𝐫}.`$
Introducing a Lagrange multiplier $`\alpha `$ to satisfy the conservation of mass, the condition that $`J`$ is an extremum is written (to first order)

(11) $`0=\delta J\alpha \delta M={\displaystyle d^3𝐫\left[\frac{\mathrm{\Phi }}{T}+\frac{k}{m}\left(\mathrm{ln}\frac{\rho }{m}+1\right)+\alpha \right]\delta \rho }.`$
This condition must be satisfied for any variations $`\delta \rho `$. This yields the Boltzmann distribution

(12)
$$\rho =Ae^{\beta \mathrm{\Phi }},$$

where we have set

(13)
$$\beta =\frac{m}{kT}.$$

Therefore, the Boltzmann distribution (12) is a critical point of free energy. This does not insure, however, that it is a maximum of $`J`$. It is not even clear that the extremum problem leading to Eq. (12) has a solution. Indeed, the gravitational potential that appears in Eq. (12) must be determined self-consistently by solving the mean field equation

(14)
$$\mathrm{\Delta }\mathrm{\Phi }=4\pi GAe^{\beta \mathrm{\Phi }},$$

obtained by substituting the density (12) in the Poisson equation (4), and relating the Lagrange multiplier $`A`$ to the constraint $`M`$. As we shall see, this problem does not always have a solution. When a solution exists, we must consider the sign of the second order variations $`\delta ^2J`$ to determine whether it is a maximum or not.

### 2.3 Absence of global maximum of free energy

It is easy to show that there is no global maximum of free energy. To prove this result, we just need to consider a homogeneous sphere of mass $`M`$ and radius $`R`$. The total energy of this sphere is

(15)
$$E=\frac{3}{2}NkT\frac{3GM^2}{5R},$$

and its free energy

(16)
$$J=\frac{3}{2}Nk\mathrm{ln}\left(\frac{2\pi kT}{m}\right)Nk\mathrm{ln}\left(\frac{3N}{4\pi R^3}\right)+\frac{3GM^2}{5TR}.$$

From the above formula it is clear that $`J`$ goes to $`+\mathrm{}`$ when we spread the density to infinity ($`R+\mathrm{}`$) or when we contract the system to a point ($`R0`$). Since the mass is conserved in this process no global maximum of free energy can exist.

In addition, there is not even a local maximum of free energy in an unbounded domain: unbounded isothermal spheres have an infinite mass. For non rotating systems, the equilibrium states are expected to be spherically symmetric. In that case, the Boltzmann-Poisson equation (14) can be rewritten

(17)
$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d\mathrm{\Phi }}{dr}\right)=4\pi GAe^{\beta \mathrm{\Phi }}.$$

This equation has been studied extensively in the context of isothermal gaseous spheres (Chandrasekhar 1942). We can check by direct substitution that the distribution

(18)
$$\mathrm{\Phi }_s(r)=\frac{1}{\beta }\mathrm{ln}(2\pi G\beta Ar^2),\rho _s(r)=\frac{1}{2\pi G\beta r^2},$$

is an exact solution of Eq. (17) known as the singular isothermal sphere (Binney & Tremaine 1987). Since $`\rho r^2`$ at large distances, the total mass of the system $`M=_0^+\mathrm{}\rho 4\pi r^2𝑑r`$ is infinite! More generally, we can show that any solution of the meanfield equation (17) behaves like the singular sphere as $`r+\mathrm{}`$ and has therefore an infinite mass (Chandrasekhar 1942).

We shall avoid the infinite mass problem by confining artificially the system within a spherical box of radius $`R`$. It is only under this simplifying assumption that a rigorous thermodynamics of self-gravitating systems can be carried out. Physically, the maximum cut-off determines the scale at which other processes intervene to limitate the spatial extent of the system. It is clear from Eq. (16) that the introduction of an upper cut-off does not remove the absence of global maximum of free energy (the system can always increase $`J`$ by collapsing to a point). However, the presence of a wall at radius $`R`$ allows the existence of local maxima of free energy (metastable states).

### 2.4 The equilibrium phase diagram

The microcanonical and canonical ensembles yield the same critical points, i.e. the critical points of entropy at fixed mass and energy and the critical points of free energy at fixed mass and temperature coincide. Only the onset of instability will differ from one ensemble to the other. The thermodynamical parameters for bounded isothermal spheres in the meanfield approximation have been calculated by Lynden-Bell & Wood (1968) and we shall directly use their results.

To that purpose, we introduce the function $`\psi =\beta (\mathrm{\Phi }\mathrm{\Phi }_0)`$ where $`\mathrm{\Phi }_0`$ is the gravitational potential at $`r=0`$. Then, the density field can be written

(19)
$$\rho =\rho _0e^\psi ,$$

where $`\rho _0`$ is the central density. Introducing the notation $`\xi =(4\pi G\beta \rho _0)^{1/2}r`$, the Boltzmann-Poisson equation (17) reduces to the standard Emden form (Chandrasekhar 1942)

(20)
$$\frac{1}{\xi ^2}\frac{d}{d\xi }\left(\xi ^2\frac{d\psi }{d\xi }\right)=e^\psi .$$

Eq. (20) has a simple analytic solution, the singular sphere

(21)
$$e^{\psi _s}=\frac{2}{\xi ^2}.$$

The regular solutions of Eq. (20) must satisfy the boundary conditions

(22)
$$\psi =\psi ^{}=0,$$

when $`\xi =0`$. These solutions must be computed numerically. However, their asymptotic behaviors are well-known

(23)
$$\psi =\frac{1}{6}\xi ^2\frac{1}{120}\xi ^4+\frac{1}{1890}\xi ^6+\mathrm{}(\xi 0),$$

(24)
$$e^\psi =\frac{2}{\xi ^2}\left\{1+\frac{A}{\xi ^{1/2}}\mathrm{cos}\left(\frac{\sqrt{7}}{2}\mathrm{ln}\xi +\delta \right)\right\}(\xi +\mathrm{}).$$

The curve (24) intersects the singular solution (21) infinitely often at points that asymptotically increase geometrically in the ratio $`1:e^{2\pi /\sqrt{7}}=1:10.74\mathrm{}`$. This property will have important physical consequences in section 5.

In the case of bounded isothermal systems, the solutions of Eq. (20) are terminated by the box at different radii given by $`\alpha =(4\pi G\beta \rho _0)^{1/2}R`$. Lynden-Bell & Wood (1968) show that the reduced temperature and reduced energy can be expressed in terms of $`\alpha `$ by

(25)
$$\eta \frac{\beta GM}{R}=\alpha \psi ^{}(\alpha ),$$

(26)
$$\mathrm{\Lambda }\frac{ER}{GM^2}=\frac{3}{2}\frac{1}{\alpha \psi ^{}(\alpha )}\frac{e^{\psi (\alpha )}}{\psi ^{}(\alpha )^2}.$$

For each value of the normalized inverse temperature $`\frac{\beta GM}{R}`$, we can solve Eq. (25) to get $`\alpha `$. Substituting the result in Eq. (26), we deduce the corresponding value of the normalized energy $`\frac{ER}{GM^2}`$. We can then determine the equilibrium diagram $`E\beta `$ (Figs. 1-2). The critical points of entropy or free energy form a spiral. This curve is parametrized by $`\alpha `$ that goes from $`0`$ to $`+\mathrm{}`$ as we spiral inward. Instead of the parameter $`\alpha `$, it may be more relevant to introduce the density contrast

(27)
$$=\frac{\rho _0}{\rho (R)}=e^{\psi (\alpha )},$$

which gives a more physical parametrization of the solutions. For $`\alpha 0`$, we can use the asymptotic behavior $`\psi \frac{1}{6}\alpha ^2`$ of the potential at the origin \[Eq. (23)\] and we find $`\frac{\beta GM}{R}\frac{\alpha ^2}{3}0`$, $`\frac{ER}{GM^2}\frac{9}{2\alpha ^2}\mathrm{}`$ and $`1`$. This corresponds to high temperatures $`T+\mathrm{}`$. In this case, self-gravity is negligible with respect to thermal motion and the system behaves like an ordinary gas with uniform density. We check, by eliminating $`\alpha `$ from the two foregoing relations, that $`E=\frac{3M}{2\beta }`$. This is indeed the equation of state for an ideal gas without self-gravity. When we proceed along the spiral, the density contrast increases monotonically and goes to $`+\mathrm{}`$ at the center of the spiral. Using the asymptotic expansion $`e^{\psi (\alpha )}\frac{2}{\alpha ^2}`$ of the potential at large distances \[Eq. (24)\], we find $`\frac{\beta GM}{R}2`$, $`\frac{ER}{GM^2}\frac{1}{4}`$ and $`\frac{\alpha ^2}{2}+\mathrm{}`$. We can check by a direct calculation that the ending point $`(1/4,2)`$ of the spiral corresponds to the singular sphere (18).

### 2.5 The Milne variables

It will be convenient in the following to introduce the Milne variables $`(u,v)`$ defined by (Chandrasekhar 1942):

(28)
$$u=\frac{\xi e^\psi }{\psi ^{}}\mathrm{and}v=\xi \psi ^{}.$$

Taking the logarithmic derivative of $`u`$ and $`v`$ with respect to $`\xi `$ and using Eq. (20), we get

(29)
$$\frac{1}{u}\frac{du}{d\xi }=\frac{1}{\xi }(3vu),$$

(30)
$$\frac{1}{v}\frac{dv}{d\xi }=\frac{1}{\xi }(u1).$$

Taking the ratio of the foregoing equations, we find that the variables $`u`$ and $`v`$ are related to each other by a first order differential equation

(31)
$$\frac{u}{v}\frac{dv}{du}=\frac{u1}{u+v3}.$$

Therefore, by using the Milne variables, the degree of the meanfield equation (20) has been reduced from two to one. As discussed extensively by Chandrasekhar (1942), this property is related to the homology invariance of the solutions of the Emden equation. The solution curve in the $`(u,v)`$ plane is well known and is ploted in Fig. 3. The curve is parametrized by $`\xi `$. Starting from the point $`(u,v)=(3,0)`$ corresponding to $`\xi =0`$, the solution curve spirals indefinitely around the point $`(u_s,v_s)=(1,2)`$, corresponding to the singular sphere, as $`\xi `$ tends to infinity. All isothermal spheres must necessarily lie on this curve. For bounded isothermal spheres, $`\xi `$ must be terminated at the box radius $`\alpha `$. Clearly, the spiral behavior of the $`(u,v)`$ curve (and also the curve of Figs. 1-2) can be ascribed to the oscillating behavior of the solution (24) as $`\xi +\mathrm{}`$.

It turns out that the normalized temperature and energy can be expressed very simply in terms of the values of $`u`$ and $`v`$ at the normalized box radius $`\alpha `$. Indeed, writing $`u_0=u(\alpha )`$ and $`v_0=v(\alpha )`$ and using Eqs. (25)-(26), we get

(32)
$$\mathrm{\Lambda }=\frac{1}{v_0}\left(\frac{3}{2}u_0\right),$$

(33)
$$\eta =v_0.$$

The relation (32) was previously noticed by Padmanabhan (1989). The intersection between the lines defined by Eqs. (32)-(33) and the spiral in the $`(u,v)`$ plane determines the value of $`\alpha `$ corresponding to a given energy or temperature. As noticed by Padmanabhan, for $`\mathrm{\Lambda }>\mathrm{\Lambda }_c=0.335`$ there is no intersection. Thus, in the microcanonical ensemble, no isothermal sphere can exist if $`\frac{ER}{GM^2}>0.335`$. This result can also be read from Fig. 1 and was discovered by Antonov (1962). For sufficiently small energies, the system is expected to collapse and overheat in agreement with the “gravothermal catastrophe” picture (Lynden-Bell & Wood 1968). Similarly, considering Eq. (33), we find that there is no intersection for $`\frac{\beta GM}{R}>v_{max}=2.52`$ (where $`v_{max}`$ corresponds to $`u=1`$). This result can also be read from Fig. 2. In the canonical ensemble, a gaseous sphere is expected to collapse below a critical temperature $`kT_c=\frac{GmM}{2.52R}`$ (Lynden-Bell & Wood 1968).

### 2.6 The stability analysis of Katz

More precisely, it is possible to prove the following stability result: when a self-gravitating system is confined within a box and maintained at a constant temperature $`T`$, local maxima of free energy exist only for $`\frac{\beta GM}{R}2.52`$; they have a density contrast $`=\rho (0)/\rho (R)<32.2`$ (and $`\alpha <8.99`$). Critical points of free energy with density contrast $`>32.2`$ (and $`\alpha >8.99`$) are unstable saddle points. For $`\frac{\beta GM}{R}>2.52`$, there are no critical points of free energy.

The stability of the solutions can be deduced from the topology of the $`E\beta `$ curve by using the method of Katz (1978) who has extended Poincaré’s theory of linear series of equilibrium. The parameter conjugate to the free energy with respect to the inverse temperature $`\beta `$ is $`E=(\frac{J}{\beta })_{M,R}`$. Then, if we plot $`E`$ as a function of $`\beta `$ (we just need to rotate Fig. 2 by $`90^o`$) we have the following results: (i) a change of stability can occur only at a limit point where $`\beta `$ is an extremum ($`dE/d\beta `$ infinite) (ii) a mode of stability is lost when the curve rotates clockwise and gained otherwise. Now, we know that for $`T`$ sufficiently large the solutions are stable because, in this limit, self-gravity is negligible and the system behaves like an ordinary gas. From point (i) we conclude that the solutions with $`=\rho (0)/\rho (R)<32.2`$ are stable. As the curve spirals inward for $`>32.2`$, more and more modes of stability are lost. In this respect, the singular sphere, at the end of the spiral, is the most unstable solution <sup>1</sup><sup>1</sup>1Another application of Katz theorem where the spiral unwinds and the stability is regained is given by Chavanis & Sommeria (1998) for self-gravitating fermions.. It can be noted that instability sets in precisely when the specific heat $`C=dE/dT`$ becomes negative (in the range $`32.2<<709`$). General considerations indicate that negative specific heats are forbidden in the canonical ensemble (Padmanabhan 1990). By contrast, thermally isolated systems can have negative specific heats: in the microcanonical ensemble, the solutions with $`32.2<<709`$ are stable (Katz 1978, Padmanabhan 1989). This clearly demonstrates that the statistical ensembles do not coincide for self-gravitating systems: the region of negative specific heats in the microcanonical ensemble is replaced by a phase transition (an “isothermal collapse”) in the canonical ensemble (Padmanabhan 1990). It is at the verge of this phase transition that the meanfield approximation ceases to be valid as demonstrated by de Vega et al. (2001a, 2001b).

### 2.7 Analogy with critical phenomena

In this section, we study some analogies with the theory of critical phenomena. Further analogies are discussed in Chavanis et al. (2001) where a simple dynamical model for studying phase transitions in self-gravitating systems is proposed.

The specific heat can be written

(34)
$$C=\frac{dE}{dT}=\frac{k}{m}\beta ^2\frac{dE}{d\beta }=Nk\eta ^2\frac{d\mathrm{\Lambda }}{d\eta }.$$

Using Eqs. (32)(33)(29)(30), we easily get

(35)
$$\frac{d\eta }{d\alpha }=\frac{v_0}{\alpha }(u_01),$$

(36)
$$\frac{d\mathrm{\Lambda }}{d\alpha }=\frac{1}{2\alpha v_0}(4u_0^2+2u_0v_011u_0+3)$$

Therefore, the specific heat is determined as a function of $`\alpha `$ by the expression

(37)
$$C=Nk\frac{4u_0^2+2u_0v_011u_0+3}{2(u_01)}.$$

Of course, $`C\mathrm{}`$ if $`d\eta =0`$ and $`C0`$ if $`d\mathrm{\Lambda }=0`$. On the other hand, for the ideal gas $`(u_0,v_0)=(3,0)`$, we recover the well-known result $`C=\frac{3}{2}Nk`$. The specific heat is ploted as a function of $`\alpha `$ on Fig. 4.

At the critical point $`\eta =\eta _c`$ at which $`C`$ diverges, we write $`\alpha =\alpha _c+ϵ`$ and expand $`\eta `$ and $`\mathrm{\Lambda }`$ in terms of $`ϵ`$, using Eqs. (35)-(36) and $`u_0(\alpha _c)=1`$, $`v_0(\alpha _c)=\eta _c`$. To lowest order, we have

(38)
$$\eta \eta _c=\left(\frac{d^2\eta }{d\alpha ^2}\right)_{\alpha _c}\frac{ϵ^2}{2}=\frac{\eta _c}{2\alpha _c^2}(2\eta _c)ϵ^2,$$

(39)
$$\mathrm{\Lambda }\mathrm{\Lambda }_0=\left(\frac{d\mathrm{\Lambda }}{d\alpha }\right)_{\alpha _c}ϵ=\frac{\eta _c2}{\alpha _c\eta _c}ϵ,$$

where $`\mathrm{\Lambda }_00.197`$ is the value of $`\mathrm{\Lambda }`$ at $`\alpha =\alpha _c`$. Eliminating $`ϵ`$ from these expressions, we find the relation between $`\eta `$ and $`\mathrm{\Lambda }`$ close to the critical point

(40)
$$\eta _c\eta =\frac{\eta _c^3}{2(\eta _c2)}(\mathrm{\Lambda }\mathrm{\Lambda }_0)^2.$$

Then, using Eq. (34), we find explicitely

(43) $`C=\{\begin{array}{cc}Nk\frac{\eta _c2}{\eta _c}(\mathrm{\Lambda }\mathrm{\Lambda }_0)^1,& \\ \pm Nk\left[\frac{\eta _c(\eta _c2)}{2}\right]^{1/2}(\eta _c\eta )^{1/2}.& \end{array}`$
The specific heat diverges as a function of the energy with an exponent $`\nu =1`$ and as a function of the temperature with an exponent $`\nu ^{}=1/2`$ (see Figs. 5-6). It becomes negative for $`\mathrm{\Lambda }_0<\mathrm{\Lambda }<\mathrm{\Lambda }_c`$. This region of negative specific heats is allowed in the microcanonical ensemble but not in the canonical ensemble.

It is also of interest to determine the behavior of the central density $`\rho _0`$ near the critical point (C. Sire, private communication). The central density is related to $`\alpha `$ by the relation $`\alpha =(4\pi G\beta \rho _0)^{1/2}R`$ or, alternatively,

(44)
$$\rho _0=\frac{\alpha ^2M}{4\pi R^3\eta }.$$

To first order in $`ϵ`$, we have

(45)
$$\rho _0\rho _0^c=\frac{2\alpha _cM}{4\pi R^3\eta _c}ϵ.$$

where $`\rho _0^c=\alpha _c^2M/(4\pi R^3\eta _c)`$ is the value of the central density at the critical point $`\eta _c`$. Using Eqs. (38)-(39), we get

(48) $`{\displaystyle \frac{\rho _0^c\rho _0}{\rho _0^c}}=\{\begin{array}{cc}\frac{2\eta _c}{\eta _c2}(\mathrm{\Lambda }\mathrm{\Lambda }_0),& \\ \pm \left[\frac{8}{\eta _c2}\left(1\frac{\eta }{\eta _c}\right)\right]^{1/2}.& \end{array}`$
Quite remarkably, this stationary result can also be derived from dynamical considerations (see Chavanis et al. 2001). Note also that, using Eqs. (27)(28), we find $`_c=\alpha _c^2/\eta _c`$, which implies $`\rho ^c(R)=M/(4\pi R^3)`$.

## 3 An extension of Padmanabhan’s stability analysis

### 3.1 The condition of thermodymamical stability

The stability of the gaseous spheres can be studied by extending Padmanabhan (1989)’s stability analysis to the canonical ensemble. This will provide the profiles of the perturbations that induce instability.

A critical point of free energy is a local maximum if, and only if, the second variations

(49)
$$\delta ^2J=\frac{k}{m}\frac{(\delta \rho )^2}{2\rho }d^3𝐫\frac{1}{2T}\delta \rho \delta \mathrm{\Phi }d^3𝐫,$$

are negative for any variation $`\delta \rho `$ that conserves mass to first order. This conservation law imposes

(50)
$$\delta \rho d^3𝐫=0.$$

We shall restrict ourselves to spherically symmetric perturbations since only spherically symmetric perturbations can induce instability for non rotating systems. Following Padmanabhan (1989), we introduce the variable $`q`$ defined by

(51)
$$\delta \rho =\frac{1}{4\pi r^2}\frac{dq}{dr}.$$

Physically, $`q`$ represents the mass perturbation $`q(r)\delta M(r)=_0^r4\pi r^{}_{}{}^{}2\delta \rho (r^{})𝑑r^{}`$ within the sphere of radius $`r`$. The boundary conditions on $`q`$ are thus

(52)
$$q(0)=q(R)=0.$$

Substituting the foregoing expression for $`q(r)`$ in Eq. (49), we obtain

(53)
$$\delta ^2J=\frac{1}{2T}_0^R\frac{dq}{dr}\delta \mathrm{\Phi }𝑑r\frac{k}{m}_0^R\frac{1}{8\pi \rho r^2}\left(\frac{dq}{dr}\right)^2𝑑r.$$

Integrating by parts and using the boundary conditions (52), we get

(54)
$$\delta ^2J=\frac{1}{2T}_0^Rq\frac{d\delta \mathrm{\Phi }}{dr}𝑑r+\frac{k}{m}_0^Rq\frac{d}{dr}\left(\frac{1}{8\pi \rho r^2}\frac{dq}{dr}\right)𝑑r.$$

Then, using Gauss theorem in the form

(55)
$$\frac{d\delta \mathrm{\Phi }}{dr}=\frac{Gq}{r^2},$$

we find

(56)
$$\delta ^2J=\frac{1}{2T}_0^R\frac{Gq^2}{r^2}𝑑r+\frac{k}{m}_0^Rq\frac{d}{dr}\left(\frac{1}{8\pi \rho r^2}\frac{dq}{dr}\right)𝑑r,$$

or, equivalently,

(57)
$$\delta ^2J=\frac{1}{2}_0^R𝑑rq\left[\frac{G}{Tr^2}+\frac{k}{m}\frac{d}{dr}\left(\frac{1}{4\pi \rho r^2}\frac{d}{dr}\right)\right]q.$$

The second variations of free energy can be positive (implying instability) only if the differential operator which occurs in the integral has positive eigenvalues. We need therefore to consider the eigenvalue problem

(58)
$$\left[\frac{k}{m}\frac{d}{dr}\left(\frac{1}{4\pi \rho r^2}\frac{d}{dr}\right)+\frac{G}{Tr^2}\right]q_\lambda (r)=\lambda q_\lambda (r),$$

with $`q_\lambda (0)=q_\lambda (R)=0`$. If all the eigenvalues $`\lambda `$ are negative, then the critical point is a maximum  of free energy. If at least one eigenvalue is positive, the critical point is an unstable saddle point. The point of marginal stability $`\eta =\eta _c`$ is determined by the condition that the largest eigenvalue is equal to zero. We thus have to solve the differential equation

(59)
$$\frac{k}{m}\frac{d}{dr}\left(\frac{1}{4\pi \rho r^2}\frac{dF}{dr}\right)+\frac{GF(r)}{Tr^2}=0,$$

with $`F(0)=F(R)=0`$. Note that this equation corresponds to that found by Padmanabhan (1989) in the microcanonical ensemble with $`V=0`$.

### 3.2 The point of marginal stability

Introducing the dimensionless variables of section 2.4, we can rewrite the foregoing equation in the form

(60)
$$\frac{d}{d\xi }\left(\frac{e^\psi }{\xi ^2}\frac{dF}{d\xi }\right)+\frac{F(\xi )}{\xi ^2}=0,$$

with $`F(0)=F(\alpha )=0`$. As shown by Padmanabhan (1989), this equation can be solved without solving explicitly the Emden equation (20). Let us introduce the differential operator

(61)
$$=\frac{d}{d\xi }\left(\frac{e^\psi }{\xi ^2}\frac{d}{d\xi }\right)+\frac{1}{\xi ^2}.$$

Using Eq. (20), it is readily established that

(62) $`(\xi ^2\psi ^{})={\displaystyle \frac{d}{d\xi }}\left({\displaystyle \frac{e^\psi }{\xi ^2}}{\displaystyle \frac{d}{d\xi }}(\xi ^2\psi ^{})\right)+\psi ^{}={\displaystyle \frac{d}{d\xi }}(e^\psi \times e^\psi )+\psi ^{}=\psi ^{}.`$
Similarly, we find

$`(\xi ^3e^\psi )={\displaystyle \frac{d}{d\xi }}\left({\displaystyle \frac{e^\psi }{\xi ^2}}{\displaystyle \frac{d}{d\xi }}(\xi ^3e^\psi )\right)+\xi e^\psi `$
(63) $`={\displaystyle \frac{d}{d\xi }}(3\xi \psi ^{})+\xi e^\psi =\xi \psi ^{\prime \prime }\psi ^{}+\xi e^\psi =\psi ^{}.`$
These identities suggest to seek a solution of Eq. (60) in the form

(64)
$$F(\xi )=c_1\xi ^3e^\psi +c_2\xi ^2\psi ^{}$$

where $`c_1`$ and $`c_2`$ are constants. Substituting the foregoing expression for $`F(\xi )`$ in the differential equation (60), we find that $`c_1+c_2=0`$. The function $`F(\xi )`$ can therefore be reexpressed as

(65)
$$F(\xi )=c_1(\xi ^3e^\psi \xi ^2\psi ^{}).$$

The boundary condition $`F(0)=0`$ is automatically satisfied (we can show furthermore that Eq. (65) is the only solution consistent with this boundary condition). The condition $`F(\alpha )=0`$ will determine the point on the spiral of Fig. 2 at which the gaseous spheres become unstable (i.e., the point of marginal stability). From Eq. (65), we get

(66)
$$\alpha ^3e^{\psi (\alpha )}\alpha ^2\psi ^{}(\alpha )=0.$$

This can be expressed in terms of the Milne variables (28) as

(67)
$$u_0=1.$$

Considering Eq. (31), we see that the condition (67) determines the points for which $`v`$ is extremum ($`dv/du=0`$). Since $`\eta =v_0`$ according to Eq. (33), we deduce that a change of stability occurs each time that the temperature is extremum (see Fig. 7) in complete agreement with the theorem of Katz (1978). In particular, the first instability (for which the largest eigenvalue $`\lambda `$ becomes positive) corresponds to $`\eta _c=v_{max}=2.52`$ (for this value $`\alpha =8.99`$). We show in the next section that the secondary instabilities occur at values of $`\alpha `$ that follow asymptotically a geometric progression. The profiles of the perturbations that induce these instabilities are discussed in section 3.4.

### 3.3 Secondary instabilities

A mode of stability is lost each time that $`u_0=u(\alpha )=1`$. Now, it is easy to have an analytical estimate of $`u(\xi )`$ for large values of $`\xi `$. Introducing the asymptotic behavior (24) of $`\psi `$ in the Milne variable $`u`$ defined by Eq. (28) we get for $`\xi +\mathrm{}`$:

(68) $`u={\displaystyle \frac{1+\frac{A}{\xi ^{1/2}}\mathrm{cos}\left(\frac{\sqrt{7}}{2}\mathrm{ln}\xi +\delta \right)}{1+\frac{A}{4\xi ^{1/2}}\left[\sqrt{7}\mathrm{sin}\left(\frac{\sqrt{7}}{2}\mathrm{ln}\xi +\delta \right)+\mathrm{cos}\left(\frac{\sqrt{7}}{2}\mathrm{ln}\xi +\delta \right)\right]}}.`$
Therefore, the condition $`u_0=1`$ corresponds to values of $`\alpha `$ satisfying

(69)
$$\mathrm{tan}\left(\frac{\sqrt{7}}{2}\mathrm{ln}\alpha +\delta \right)=\frac{3}{\sqrt{7}}(\alpha +\mathrm{}),$$

or equivalently

(70) $`{\displaystyle \frac{\sqrt{7}}{2}}\mathrm{ln}\alpha _n+\delta =\mathrm{arctan}(3/\sqrt{7})+n\pi (n\mathrm{integer}).`$
This determines a succession of values that follow the geometric progression

(71)
$$\alpha _ne^{\frac{2\pi n}{\sqrt{7}}}=[10.74\mathrm{}]^n(n+\mathrm{},\mathrm{integer}).$$

At these points a new eigenvalue becomes positive implying secondary instabilities. Clearly, this geometric progression can be traced back to the ‘curious’ asymptotic behavior of the solutions governing the isothermal gas sphere (see Eq. (24)) which intersects the singular solution (21) infinitely often at points that asymptotically increase geometrically in the ratio $`1:10.74`$.

A mode of stability is lost each time that $`\alpha `$ achieves one of the values given by Eq. (71). If one prefers, the same criterion can be expressed in terms of the density contrast $``$. Using Eq. (27) and the asymptotic expansion (24) we get $`=e^{\psi (\alpha )}\frac{\alpha ^2}{2}`$ for $`\alpha +\mathrm{}`$. Hence, the result (71) is translated in

(72)
$$_ne^{\frac{4\pi n}{\sqrt{7}}}=[115.5\mathrm{}]^n(n+\mathrm{},\mathrm{integer}).$$

### 3.4 Profiles of the perturbation at the critical points

In this section, we study the form of the perturbations that trigger instability at the critical points. Following the method of Padmanabhan (1989), it is possible to describe the behavior of these perturbation profiles without numerical work. Indeed, according to Eq. (51), the eigenfunction associated with the eigenvalue $`\lambda =0`$ can be written

(73)
$$\frac{\delta \rho }{\rho _0}=\frac{1}{4\pi \xi ^2}\frac{dF}{d\xi }$$

where $`F(\xi )`$ is given by Eq. (65). Substituting this result in Eq. (73) and simplifying the derivative with the aid of Eq. (20) we obtain

(74)
$$\frac{\delta \rho }{\rho }=\frac{c_1}{4\pi }(2\xi \psi ^{}).$$

Introducing the Milne variable $`v`$ defined by Eq. (28), we can rewrite the foregoing equation in the form

(75)
$$\frac{\delta \rho }{\rho }=\frac{c_1}{4\pi }(2v).$$

The density perturbation $`\delta \rho `$ becomes zero at the point(s) $`\xi _i`$ such that $`v(\xi _i)=2`$. The number of zeros is therefore given by the number of intersections between the spiral in the $`(u,v)`$ plane and the line $`v=2`$ (see Fig. 8). For the first mode of instability $`(\alpha _1=8.99)`$, we need to follow the spiral until the first extremum of $`v`$, which corresponds to the point at which the box terminates (recall that $`u(\alpha _1)=1`$, see section 3.2). In that case there is only one intersection $`\xi _1`$ with the line $`v=2`$. Therefore, the first mode of instability does not present a “core-halo” structure (Fig. 9), contrary to what happens in the microcanonical ensemble (Padmanabhan 1989). For the second mode of instability ($`\alpha _2=98.0`$), we need to follow the spiral until the second extremum of $`v`$. The intersection with the line $`v=2`$ now determines two zeros $`\xi _1`$ and $`\xi _2`$, the first one being the same as for the first mode of instability. The profile of the second mode of instability therefore has a “core-halo” structure (Fig. 10). When we spiral inward, more and more intersections are obtained. Therefore the high order modes of instability present more and more oscillations. For these modes, it is easy to determine the asymptotic positions of their zeros. From Eqs. (28) and (24), we have for $`\xi +\mathrm{}`$:

(76) $`v=2+{\displaystyle \frac{A}{2\xi ^{1/2}}}\times \left[\sqrt{7}\mathrm{sin}\left({\displaystyle \frac{\sqrt{7}}{2}}\mathrm{ln}\xi +\delta \right)+\mathrm{cos}\left({\displaystyle \frac{\sqrt{7}}{2}}\mathrm{ln}\xi +\delta \right)\right].`$
Thus, substituting $`v=2`$ in Eq. (76), we get

(77) $`{\displaystyle \frac{\sqrt{7}}{2}}\mathrm{ln}\xi _i+\delta =\mathrm{arctan}(1/\sqrt{7})+i\pi (i\mathrm{integer}).`$
Therefore, the zeros follow asymptotically the geometric progression <sup>2</sup><sup>2</sup>2It might be recalled, parenthetically, that the positions of the planets in the solar system also follow a geometric progression, but with a different ratio $`q2`$ (Titius-Bode law). This is recognized to be an effect of scale invariance (Graner & Dubrulle 1994, Dubrulle & Graner 1994). There is, however, no real connexion with the present work and planets may have formed quite differently (see Chavanis (2000) and references therein).

(78)
$$\xi _ie^{\frac{2\pi i}{\sqrt{7}}}=[10.74\mathrm{}]^i(i+\mathrm{},i\mathrm{integer}).$$

Note that the perturbation that destabilizes the singular sphere (at the end of the spiral) has an infinite number of oscillations!

## 4 Dynamical stability of isothermal gaseous spheres

We shall now investigate the dynamical stability of isothermal spheres described by Navier-Stokes equations and compare the results with the thermodynamical approach. This problem was previously addressed by Yabushita (1968) in the isobaric ensemble and by Semelin et al. (2001) in the canonical ensemble. We shall prove analytically and with no approximation the equivalence between dynamical and thermodynamical stability and the fact that the stability of isothermal spheres does not depend on the viscosity. We shall also describe the behavior of the velocity profiles at the critical points.

The equations of the problem are the equation of continuity, the equation of motion and the Poisson equation

(79)
$$\frac{\rho }{t}+(\rho 𝐯)=0,$$

(80) $`{\displaystyle \frac{𝐯}{t}}+(𝐯)𝐯={\displaystyle \frac{1}{\rho }}p\mathrm{\Phi }+{\displaystyle \frac{\eta }{\rho }}\mathrm{\Delta }𝐯+{\displaystyle \frac{1}{\rho }}(\zeta +{\displaystyle \frac{\eta }{3}})(𝐯),`$
(81)
$$\mathrm{\Delta }\mathrm{\Phi }=4\pi G\rho .$$

They must be completed by an equation of state that we take of the form $`p=\rho \frac{k}{m}T`$ where $`T`$ is constant. We also recall that we work within a box of radius $`R`$. For astrophysical fluids, the Reynolds numbers are so huge that we can neglect the viscosity in the Navier-Stokes equation. We shall take therefore $`\eta =\zeta =0`$ in most of the calculations but we will also prove that viscosity does not change the onset of instability.

Clearly, the stationary solutions of Eqs. (79)(80)(81) correspond to the isothermal gaseous spheres studied in the previous sections. It should be recalled that the condition of mechanical equilibrium corresponds to the balance between the pressure force and the gravitational force. For an isothermal gas, this yields

(82)
$$T\frac{d\rho }{dr}+\rho \frac{d\mathrm{\Phi }}{dr}=0.$$

We now consider a small perturbation around a stationary solution and write

(83)
$$𝐯=\delta 𝐯(𝐫,t),\rho =\overline{\rho }+\delta \rho (𝐫,t),$$

(84)
$$p=\overline{p}+\delta p(𝐫,t),\mathrm{\Phi }=\overline{\mathrm{\Phi }}+\delta \mathrm{\Phi }(𝐫,t),$$

where the bar refers to the stationary solution (in the following we shall drop the bar). The linearized equations for the perturbations are

(85)
$$\rho \frac{\delta 𝐯}{t}=\frac{kT}{m}\delta \rho \delta \rho \mathrm{\Phi }\rho \delta \mathrm{\Phi },$$

(86)
$$\frac{\delta \rho }{t}+(\rho \delta 𝐯)=0,$$

(87)
$$\mathrm{\Delta }\delta \mathrm{\Phi }=4\pi G\delta \rho .$$

Since we work with exact stationary solutions of the fluid equations, that are inhomogeneous, we do not have to advocate the Jeans swindle (Binney & Tremaine 1987). In this sense, our approach is more rigorous than the classical treatment of Jeans starting from an infinite and homogeneous medium <sup>3</sup><sup>3</sup>3It should be recalled, however, that the Jeans procedure is justified in a cosmological context if we take into account the expansion of the universe..

In the following, we restrict ourselves to radial perturbations. It is known that non-radial perturbations do not lead to new instabilities. Writing the time dependance of the perturbation in the form $`\delta ve^{\lambda t}`$, $`\delta \rho e^{\lambda t}`$,…, the equations of the problem become

(88)
$$\lambda \rho \delta v=\frac{kT}{m}\frac{d\delta \rho }{dr}\delta \rho \frac{d\mathrm{\Phi }}{dr}\rho \frac{d\delta \mathrm{\Phi }}{dr},$$

(89)
$$\lambda \delta \rho +\frac{1}{r^2}\frac{d}{dr}(\rho r^2\delta v)=0,$$

(90)
$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d\delta \mathrm{\Phi }}{dr}\right)=4\pi G\delta \rho .$$

Introducing the notation (51) and using the Gauss theorem (55) we see that the Poisson equation (90) is automatically satisfied. The continuity equation (89) becomes

(91)
$$\frac{\lambda }{4\pi r^2}\frac{dq}{dr}+\frac{1}{r^2}\frac{d}{dr}(\rho r^2\delta v)=0.$$

This equation is readily integrated. Using the boundary condition $`q(0)=0`$, we get

(92)
$$\delta v=\frac{\lambda }{4\pi \rho r^2}q.$$

Substituting this result back into Eq. (88), we obtain

(93)
$$\frac{\lambda ^2}{4\pi r^2}q=\frac{kT}{m}\frac{d}{dr}\left(\frac{1}{4\pi r^2}\frac{dq}{dr}\right)+\frac{1}{4\pi r^2}\frac{dq}{dr}\frac{d\mathrm{\Phi }}{dr}+\frac{G\rho }{r^2}q.$$

Using the condition of mechanical equilibrium (82), the foregoing equation can be rewritten

(94)
$$\frac{k}{m}\frac{d}{dr}\left(\frac{1}{4\pi \rho r^2}\frac{dq}{dr}\right)+\frac{Gq}{Tr^2}=\frac{\lambda ^2}{4\pi \rho Tr^2}q.$$

This is similar to the eigenvalue equation (58) associated with the second variations of the free energy. In particular they coincide for marginal stability $`\lambda =0`$. Therefore, dynamical and thermodynamical instability occur at the same point in the series of equilibrium. This was not obvious a priori since the Navier-Stokes equations without viscosity conserve the free energy. Taking into account a finite viscosity, we obtain instead of Eq. (94),

(95) $`{\displaystyle \frac{k}{m}}{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{1}{4\pi \rho r^2}}{\displaystyle \frac{dq}{dr}}\right)+{\displaystyle \frac{Gq}{Tr^2}}={\displaystyle \frac{\lambda ^2}{4\pi \rho Tr^2}}q{\displaystyle \frac{1}{\rho T}}\left(\zeta +{\displaystyle \frac{4\eta }{3}}\right){\displaystyle \frac{d}{dr}}\left[{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{\lambda q}{4\pi \rho }}\right)\right].`$
For $`\lambda =0`$ the viscous term cancels out. Therefore, viscosity does not change the onset of instability, nor does it alter the profile of the perturbation that triggers the instability.

The profiles of the density perturbation have already been discussed in section 3.4. The profiles of the velocity perturbation at the critical points are given by Eq. (92). In order to work with dimensionless variables, we introduce the speed of sound

(96)
$$c_s^2=\frac{dp}{d\rho }=\frac{k}{m}T=\frac{1}{\beta }.$$

Note that $`c_s`$ can be written $`c_s=\frac{1}{3}L_K/t_{dyn}`$ where $`t_{dyn}=(4\pi G\rho _0)^{1/2}`$ is a time scale comparable with the dynamical time (Binney & Tremaine 1987). Then

(97)
$$\frac{\delta v}{c_s}=\frac{\lambda ^{}}{4\pi \xi ^2}e^\psi F(\xi ),$$

where we have set $`\lambda ^{}=\lambda t_{dyn}`$. Using Eqs. (65)-(28), we find

(98)
$$\frac{\delta v}{c_s}=\frac{\lambda ^{}}{4\pi }c_1\psi ^{}e^\psi (u1).$$

We assume that we are just at the onset of the instability ($`\lambda ^{}=0^+`$) so that Eq. (98) is applicable with $`\lambda ^{}>0`$ (the velocity profile at marginal stability, i.e. $`\lambda ^{}=0`$, is simply $`\delta v=0`$). Since $`\psi ^{}(0)=0`$, the velocity always vanishes at the center of the sphere $`\delta v(\xi =0)=0`$. The other zero(s) are determined by the condition $`u(\xi _i)=1`$. Since the critical points of instability are also characterized by $`u(\alpha )=1`$, see section 3.2, we deduce that the velocity perturbation always vanishes at the wall: $`\delta v(\xi =\alpha )=0`$. This result is to be expected since the mass does not leave the sphere. The other zeros are determined by the intersections between the spiral in the $`(u,v)`$ plane and the line $`u=1`$ (see Fig. 8). They correspond precisely to the values of $`\alpha `$ at which a new mode of instability occurs, i.e. $`\xi _i=\alpha _i`$. It is therefore straightforward to determine the number of zeros in the velocity profile. For the first instability there are only two zeros at $`\xi =0`$ and $`\xi =\alpha _1`$. For the second instability there are three zeros at $`\xi =0`$, $`\xi =\alpha _1`$ and $`\xi =\alpha _2`$. For the $`n^{th}`$ instability there are $`n+1`$ zeros $`\xi =0`$, $`\xi =\alpha _1`$, $`\xi =\alpha _2`$,…, $`\xi =\alpha _n`$. The exact velocity profiles given by Eq. (98) are displayed on Figs. 11-12 for $`n=1,2`$. They are in agreement with the numerical results of Semelin et al. (2001). For high order modes of instability, the zeros follow asymptotically the geometric progression (71). Returning to dimensional variables, the zeros of the $`n^{th}`$ mode of instability are located at $`r=0`$ and $`r_i=(\alpha _i/\alpha _n)R(i=1,\mathrm{},n)`$.

## 5 Fragmentation and fractal structure of an isothermal self-gravitating gas

We shall now express the previous results in a more physical form and make some speculations about the fragmentation and the fractal structure of an isothermal self-gravitating gas.

### 5.1 The King’s length

Introducing the King’s radius (Binney & Tremaine 1987)

(99)
$$L_K=\left(\frac{9}{4\pi G\beta \rho _0}\right)^{1/2},$$

the parameter $`\alpha `$ defined in section 2.4 can be rewritten

(100)
$$\alpha =3\frac{R}{L_K}.$$

The King’s radius gives a good estimate of the core radius of an isothermal sphere. In the following, we shall consider that the core radius $`L_K`$ is fixed and that the domain size $`R`$ is increased. Clearly, this amounts to increasing the parameter $`\alpha `$ along the spiral. From the study of section 3.2, we find that the onset of instability corresponds to

(101)
$$R>\frac{8.99}{3}L_K.$$

For larger radii there are still critical points of free energy but they are more and more unstable. According to Eq. (71), a new mode of instability occurs at radii that follow asymptotically the geometric progression

(102)
$$R_n[10.74\mathrm{}]^nL_K.$$

This result is clearly related to that of Semelin et al. (1999) who first noticed the appearance of a hierarchy of scales connected to the instability of isothermal gaseous spheres. However, the connection with their work is not straightforward. Semelin et al. (1999) approximate the regular solutions of Eq. (17) by the analytical profile (18) and introduce an arbitrary short cut-off distance $`\delta `$ to avoid the central singularity. By contrast, there is no cut-off at small scales in our theory since we consider regular spheres that are exact solutions of the Emden equation (although we do not have analytical expressions for them). In this sense, there are no approximation nor indetermination in our theory. In fact, an “effective” cut-off is played by the King’s length (or core radius) but it is not arbitrary and depends on physical parameters, the central density $`\rho _0`$ and the temperature $`T`$. For $`rL_K`$, the density profile behaves like $`r^2`$ so that regular gaseous spheres can indeed be approximated by a $`r^2`$ profile plus a cut-off $`\delta `$ at small scales as done by Semelin et al. (1999). Our approach suggests to take $`\delta `$ as the King’s length. We have therefore given a new derivation of their result which is considerably simpler and does not introduce sophisticated renormalization group technics. It also avoids the introduction of an arbitrary cut-off. This provides, hopefully, an easier interpretation of their intriguing results.

According to section 3.4, if we increase the domain size along the series of equilibrium (i.e. maintaining a fixed core radius) smaller and smaller regions become unstable. The modes that trigger these instabilities present numerous oscillations whose zeros also follow a geometric progression. These oscillations can be regarded as sort of ‘germs’ leading to the formation of ‘clumps’ preceding the fragmentation of the isothermal gas. It is expected that these ‘clumps’ will evolve by achieving higher and higher density contrasts, and finally fragment in turn into substructures. This yields a hierarchy of structures fitting one into each other in a self-similar way. This picture is given further support by the fact that both the domain sizes that induce instability and the zeros of the perturbation in a given domain follow a geometric progression with the same ratio (see Eqs. (102)(78)). This property may explain in a quite natural fashion the fractal structure that a self-gravitating medium (e.g., the interstellar medium, the large scale structures of the universe,…) can build under certain conditions.

The hierarchy of scales that appears in an isothermal gas is due intrinsically to the oscillating nature of the equilibrium phase diagram (see Figs. 1,2,3). This curious behavior has been known for a long time in the context of stellar structure (Chandrasekhar 1942) but it had not been apparently related to the geometric progression of scales inducing instability.

### 5.2 The Jeans length

We can also introduce another length scale defined with the average density $`\overline{\rho }=\frac{3M}{4\pi R^3}`$ of the system instead of the central density. This is the Jeans length (Binney & Tremaine 1987)

(103)
$$L_J=\left(\frac{9}{4\pi G\beta \overline{\rho }}\right)^{1/2}.$$

The reduced temperature $`\eta `$ can be expressed in terms of the Jeans length as (de Vega & Sanchez 2001a)

(104)
$$\eta =\frac{\beta GM}{R}=3\left(\frac{R}{L_J}\right)^2.$$

In the following, we shall consider that the Jeans length $`L_J`$ is fixed and that the domain size $`R`$ is increased. Clearly, this amounts to increasing the parameter $`\eta `$ along the vertical axis of Fig. 2. Using the results of section 2.6, we find that the gaseous sphere becomes unstable when

(105)
$$R>\sqrt{\frac{2.52\mathrm{}}{3}}L_J.$$

For larger radii the system will collapse because there are no critical point of free energy (i.e. hydrostatic equilibrium). In that case there are clearly no secondary instabilities. Indeed, when the condition (105) is realized we leave the spiral and loose the geometric hierarchy of scales associated with its turning points. Since the mode of instability at the critical point does not oscillate (see section 3.4), we expect that the system will collapse without fragmenting, probably like in the study of Penston (1969).

### 5.3 Fragmentation or not fragmentation?

The fragmentation of the isothermal self-gravitating gas and the fractal structure that it can generate were previously discussed by de Vega et al. (1996a, 1996b, 1998) and we here improve and complete their physical arguments. In particular, the analytical calculation of the modes of instablity (without approximation), the geometric progression of both the domain sizes and the ‘clumps’ sizes, and the precise physical picture that emerges from these results are new. In addition, we give arguments explaining why the instability of a self-gravitating gas does not always lead to a fractal structure (indeed globular clusters or elliptical galaxies, for example, do not show a hierarchical structure). In our model, the selection of the regime is related to how the size of the domain compares with the Jeans length or the King’s length. If $`R>0.916\mathrm{}L_J`$, we expect a collapse without fragmentation but if $`R>2.996\mathrm{}L_K`$ (with a sufficiently large ratio) the medium is expected to fragment into several ‘clumps’.

The distinction between the Jeans scale and the King’s scale is therefore essential in the present context and had not been noticed previously. In particular it is not the Jeans lengths but the King’s lengths that follow a geometric progression \[see formula (102)\] contrary to the claim of Semelin et al. (1999). These two scales only coincide for a perfect gas with no self-gravity that is homogeneous. Indeed, for $`\alpha 0`$, we have $`L_JL_K\frac{3R}{\alpha }+\mathrm{}`$. By contrast, for $`\alpha +\mathrm{}`$ (corresponding to the singular sphere) we get $`L_K\frac{3R}{\alpha }0`$ and $`L_J\sqrt{\frac{3}{2}}R`$ (we have used the asymptotic expansions of section 2.4). To appreciate the difference between the two length scales, it is convenient to introduce a dimensionless number

(106)
$$𝒩=\frac{L_J}{L_K}.$$

From the above formulae, we have the equivalent expressions

(107)
$$𝒩=\left(\frac{\rho _0}{\overline{\rho }}\right)^{1/2}=\frac{\alpha }{\sqrt{3\eta }}=\left(\frac{\alpha }{3\psi ^{}(\alpha )}\right)^{1/2}.$$

For $`\alpha 0`$ (perfect gas), $`𝒩1`$ and for $`\alpha +\mathrm{}`$ (singular sphere), $`𝒩\frac{\alpha }{\sqrt{6}}+\mathrm{}`$. Therefore the number $`𝒩`$ quantifies the inhomogeneity of the self-gravitating gas. The evolution of $`𝒩(\alpha )`$ along the spiral is ploted on Fig. 13. From Eq. (107), $`𝒩`$ is just the ratio between the central density and the average density to the power $`1/2`$.

## 6 Conclusion

In the first part of the paper, we have investigated the stability of isothermal spheres in the canonical ensemble using a thermodynamical approach. This is a natural extension of the work of Padmanabhan (1989) in the microcanonical ensemble. As is well known, there is no hydrostatic equilibrium for an isothermal gas below a critical temperature $`kT_c=\frac{GmM}{2.52R}`$ (Lynden-Bell & Wood 1968). This corresponds to a situation in which the domain size becomes comparable with the Jeans length, i.e. $`R>0.916\mathrm{}L_J`$. In that case, we expect the system to collapse. By analyzing the second variations of free energy, we have found that the point of minimum temperature is precisely the point at which the free energy ceases to be a local maximum and becomes a saddle point in the series of equilibrium. This vindicates the theorem of Katz (1978) and provides in addition the form of the perturbation that induces instability. Contrary to the microcanonical ensemble, it has not a “core-halo” structure, nor multiple oscillations. We expect therefore the system to collapse without fragmenting. On the other hand, if we start from isothermal gaseous spheres with high density contrasts $`>32.2`$ (this corresponds to $`R>2.996\mathrm{}L_K`$), such spheres are unstable saddle points and probably fragment into small ‘clumps’ (although they would be said to be stable by a naive Jeans argument since $`R<0.916\mathrm{}L_J`$). The number of clumps, presumably associated with the number of zeros in the mode of instability, increases with the density contrast or, equivalently, as the domain size becomes larger and larger with respect to the King’s radius. The distinction between the Jeans length and the King’s length is essential in the present context and had not been noticed previously. The Jeans length is related to the parameter $`\eta `$ determining whether an hydrostatic equilibrium state exists or not (for an isothermal sphere) while the King’s length is related to the parameter $`\alpha `$ parametrizing the series of equilibrium. They are very different parameters.

In the second part of the paper, we have extended Jeans instability criterion to the case of inhomogeneous bounded gaseous spheres described by Navier-Stokes equations. This treatment avoids the well-known difficulties associated with an infinite and homogeneous medium (Jeans swindle). Quite remarkably, we have provided an elegant analytical solution of this more general problem without making any approximation.

Now, the question that naturally emerges is what happens to a gravitationally unstable gaseous configuration. The collapse of isothermal gas spheres has been investigated by Penston (1969) using Navier-Stokes equations (without viscosity). He found an analytical solution for which the collapse is self-similar and develops a finite time singularity (i.e., the central density becomes infinite in a finite time). This only partially answers the above question. Indeed, his solution is just a particular solution of the fluid equations and other solutions (e.g., fragmentation) may be possible. In particular, it is not clear whether his solution is stable with respect to non spherically symmetric perturbations. In addition, Penston (1969) works in an infinite domain for which the isothermal spheres have infinite mass and are always unstable. Therefore, he could not evidence a phase transition between an equilibrium “gaseous” state and a “collapsed” state, depending on the value of the temperature.

The thermodynamics of self-gravitating systems is so intriguing that we have been tempted to explore this problem anew from a dynamical point of view (Chavanis et al. 2001). We have introduced a numerical algorithm in the form of a relaxation equation (Chavanis 1996, 2001a; Chavanis et al. 1996) constructed so as to increase monotonically entropy at fixed energy in the microcanonical ensemble and free energy at fixed temperature in the canonical ensemble. With this numerical algorithm, we have covered the whole bifurcation diagram in parameter space and we have checked, by an independant method, the stability limits of Katz (1978). We have verified that the density profile that induces instability at the critical point has a “core-halo” structure in the microcanonical ensemble but not in the canonical ensemble, in agreement with the study of Padmanabhan (1989) and the present study. We have also checked that the number of oscillations in the perturbation profile increases when we start from unstable isothermal spheres with high density contrasts. With this algorithm, we have explored the “bassin of attraction” of the stable isothermal spheres. Since they only are local maxima of entropy or free energy, an initial condition can either relax towards a metastable equilibrium state or collapse, depending on its topology. Hence, the control parameter $`\mathrm{\Lambda }`$ or $`\eta `$ is not sufficient to characterize completely the final state of the system even in the case when an equilibrium exist. This depends whether the initial condition lies in the “bassin of attraction” of the equilibrium solution or not. The complete characterization of this bassin of attraction is difficult but we have given some examples. With our numerical algorithm, it will be possible to include rotation in order to compute entropy maxima that are not spherically symmetric and determine the influence of angular momentum conservation on the collapse of a rotating self-gravitating system (Chavanis & Rosier, in preparation). Some interesting works have been done in that direction (Lagoute & Longaretti 1996, Laliena 1999, Fliegans & Gross 2001 and Lynden-Bell 2000) but the complete description of rotating isothermal gaseous spheres is still in its infancy stage.

It is also of first interest to test numerically the speculations that we have made concerning the development of the instability, i.e. if the system will collapse as a whole or if it will fragment and break into substructures. These two regimes are clearly separated, at least in the theoretical model that we have considered here (Antonov model). We have suggested that if we start from a solution of the Emden equation with a high density contrast, it should break in a series of ‘clumps’ associated with the oscillations of the density perturbation that we have calculated in the linear regime. Obviously, if we want to follow the nonlinear evolution correctly, we must relax the hypothesis of spherical symmetry.

In this paper, we have assumed that the particles can be treated by classical mechanics and that they are point-like. The extension of our work to general relativity (for, e.g., neutron stars) is almost straightforward (Chavanis 2001b). We find that there is no hydrostatic equilibrium if the system size is smaller than a multiple of the Schwarzschild radius. We also show analytically that the point of smallest radius corresponds precisely to the point at which the series of equilibrium becomes unstable. Secondary modes of instability are found; they follow a geometric progression like in the Newtonian case but with a different ratio. When the speed of light goes to $`+\mathrm{}`$, we recover the results of the present study. It is seen that relativistic effects favour gravitational instability, i.e. instability occurs sooner than in the Newtonian limit. On the other hand, the case of self-gravitating fermions enclosed within a box has been considered by Chavanis & Sommeria (1998). Their study was formulated in the context of the “violent relaxation” of collisionless stellar systems introduced by Lynden-Bell (1967) but the results also apply to quantum particles such as massive neutrinos in Dark Matter models (Ingrosso et al. 1992) or degenerate electrons in White Dwarf stars (Chandrasekhar 1959). When degeneracy is accounted for, there exists a global entropy maximum for all values of energy (Chavanis & Sommeria 1998, Robert 1998). It has a “core-halo” structure with a degenerate core and a dilute atmosphere. Depending on the degeneracy parameter, there can be a “gravothermal catastrophe” at $`E=0.335GM^2/R`$ but the core ceases to shrink when it becomes degenerate. Considering a classical gas with a short distance cut-off $`a`$ essentially leads to the same results (see, e.g., Padmanabhan 1990, Follana & Laliena 2000).

For astrophysical purposes, it is still a matter of debate to decide whether collisionless stellar systems like elliptical galaxies are degenerate (in the sense of Lynden-Bell) or not. Since degeneracy can stabilize the system without changing its overall structure at large radii, we have suggested that degeneracy could play a role in galactic nuclei (Chavanis & Sommeria 1998). The recent simulations of Leeuwin and Athanassoula (2000) and the theoretical model of Stiavelli (1998) seem to go in that direction: if the nucleus of elliptical galaxies contains a (primordial) black hole, degeneracy must be taken into account and can explain the cusps observed in the center of galaxies. This form of degeneracy is also relevant for massive neutrinos in Dark Matter models where it competes with quantum degeneracy (Kull et al. 1996).

All together, these results suggest that the statistical mechanics and the thermodynamics of self-gravitating systems can have quite relevant astrophysical applications at different scales of the universe: stars (Chandraskhar 1959), interstellar medium (de Vega et al. 1996a, 1996b), globular clusters (Lynden-Bell & Wood 1968), elliptical galaxies (Lynden-Bell 1967, Hjorth & Madsen 1993, Chavanis & Sommeria 1998), dark matter (Ingrosso et al. 1992, Kull et al. 1996), cosmology (Saslaw & Hamilton 1984, de Vega et al. 1998)…The same ideas of statistical mechanics have been introduced in two-dimensional turbulence to explain the formation and maintenance of coherent vortices, like Jupiter’s Great Red Spot, which are common features of large-scale geophysical and astrophysical flows (see Bouchet & Sommeria 2000, Bouchet et al. 2001 and references therein). The analogy between the statistical mechanics of two-dimensional vortices and stellar systems has been described extensively by Chavanis (1996, 1998b, 2001c). This analogy concerns not only the equilibrium state (the formation of large-scale structures) but also the relaxation towards equilibrium (Chavanis et al. 1996, Chavanis 2001d) and the statistics of fluctuations (Chavanis & Sire 2000). These ideas may also have applications in the context of planet formation where large-scale vortices present in the solar nebula could efficiently trap dust particles to form the planetesimals and the planets (Barge & Sommeria 1995, Tanga et al. 1996, Bracco et al. 1999, Chavanis 2000). Accordingly, the statistical mechanics of these nonlinear media seems to be able to account for the fascinating process of self-organization in nature.

## 7 Acknowledgments

I acknowledge interesting discussions with L. Athanassoula, F. Combes, O. Fliegans, J. Hjorth, H. Kandrup, V. Laliena, F. Leeuwin, D. Lynden-Bell, T. Padmanabhan, R. Robert, B. Semelin and J. Sommeria on these exciting problems of thermodynamics. This work was initiated during my stay at the Institute for Theoretical Physics, Santa Barbara, during the program on Hydrodynamical and Astrophysical Turbulence (February-June 2000). This research was supported in part by the National Science Foundation under Grant No. PHY94-07194.


format: text

text:
# Nature of widely separated ultraluminous infrared galaxies

## 1 Introduction

Ultraluminous infrared galaxies (ULIRGs) are known to emit the bulk of their energy in the infrared and their infrared luminosities ($`L_{IR}10^{12}L_{}`$) are equivalent to the bolometric luminosities of optically selected quasars. For the past decade, there have been extensive studies of the nature of these ULIRGs (Sanders & Mirabel 1996). Imaging studies at optical and near-infrared wavelengths show that most of these galaxies are either interacting or merging (Sanders et al. 1988, Zou et al. 1991, Leech et al. 1994, Murphy et al. 1996, Borne et al. 2000, Kim et al. 2001). The fraction of ULIRGs in interacting or merging ranges from at least 61% (25 objects out of 41, Zou et al. 1991) to possibly as high as 100% for all 10 nearby ULIRGs (Sanders et al. 1988). The rather small fraction of interacting systems derived by Zou et al. (1991) resulted from the use of POSS images in which faint tidal tails are almost impossible to detect. Recent imaging data by Kim et al. (2001) showed that close to 100% of ULIRGs in the complete 1 Jy sample have faint tidal tails, a signature of recent interaction and merger.

One of the most basic issues of ULIRGs is the energy source powering the enormous amount of their infrared radiation. Spectroscopic studies in the optical, near-infrared, and mid-infrared suggest that about 70% and 30% of these galaxies are powered by starburst activities and active galactic nuclei (AGNs) respectively (Veilleux et al. 1999a, Veilleux et al. 1999b, Genzel et al. 1998a). The huge infrared luminosity of the ULIRGs suggests that they contain a large reservoir of molecular gas to provide fuel for the star formation activity and/or a central AGN. Extensive searches for CO emission have indeed shown that all ULIRGs possess a huge amount of molecular gas ($``$ 1010 M) and that there is a trend of increasing molecular gas content with the far-infrared emission (Downes et al. 1993, Solomon et al. 1997). Recent interferometric observations of nearby ULIRGs (Downes & Solomon 1998) and luminous infrared galaxies (Bryant & Scoville 1999) have determined that the molecular gas is concentrated in a circumnuclear disk of typically 1 kpc in size. From the interferometric observations of double-nucleus ULIRGs with projected nuclear separations of 3-5 kpc, Evans et al. (1999, 2000) find that one of the interacting pair contains the majority of molecular gas in warm ULIRGs whereas both components contain molecular gas in cool ULIRGs. Based on observational evidence, Sanders et al. (1988) proposed that ULIRGs originate from the merger of two gas-rich spiral galaxies. The large amount of molecular gas accumulated in the center of the galaxies due to the gravitational interaction can sustain strong star formation activities as well as provide fuel to the central AGN. As the merger proceeds, the molecular gas reservoir is eventually depleted by the star formation (Gao & Solomon 1999) and dispersed by the AGN activities and supernovae. Finally the AGN reveals itself as an optical quasar. The estimated lifetime of the ULIRG phase is about 107 to 108 yr. Genzel et al. (1998b) reviewed most recent observational evidence in favor of this scenario.

Numerical simulations of interacting galaxies have shown that violent star formation activity can be triggered in the merger of two massive disk galaxies (Barnes & Hernquist 1996, Mihos & Hernquist 1996, Springel 2000). Furthermore, star formation activities are found to be dependent mainly on the galaxy structure and the orbital geometry plays only a modest role (Mihos & Hernquist 1994, 1996). They suggest that galaxies with dense central bulges show strong star formation activities at the final stages of merging, while bulgeless galaxies show weaker star formation activity at the early stage of interacting and do not show any subsequent strong star formation activities. Thus, according to the numerical simulations, the majority of ULIRGs would be in the final stages of merger.

Recently, Kim et al. (2001) have completed an optical and near-infrared imaging survey of the IRAS 1 Jy sample of ULIRGs (Kim & Sanders 1998) and find that all of the ULIRGs show signs of tidal interaction and therefore are in a merging or interacting process (Kim et al. 2001). More interestingly, we also find that 6 out of 118 ULIRGs have projected nuclear separation larger than 20 kpc. The large nuclear separation observed in these ULIRGs seems to contradict the conventional view of the origin of ULIRGs and the numerical simulation results outlined above. To become an ULIRG in the early phase of interaction, the galaxies must have accumulated a large amount of molecular gas and it is still unclear how the star formation activities are triggered in these cases.
The position uncertainty associated with the poor resolution of IRAS in the far-infrared did not allow us to determine the distribution of far-infrared emission as well as the molecular gas in these galaxies. Without that information, it is difficult to clarify the nature of these unusual ULIRGs. In this paper we present high angular resolution CO(J=1-0) observations and new K’-band imaging to study the widely separated ULIRGs. CO emission and the molecular gas are known to be associated with the source of far-infrared emission. Thus observations of CO emission will allow us to determine the location and the distribution of far-infrared luminosity of these ULIRGs, in addition to the distribution of molecular gas. Also making use of the optical imaging data (Kim et al. 2001) and our new K’-band images, we discuss the nature of these widely separated ULIRGs. We will use $`q_0=0.5`$ and $`H_0=75`$ km s-1 Mpc-1 throughout this paper.

## 2 Observations

We used the BIMA (Berkeley-Illinois-Maryland Association) array (Welch et al. 1996) to search for the CO(J=1-0) emission line from the selected ULIRGs. Each galaxy in our sample was observed in various runs during 1998 and 1999 (Table. 1) in the B-configuration with baselines ranging from 5 k$`\lambda `$ to 80 k$`\lambda `$, and the C-configuration with baselines ranging from 3 k$`\lambda `$ to 30 k$`\lambda `$. The digital correlator was configured to cover a velocity range of 1600 km s-1 centered at the velocity of the ULIRGs calculated from their redshift. Phase calibration was achieved by observing nearby quasars every 20 to 30 minutes. We also observed planets to determine the absolute flux scale.

We used the MIRIAD data reduction package (Sault et al. 1994) to calibrate and process our data. After phase calibration and flux calibration the visibilities were Fourier-transformed to make the images. These images were then cleaned using the standard task CLEAN of the MIRIAD package. The final maps have angular resolutions of about 4 arcsec for B-configuration and 8 arcsec for C-configuration, respectively.

The K’-band (2.1 $`\mu m`$) images of 4 galaxies IRAS F11180+1623, IRAS F14394+5332, IRAS F17028+5817, IRAS F23327+2913 were obtained with the QUIRC 1024 x 1024 near-infrared camera on July 15, 2000, at the UH 88” telescope on Mauna Kea. The plate scale of the QUIRC is 0.1886′′/pixel with f/10 and the seeing is estimated about 0.5′′ to 0.7′′. The sky frames were constructed by median filtering of 9 dithered (20′′ rectilinear motion) object frames. Dome flats were obtained by turning a lamp in the dome on and off at the beginning and end of each night. The dark subtracted and normalized flats were used for flat-fielding the sky-subtracted images. Flat-fielded images were then coadded after shifting in fractional pixel units with respect to the object center. The final coadded images were calibrated from observations of infrared standard stars (Elias et. al. 1982). K’-band images of two other galaxies (IRAS F10485-1447, IRAS F12359-0725) together with R-band images of all the galaxies in our sample were taken from the work of Kim et al. (2001). The optical spectra of the nothern object of IRAS F23327$`+`$2913 were taken in October 1995 at the UH 88” telescope with spectral resolution of 8 $`\AA `$ in the spectral range of $`\lambda \lambda `$ 4500-8500$`\AA `$. The spectra were processed with the standard data reduction procedures.

## 3 Results

We detected CO(J=1-0) emission from 4 galaxies IRAS F11180$`+`$1623, IRAS F14394$`+`$5332, IRAS F17028$`+`$5817, IRAS F23327$`+`$2913 at $``$ 3$`\sigma `$ in the channel maps. In IRAS F12359-0725 we found possible CO(J=1-0) emission of 3$`\sigma `$ with channel width of 200 km s-1 which is coincident within the synthesized beam with one of the galaxy in the system. Fig 1. shows the integrated intensity maps for the galaxies observed in our program together with their R-band and K’-band images. The CO emission line profiles of the detected nuclei are shown in Fig. 2. Some of the properties of these ULIRGs are provided in Table. 2. The CO emitting regions are not resolved by our synthesized beams. In the case of IRAS F14394$`+`$5332, which is the only ULIRG mapped in B-array, we can infer an upper limit to the size of the CO emitting region to be $``$ 1” or $``$ 2 kpc in linear scale.

To estimate the molecular gas content from the total flux of CO(J=1-0) line we adopt the same conversion factor as for the Milky Way galaxy (Bryant & Scoville 1999):

$$\mathrm{M}_\mathrm{g}=1.20\times 10^4\mathrm{F}_{\mathrm{CO}}\mathrm{D}_\mathrm{L}^2(1+z)^1$$
(1)
where Mg is the molecular gas mass in M, FCO is the integrated flux in Jy km s-1, DL is the luminosity distance in Mpc, and $`z`$ is the redshift. The molecular gas mass derived for the galaxies in our sample is similar to that found in other ultraluminous galaxies (Solomon et al. 1997) and in some high redshift galaxies (Combes et al. 1999) using the same Galactic conversion factor as ours. The 1$`\sigma `$ upper limit of the CO flux is rms\*$`\mathrm{\Delta }V_{\mathrm{FWHM}}`$/$`\sqrt{\mathrm{\Delta }V_{\mathrm{FWHM}}/\mathrm{\Delta }V_{\mathrm{res}}}`$ where rms is the noise level in the channel maps, $`\mathrm{\Delta }V_{\mathrm{FWHM}}`$ is the linewidth determined from the gaussian fitting to the line profiles and $`\mathrm{\Delta }V_{\mathrm{res}}`$ is the channel width. The 1$`\sigma `$ upper limit to the molecular gas mass of the companion is then derived following the above equation. For the ULIRGs with detected CO emission, the companion may possess a molecular gas content of up to 10% to 20% the molecular gas content of the main component (Table. 2). The upper limit is in the range of several 109 M which implies the companions are not as gas rich as the main component but may possess a molecular gas content comparable to that of normal spiral galaxies.

Gaussian fitting to the line profiles gives the FWHM of 200 to 400 km s-1 except for the case of IRAS F11180$`+`$1623. The weak CO(J=1-0) emission from IRAS F11180$`+`$1623 did not allow us to determine accurately its linewidth. Thus, the value of 600 km s-1 should be regarded as upper limit to the true linewidth of IRAS F11180$`+`$1623.

Although we did not succeed to detect CO(J=1-0) emission from IRAS F10485$``$1447, the 1$`\sigma `$ upper limit to the molecular gas mass of $``$ 4$`\times `$109 M is still quite high. In addition, this galaxy has been observed only with the B-array which provides long baselines data with high phase scattering due to the atmospheric fluctuation. It would be very useful to search again for CO emission with improved sensitivity.
Using optical spectra, Veilleux et al. (1999a) determined the spectral type of all ULIRGs in the 1 Jy sample. In all cases, the galaxy detected in CO(J=1-0) line has LINER or Seyfert II spectral type while the companion has no emission line except for IRAS F17028$`+`$5817E galaxy, which has spectral type HII, typical of starburst galaxies (Table 3.).
The magnitudes of each component of the ULIRGs are derived from our K’-band images and presented in Table. 3. In most cases the galaxies have compact light distribution in K’-band (Fig. 1). However, K’-band image of IRAS F14394$`+`$5332 shows the extended emission with distortion in the E component, suggestive of a recent merger.

## 4 Detailed Description of Individual System

IRAS F10485$``$1447 \- This system is close to two bright foreground stars which are visible on the upper left of both R-band and K’-band images (Fig. 1). On the R-band image, the E component shows short tidal tails in two opposite directions and the W component has a long curved tidal tail of about 10 kpc. The W component also has a distortion toward SW direction. The R-band luminosity of the W component is only slightly brighter than that of the E component, whereas in the K’-band the W component is more than one magnitude brighter than the E component. The optical emission lines were detected from the W component with LINER spectral type (Veilleux et al. 1999a). The presence of absorption lines in the spectrum of the E galaxy (S. Veilleux, private communication) indicates that both galaxies are at the same redshift and are interacting. Though we did not detect any CO emission in this system, the general trend of finding molecular gas in the active component of the widely separated ULIRGs suggests that most of the gas could reside in the W component.
IRAS F11180$`+`$1623 \- The optical morphology of the E component (Fig. 1) shows two tidal tails in two opposite directions whereas the W component does not show any strong disturbances. The spectral type of the E component is LINER (Veilleux et al. 1999a) and no detectable optical emission lines were observed in the W component. The CO emission was detected only on the E component. The E component may itself be a merger because it has two opposite tidal tails. One of the tidal tails on the E component tends to be curved and directed to the W component. There is also a small distortion on the east side of the W component, suggesting that the two galaxies are interacting. The lack of CO emission as well as the non-detection of optical emission lines in the W component imply that this galaxy contains only a small amount of gas and is in the early stage of tidal interaction with the E component. Because there is no measured redshift for the W component, we classify this pair as a tentative merger system. There are 3 additional galaxies around the system, and if they are at the same redshift as the two galaxies considered, then this whole system is likely an example of the Hickson compact groups (Hickson 1982).
IRAS F12359$``$0725 \- On the R-band image (Fig. 1), the S component shows a curved tidal tail whereas the N component has no visible tail except a small depression on the western part of the galaxy. The two galaxies are found to be at the same redshift and are likely in the early phase of interaction. A foreground star is present to the west of the system on both R-band and K’-band images. The R-band magnitudes of both components are similar, whereas K’-band magnitude of the N component is 1.1 magnitude brighter than that of the S component. The N galaxy with LINER spectral type is more active than the S component which has HII spectral type. We found possible CO(J=1-0) emission coincident with the N component within the synthesized beam of 10”$`\times `$6”. As the N component is the most active in the system, the detection follows the same trend for other widely separated ULIRGs in our sample.
IRAS F14394$`+`$5332 \- This system has large nuclear separation and shows long tidal tail ($``$ 35 kpc, Fig. 1). The R-band luminosity of the E component is slightly higher than that of the W component. On K’-band, the E component is about 1.1 magnitudes brighter than the W component. Furthermore, the E component shows strong optical emission lines and has a Seyfert II spectral type. In contrast, the W component shows no detectable emission lines. We detected CO(J=1-0) emission from the E component only.

The E component has Seyfert II spectra and bright near-infrared nucleus. Radio continuum observations at 20 cm with the VLA - B array (4 arcsec resolution) detected a fairly strong continuum source of 38 mJy (Crawford et al. 1996) at the same position of this galaxy. This result is expected as Seyfert galaxies are known to be strong continuum sources. An obscured AGN at the center of the E component is likely the source of most of the continuum emission.

Our new K’-band image of IRAS F14394$`+`$5332 (Fig. 1) under good seeing condition (0.5′′) reveals that the E component has two nuclei in the center with nuclear separation of 1.29′′. Thus, K’-band morphology suggests that this galaxy is already a merging system.
IRAS F17028$`+`$5817 \- On the R-band image (Fig. 1), the W component is larger than the E component, but it has a more diffuse light distribution. Faint emission resembling tidal tails can be seen in the W component, in the R-band images of Kim et al. (2001) and Murphy et al. (1996). On K’-band image, the W component is much brighter than the E component and shows a compact light distribution. The CO emission was detected in this W component. The W component also has a very red nuclear color as observed in the typical merger ULIRGs (Kim et al. 2001). Strong optical emission lines detected from both components indicate that they are at the same redshift. The spectral types of the W & E components are LINER and HII region, respectively.
IRAS F23327$`+`$2913 \- On the R-band image (Fig. 1), the S component resembles a normal spiral with very thick bar structure, whereas the N component looks much more disturbed with a big plume on the northern direction and has small protrusion on the eastern direction. Both components have comparable R-band magnitudes. In K’-band the S component shows a compact distribution and is brighter than the N component. Interestingly, the optical emission lines were detected only from the S component and not from the much disturbed N component, except for the strong stellar absorption line of Mg Ib $`\lambda `$5175Å (EW (Mg Ib) = 5.5Å, Fig. 3). The spectral type of the S component is LINER and the CO emission line was detected only from this component. Radio continuum emission at 20 cm and 6 cm has been detected in this component (Crawford et al. 1996). The continuum emission at 6 cm was resolved by the 4 arcsec beam of the VLA.
The large equivalent width of the Mg Ib absorption line which originates from the intermediate stellar population ($`10^9`$ yr, Bica, Alloin, & Schmidt 1990), and the lack of Balmer absorption lines (originating from stellar population of age $`10^810^9`$ yr) both suggests that star formation activity in this galaxy has already ceased some time ago. It is unclear whether the big northern plume on the N component is a remnant of past tidal interaction or caused by recent tidal interaction with the S component. The latter is more likely, since if the plume was caused by the past tidal interaction more than 109 years ago then, it should already have been relaxed at this time and not be observable.

## 5 Discussion

### 5.1 Molecular Gas Mass

The most interesting result of our observations is the detection of CO emission from only one component in the pair of the selected ULIRGs. The molecular gas content in this component is quite high (several times 1010 M) if we adopt the Milky Way conversion factor. However, the value of the conversion factor has been suggested to be lower for ULIRGs (Solomon et al. 1997). By modelling the high resolution observations of CO lines from nearby ULIRGs, Downes & Solomon (1998) derived molecular gas masses which are a factor of 5 lower.

From the IRAS fluxes at 60 $`\mu m`$ and 100 $`\mu m`$ we have estimated the dust temperatures in the ULIRGs detected with CO emission. The dust temperatures are roughly the same in these galaxies, $``$ 40 K (Table. 4). If we assume that the gas is in thermal equilibrium with the dust, namely Tgas = Tdust, we can estimate the radius of the CO emitting region from the observed CO flux FCO (Solomon et al. 1997):

$`L_{\mathrm{CO}}^{^{}}`$ $`=`$ $`3.25\times 10^7F_{\mathrm{CO}}\nu _{\mathrm{obs}}^2D_L^2(1+z)^3`$ (2)
$`R_{\mathrm{CO}}`$ $`=`$ $`(L_{\mathrm{CO}}^{^{}}/\pi T_{\mathrm{gas}}f_\mathrm{V}\mathrm{\Delta }V_{\mathrm{FWHM}})^{0.5}`$ (3)
where fV is the velocity filling factor, $`\mathrm{\Delta }`$VFWHM is the linewidth and $`\nu _{\mathrm{obs}}`$ is the observed frequency in GHz, DL is the luminosity distance in Mpc. For simplicity, we assume that the velocity filling factor fV is equal to unity. We also assume the thermalization and optical thickness of CO(J=1-0) line, with the CO(J=1-0) brightness temperature equal to the gas temperature. The dynamical (virial) mass (Table. 4) contained inside the CO emitting region is R$`{}_{\mathrm{CO}}{}^{}\mathrm{\Delta }V_{\mathrm{FWHM}}^2`$/G, where G is the gravitational constant. We find that the dynamical masses are a factor of 3 to 4 lower than the molecular gas masses. As a result, the H2-to-CO luminosity conversion factor for ULIRGs must be lower than the Galactic value. It was suggested by Downes et al. (1993) and Solomon et al. (1997) that in ULIRGs the dynamical mass is dominated by the molecular gas mass. In this case, the dynamical masses estimated above are the upper limit to the actual molecular gas mass.

### 5.2 Far Infrared and Radio Continuum Emission

The component containing most of the molecular gas is also the most active galaxy in the pair, having either Seyfert II or LINER-like spectral type. The other component is less active having either a HII spectral type typical of starburst galaxies or no detectable emission lines. Since there is a correlation between CO emission and the 100 $`\mu m`$ flux from luminous and ultraluminous galaxies (Solomon et al. 1997), the galaxy detected in CO is likely to be the source of most of the far-infrared luminosity coming from the system and presumably the site of vigorous star formation activity. Our result is consistent with the Evans et al. (1999, 2000) who found most of the molecular gas in several warm double-nucleus ULIRGs to be associated with only the Seyfert nuclei. All ULIRGs in our sample are also found to be associated with 20 cm continuum sources (Becker et al. 1994, Condon et al. 1997). The continuum sources have flux densities ranging from a few mJy to $``$ 40 mJy. For two galaxies IRAS F14394$`+`$5332 and IRAS F23327$`+`$2913, high angular continuum observations have been reported (Becker et al. 1994, Crawford et al. 1996). The continuum sources in these two galaxies coincide with the position of CO(J=1-0) emission. In particular, the spectral index of radio continuum emission derived from 20 cm and 6 cm observations of IRAS F23327$`+`$2913 is -0.52 (Crawford et al. 1996), suggesting a non-thermal nature of the emission. For other ULIRGs in our sample, due to the lack of observations at other frequencies, we could not determine the nature of the continuum emission. However, IRAS F14394$`+`$5332E, a Seyfert II galaxy, appears to have stronger radio continuum flux relative to CO luminosity, compared to other ULIRGs in this sample. Significant fraction of the radio continuum could arise from an AGN in this galaxy.

### 5.3 What is the origin of the widely separated ULIRGs ?

From the optical and K’-band images we found evidence that the active component of IRAS F14394$`+`$5332 possesses two nuclei. IRAS F14394$`+`$5332 and 3 other systems identified by Kim et al. (2001) in the 1 Jy sample of ULIRGs by imaging and spectroscopic data can be classified as multipled mergers. Recent survey with the HST images of 113 ULIRGs by Borne et al. (2000) found that the fraction of multiple merger among ULIRGs might be as high as 20%. However, this percentage could be overestimated since Borne et al. (2000) did not identify real interacting members with spectroscopic observations. Studies of the ULIRGs and luminous infrared galaxies (LIG), which are less luminous in the infrared than ULIRGs, by Gao (1996, 1997) also found that ULIRGs and LIGs have strong clustering. They suggested that some LIGs and ULIRGs were formed by the merger of groups of gas-rich galaxies.

Other systems studied in our paper appear to be in the early stage of interaction. Thus, these systems do not fit into the standard interpretation of the ULIRGs because the star-formation activity is expected to happen at the late stage of interaction. Murphy et al. (1996) also found 7 ULIRGs in their sample with separation greater than 10 kpc. They suggest several possible scenarios to explain these particular systems such as the presence of a third nucleus from previous encounter or the ULIRG phenomenon can occur early in the interaction. If the second scenario applies to the ULIRGs considered here, our detection of molecular gas in only one component of the system indicates that the ULIRG phenomenon might be triggered during the early phase of interaction in a system comprised of a gas-rich galaxy and a much less gas-rich galaxy.

The IRAS F23327$`+`$2913 system is unique in the whole 1 Jy sample of the ULIRGs compiled by Kim & Sanders (1998). The presence of absorption features in the optical spectrum of the northern component suggests that the N component did not experience star formation episode for the last $`10^9`$ yrs at least. The molecular gas is detected in the Southern component which has a LINER spectrum. Thus, spectroscopic evidence suggests that the S component is more active than the N companion. More surprisingly, in the R-band image the S component appears as a spiral galaxy with a large bar-like structure at the center. Therefore, the morphology of the S component is very different from other ULIRGs which have either elliptical appearance or strongly disturbed morphology (Sanders et al. 1988, Borne et al. 2000, Kim et al. 2001). Numerical simulations of the merger involving two massive disk galaxies (Mihos & Hernquist 1996, Springel 2000) showed that the two galaxies are strongly disturbed during the interacting phase and the final product of the merger resembles an elliptical galaxy. Therefore, IRAS F23327$`+`$2913 does not fit into the commonly accepted scenario of merger between two disk galaxies for ULIRGs. Other mechanisms should be considered to explain the presence of the large reservoir of molecular gas as well as the huge far-infrared luminosity of the S component. The extended 6 cm radio continuum emission reported by Crawford et al. (1996) indicates that the star formation activity is not concentrated in the nuclear region as in other ULIRGs (Downes & Solomon 1998, Bryant & Scoville 1999) but spread over the galaxy. The existence of the thick bar in a system of two merging disk galaxies can be observed in the early stage of the tidal interaction in the numerical simulations (eg. Noguchi 1988, Barnes 1992, Barnes & Hernquist 1996). If this is the case, then we are witnessing the earliest stage of interaction in ULIRGs. In order to become an ULIRG at this early phase of the interaction, the S component should have a large amount of gas initially or the timescale to collect gas into the nuclear region is very short. In addition, the presence of a bar-like structure suggests that the molecular gas might be driven to the center of the galaxy under the influence of a non-axisymmetric gravitational potential. However, Ho et al. (1997) and Mulchaey & Regan (1997) found that Seyfert and LINER galaxies do not show any significant preference for bars. As a result, it’s still unclear whether a bar is effective enough to drive the large amount of molecular gas currently found in the S component. Higher resolution observations of CO emission lines and radio continuum are needed to better understand the origin of this galaxy.

## 6 Conclusion

We have detected CO(J=1-0) emission line in 5 out of 6 widely separated ULIRGs selected from the complete 1 Jy sample. The molecular gas mass present in these ULIRGs indicates that they are gas-rich. Our high angular resolution observations also reveal the concentration of the molecular gas in only one galaxy of the pair. Optical spectra suggest that galaxy is also the most active in the pair, having either LINER or Seyfert II spectral type. Using optical and K’-band images of the ULIRGs, we found evidence in support of the multiple merger scenario for IRAS F14394+5332. Other systems are likely in the early phase of interaction. In particular, optical image of IRAS F23327$`+`$2913 reveals that a large amount of molecular gas is concentrated in a spiral-like galaxy with no evidence for strong tidal interaction. Our observation results suggest that the conventional view of ULIRG phenomenon as the final phase of the merger of two massive disk galaxies probably does not apply to the widely separated ULIRGs. Other mechanisms are needed to explain their origin. Future sensitive and high resolution observations will help to clarify the nature of these unusual ULIRGs.

We would like to express our thanks to an anonymous referee for critical comments which are very helpful in improving our paper. We also thank Dr. S. Veilleux for providing us the data prior to publication.