text:
# 1 Introduction

## 1 Introduction

One of the outstanding problems of modern cosmology is the understanding of structure formation in the universe . Visible matter at astronomical scales appears organized in a hierarchy of galaxies, cluster and supercluster of galaxies, which tend to be found in filamentary aggregates or two-dimensional sheets that encompass large regions with much lower density of matter and structures. These regions, called voids, occupy a large fraction of the observed universe, and are also organized in a hierarchy . The theoretical investigation of void formation is receiving increasingly attention .

There is general agreement that the galaxy two point correlation function is scale invariant (self-similar), at least at not too large scales, decaying as a power of the distance, $`1/r^\gamma `$. The exponent $`\gamma `$ as well as the scale of homogeneity is still controversial. Pietronero and coworkers found $`\gamma 1`$ and claim that the power law is obeyed up to the deepest distances . Other authors gave a different exponent ($`\gamma 1.8`$) and support the existence of an observed scale of homogeneity . Self-similar behaviour has also been found in the interstellar medium .

Although certainly the dynamics must be very important in order to explain these facts, and many physical effects might play a prominent role, it is well possible that many aspects of the observed structures may be understood in terms of the equilibrium states of self-gravitating matter alone, as claimed by the authors of Refs. . Indeed, it has been tried to apply thermodynamics to astrophysical systems since a long time (see for instance Ref. ). In this paper, we will analyze the general phase structure of a non-relativistic self-gravitating system at thermal equilibrium. The results are general and can be applied to any such system that can be considered to be in these conditions. We will not be concerned here with the very interesting questions of whether thermal equilibrium can be reached in such systems, or the way it is attained. We just assume that the self-gravitating system is at thermal equilibrium.

Let us consider a system of $`N`$ classical particles of mass $`m`$ confined in a region of volume $`V`$ and interacting each other via the Newtonian gravitational potential. Its hamiltonian reads

$$=\frac{1}{2m}\underset{i=1}{\overset{N}{}}𝒑_i^2\frac{Gm^2}{2}\underset{ij}{}\frac{1}{|𝒓_i𝒓_j|}.$$
(1)
The thermodynamics of such a system is ill defined: the entropy does not exist due to the singularity of the potential at short distances . Furthermore, even if one modifies the potential to remove the short distance singularity, the usual thermodynamical limit does not exist, since the thermodynamic potentials are not extensive, due to the long range gravitational force. In such cases, the microcanonical specific heat can become negative and different statistical ensembles are not equivalent (see Refs. ). In particular, the grand canonical ensemble is dominated by completely collapsed configurations, whatever the chemical potential. To use the ordinary thermodynamical tools, the potential must also be modified at long distances111It has been proposed that a kind of thermodynamical limit for systems with potentials decaying as $`1/r`$ could be taken by considering the so called dilute regime , but it can be shown that this statement cannot hold ..

The thermodynamics of self-gravitating system has been studied since a long time by confining a finite system on a finite box and using mean field theory. The approach developed in this paper is different: we will describe the thermodynamical properties of a self-gravitating system as the limiting case of a family of well behaved, short range systems, the interactions of which decay with distance more and more slowly, thus resembling more and more the newtonian $`1/r`$ potential. For these systems the usual thermodynamical limit, which is considered in this paper, does exist. To carry out this investigation, we will take advantage of the fact that the statistics of a self-gravitating ensemble of particles can be related to an euclidean field theory of a single scalar field , in a similar way to the relation between the Coulomb gas and the Sine-Gordon field theory . Hence, the remaining of the paper will rely on the techniques of euclidean field theory.

## 2 Field theoretical description of a self-gravitating system

The statistical mechanics of a self-gravitating system can be studied in the usual way if the gravitational potential is regularized at short and long distances. Let us choose as a regularized gravitational potential, denoted by $`U_\text{R}(r)`$, an attractive Yukawa potential, with range $`1/m_0`$, endowed with a hard core of size $`r_0`$ at short distances: $`U_\text{R}(r)=Gm^2\mathrm{e}^{m_0r}/r`$ for $`r>r_0`$ and $`U_\text{R}(r)=\mathrm{}`$ for $`r<r_0`$. After integrating out the momenta, the grand canonical partition function for chemical potential $`\mu `$ and temperature $`T`$ can be written as

$$𝒵_{\text{GC}}=\underset{N}{}\mathrm{e}^{\mu N}\left(\frac{ma^2}{2\pi \beta \mathrm{}^2}\right)^{3N/2}\overline{𝒵}_N,$$
(2)
where $`\beta =1/k_BT`$, $`a`$ is a constant with units of length, to be specified below, and

$$\overline{𝒵}_N=\frac{a^{3N}}{N!}\underset{i=1}{\overset{N}{}}d^3r_i\mathrm{exp}\{\frac{\beta }{2}\underset{ij}{}U_\text{R}(|𝒓_i𝒓_j|)\}.$$
(3)
The partition function $`\overline{𝒵}_N`$ can be approximated by dividing the space volume in $`\mathrm{\Lambda }=V/a^3`$ cells of size $`a`$ (the unit length introduced above), and replacing the integrals by appropriate Riemann sums, which can be reorganized as a sum over cell occupation numbers $`S_i`$, with the index $`i`$ running from $`1`$ to the number of cells, $`\mathrm{\Lambda }`$. If $`a`$ is of the order of the hard core size, each cell can be either void or occupied by one particle: $`S_i=0,1`$. This is similar to consider self-gravitating fermions . We get

$$\overline{𝒵}_N=\underset{S_1=0}{\overset{1}{}}\mathrm{}\underset{S_\mathrm{\Lambda }=0}{\overset{1}{}}\mathrm{exp}\left\{\frac{b^2}{2}\left[\underset{n,n^{}}{}U_\text{L}(nn^{})S_nS_n^{}U_\text{L}(0)N\right]\right\}\delta (\underset{n}{}S_nN),$$
(4)
where $`b^2=4\pi \beta Gm^2`$, $`U_\text{L}(nn^{})`$ is a proper lattice version of $`U_\text{R}(𝒓𝒓^{})/(4\pi Gm^2)`$, the factor $`1/N!`$ has been canceled due to the distinguishability of classical particles, and the Kronecker delta takes into account that the number of particles is fixed .

Since the Yukawa potential is the Green function of $`^2+m_0^2`$, we can choose its lattice version as the lattice Green function of $`_\text{L}^2+m_0^2`$, where $`_\text{L}^2`$ is a discretization of the laplacian:

$$U_\text{L}(nn^{})=\frac{1}{a^3}\left(\frac{1}{_\text{L}^2+m_0^2}\right)_{n,n^{}}.$$
(5)
Using (5) and the Hubbard-Stratonovich formula, we can write the following identity:

$`\mathrm{exp}\left\{{\displaystyle \frac{b^2}{2}}{\displaystyle \underset{n,n^{}}{}}U_\text{L}(nn^{})S_nS_n^{}\right\}=`$
$`C{\displaystyle [d\psi ]\mathrm{exp}\left\{\frac{1}{2}\underset{n,n^{}}{}a^3\psi _n\left(_\text{L}^2+m_0^2\right)_{nn^{}}\psi _n^{}+b\underset{n}{}\psi _nS_n\right\}},`$ (6)
where $`C`$ is a number independent of $`S_n`$. Plugging Eq. (6) in (4) and the resulting equation for $`\overline{𝒵}_N`$ in the expression for $`𝒵_{\text{GC}}`$, and performing the summation over $`S_i`$, we can write the grand canonical partition function of the (regularized) self-gravitating system in terms of a local euclidean field theory with a single scalar field $`\psi _n`$:

$$𝒵_{\text{GC}}=C[d\psi ]\mathrm{exp}[𝒮],$$
(7)
the action of which is given by

$$𝒮=\frac{1}{2}\underset{nn^{}}{}a^3\psi _n𝒟_{nn^{}}\psi _n^{}\underset{n}{}\mathrm{ln}\left(1+g\mathrm{e}^{b\psi _n}\right),$$
(8)
where $`g=\mathrm{e}^\mu \mathrm{e}^{\frac{b^2}{2}U_\text{L}(0)}(\frac{ma^2}{2\pi \beta \mathrm{}^2})^{3/2}`$ and $`𝒟=_\text{L}^2+m_0^2`$. Notice the unusual form of the action: the interaction term, $`\mathrm{ln}(1+g\mathrm{e}^{b\psi })`$, is unbounded from below. It behaves as $`b\psi `$ for $`\psi \mathrm{}`$. The action, however, is bounded from below if $`m_0^2>0`$. This reflects the fact that the grand canonical ensemble does not exist for non-extensive (long range and/or purely attractive) systems. The infrared and ultraviolet cutoffs, $`m_0`$ and $`a`$ respectively, make the system short ranged and repulsive at short distances (hard core of size of the order of $`a`$), and thermodynamics is well defined. In the naive continuum limit, in which $`a0`$ ($`g0`$) and the relevant field configurations are smooth, we recover the action of Ref. . This continuum action is not bounded from below and hence the corresponding functional integral diverges. Thus, the results of Ref.  are at most formal and the conclusions different from ours. The regulators, $`a`$ and $`m_0`$, play an essential role, although the conclusions, as we will see, are independent of them (provided they are small enough).

The action contains four parameters: the lattice spacing, $`a`$, which obviously has dimension of length; the inverse of the potential range, $`m_0`$, which has dimension of $`(\text{length})^1`$; $`b`$, with dimension of $`(\text{length})^{1/2}`$; and $`g`$, which is dimensionless. The field is canonically normalized, so that its dimension is $`(\text{length})^{1/2}`$. It is convenient to redefine the field such that $`\varphi _n=b\psi _n+\mathrm{ln}g`$, and to work in terms of the dimensionless quantities $`\kappa =a/b^2`$, $`\overline{m}_0^2=a^2m_0^2`$, and $`h=1/2+\kappa \overline{m}_0^2\mathrm{ln}g`$. Ignoring constant terms, the action reads

$$𝒮=\frac{\kappa }{2}\underset{nn^{}}{}\varphi _n\overline{𝒟}_{nn^{}}\varphi _n^{}\underset{n}{}\mathrm{ln}\mathrm{cosh}\frac{\varphi _n}{2}h\underset{n}{}\varphi _n,$$
(9)
where $`\overline{𝒟}=a^2_\text{L}^2+\overline{m}_0^2`$. An action of the same form, with $`\overline{m}_0^2=0`$, has been obtained in Ref.  in a similar way. In that case, however, the parameter $`h`$ has a different meaning. Notice that for $`h=0`$ the action (9) has the symmetry $`\varphi _n\varphi _n`$. Since the term linear in $`\varphi _n`$ breaks this symmetry, we will have $`\varphi >0`$ for $`h>0`$ and $`\varphi <0`$ for $`h<0`$. At $`h=0`$ we will have $`\varphi =0`$ unless the symmetry is spontaneously broken.

The average of particle and energy densities, $`\rho `$ and $`ϵ`$, can be obtained as derivatives of the grand canonical partition function. This allows to express them in terms of the field $`\varphi `$ in the following way:

$`\rho `$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{n}{}}{\displaystyle \frac{\mathrm{e}^{\varphi _n}}{1+\mathrm{e}^{\varphi _n}}},`$ (10)
$`ϵ`$ $`=`$ $`k_BT\left[{\displaystyle \frac{3}{2}}\rho {\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{V}}{\displaystyle \underset{n}{}}[\varphi _n\mathrm{ln}g\kappa ^1\overline{𝒟}_{nn}^1]{\displaystyle \frac{\mathrm{e}^{\varphi _n}}{1+\mathrm{e}^{\varphi _n}}}\right].`$ (11)
One can easily identify the first contribution to the energy in Eq. (11) as the kinetic energy and the second one as the potential energy. The later has the form $`_n(\varphi _n\mathrm{ln}g\kappa ^1\overline{𝒟}_{nn}^1)\rho _n`$, where $`\rho _n=\mathrm{e}^{\varphi _n}/(1+\mathrm{e}^{\varphi _n})`$ is the local particle density, what implies that the field $`\varphi _n`$ represents, up to an additive constant, $`\mathrm{ln}g+\kappa ^1\overline{𝒟}_{nn}^1`$, the local gravitational potential222Note that $`\overline{𝒟}_{nn}^1`$ is independent of $`n`$ due to the translational invariance of the Laplacian..

After some algebraic manipulations, which are outlined in the appendix since they are not completely straightforward, we can write the following exact equation for the average energy, which has the expected form:

$$ϵ=k_BT\left[\frac{3}{2}\rho \frac{1}{2}\frac{1}{V}\underset{n}{}\underset{r0}{}\kappa ^1\rho _n\overline{𝒟}_{n,n+r}^1\rho _{n+r}\right].$$
(12)
In the above equation $`ϵ`$ and $`\rho `$ stand for the average energy and density, respectively. Remember that $`\overline{𝒟}_{n,n+r}^1\mathrm{e}^{\overline{m}_0r}/r`$. Notices that the fact that the term with $`r=0`$ is excluded from the sum over $`r`$ in Eq. (12) implies that the contact term does not contribute to the energy. This, and the bounds $`0<\rho _n<1`$, are manifestations of the particle hard core.

## 3 Phase diagram

Perturbation theory in euclidean field theory starts by identifying the minimum of the classical action and assuming that the relevant field configurations are small fluctuations around this minimum. For the action (9) this is a good approximation if $`\kappa `$ is large. We shall argue below that it is indeed a good approximation whatever $`\kappa `$ if $`\overline{m}_0^2`$ is small.

The minimum of the action obviously corresponds to the constant field that minimizes the classical potential

$$𝒰=\frac{\kappa \overline{m}_0^2}{2}\varphi ^2\mathrm{ln}\mathrm{cosh}\frac{\varphi }{2}h\varphi ,$$
(13)
and, therefore, satisfies the equation

$$\kappa \overline{m}_0^2\varphi =\frac{1}{2}\mathrm{tanh}\frac{\varphi }{2}+h.$$
(14)
The above equation has one, two or three solutions, depending on the values of $`\kappa \overline{m}_0^2`$ and $`h`$. When it has one solution, it correspond to the global minimum; if there are two solutions, one is the global minimum and the other one a point of inflection; three solutions correspond to a local maximum, a local minimum, and the global minimum. In the last case, the minima can only be degenerate if $`h=0`$, and then symmetry implies that Eq. (14) has either one or three solutions: for $`\kappa \overline{m}_0^2>1/4`$ the only solution is $`\varphi =0`$, while for $`\kappa \overline{m}_0^2<1/4`$ we have two solutions, $`\varphi =\pm \varphi _00`$ (with $`\varphi _0`$ positive), which are the two degenerate global minima, besides the symmetric solution $`\varphi =0`$, which is a local maximum. We will denote the two global minima by $`\varphi _h=+\varphi _0`$ and $`\varphi _l=\varphi _0`$. The densities of each phase are respectively

$$\rho _l=\frac{1}{a^3}\frac{\mathrm{e}^{\varphi _l}}{1+\mathrm{e}^{\varphi _l}},\rho _h=\frac{1}{a^3}\frac{\mathrm{e}^{\varphi _h}}{1+\mathrm{e}^{\varphi _h}}.$$
(15)
Hence, the solutions $`\varphi _l`$ and $`\varphi _h`$ describe phases of low and high density respectively.

The phase diagram in the plane $`(\kappa \overline{m}_0^2,g)`$, at the classical level, is displayed in figure 1. For $`\kappa \overline{m}_0^2>1/4`$, the transition between the low (small $`g`$) and high (large $`g`$) density regimes is smooth. For $`\kappa \overline{m}_0^2<1/4`$, the low density phase is separated from the high density phase by a first order transition which takes place at $`g_{\text{FO}}=\mathrm{exp}[1/(2\kappa \overline{m}_0^2)]`$ (i.e., $`h=0`$). This first order line ends at the critical point $`\kappa _c\overline{m}_0^2=1/4`$, $`g_c=1/\mathrm{e}^2`$ ($`h=0`$). The order parameters, $`\delta \varphi =\varphi _h\varphi _l`$ and $`\delta \rho =\rho _h\rho _l`$, vanish at the critical point with the classical (mean field) exponent $`1/2`$. This mean field critical behavior is a consequence of the classical approximation and could be modified by the neglected fluctuations.

For small $`\kappa \overline{m}_0^2`$ and $`h=0`$, Eq. (14) gives $`\varphi _h1/(2\kappa \overline{m}_0^2)`$, $`a^3\rho _h1\mathrm{e}^{1/(2\kappa \overline{m}_0^2)}`$, $`\varphi _l1/(2\kappa \overline{m}_0^2)`$, and $`a^3\rho _l\mathrm{e}^{1/(2\kappa \overline{m}_0^2)}`$. Thus, the low density phase is very dilute, while the high density phase is extremely dense. Most of the densities (including presumably those of physical systems) are between $`\rho _l`$ and $`\rho _h`$, and they correspond, therefore, to the phase coexistence region.

## 4 Critical behavior and correlations

The corrections to the classical approximation can be obtained perturbatively, by expanding the action around the corresponding minimum. For each phase we write $`\varphi _n=\varphi _h+\phi _n`$ and $`\varphi _n=\varphi _l+\phi _n`$ respectively, and, ignoring again constant terms, the corresponding actions can be written as

$$𝒮=\frac{\kappa }{2}\underset{nn^{}}{}\phi _n\overline{𝒟}_{nn^{}}\phi _n^{}\pm \lambda \underset{n}{}\phi _n\underset{n}{}\mathrm{ln}\left[1+\lambda (\mathrm{e}^{\pm \phi _n}1)\right],$$
(16)
where the signs $`+`$ and $``$ correspond to the low and high density phases respectively, and $`\lambda =1a^3\rho _h=a^3\rho _l`$.

When $`\kappa \overline{m}_0^2`$ is very small, $`\lambda `$ is much smaller and the actions for either phase (16) are very close to the gaussian critical point, $`\lambda =0`$. Therefore, the classical approximation we are using is very good. In addition, we have a very interesting situation: a very sharp first order transition that separates two phases which, in turn, are close to criticality. Hence, the actions (16) produce a very large correlation length and critical (self-similar) behavior over a vast range of scales. The question is to which universality class corresponds such critical behavior. In three dimensions, the gaussian fixed point is infrared unstable under perturbations of relevant operators. This means that the critical behavior of actions defined in the neighborhood of the gaussian fixed point can be governed by a different (non-gaussian) fixed point. For this to happen, the couplings of the relevant operators of dimension larger than one must be of the order of $`\kappa \overline{m}_0^2`$. In our case $`\lambda \kappa \overline{m}_0^2`$ and we are in the opposite case: the actions (16) lie very close to the renormalized trajectory of the gaussian fixed point, given by $`\lambda =0`$, $`\kappa \overline{m}_0^2>0`$. Hence, the critical behavior is governed by the gaussian fixed point, and thus belongs to the mean field universality class.

The correlations can be computed to leading order in $`\lambda `$ with good accuracy. For the field we get $`\varphi _𝒓\varphi _𝒓^{}\varphi _𝒓\varphi _𝒓^{}=G(𝒓𝒓^{})`$, with

$$G(𝒓𝒓^{})=\frac{d^3p}{(2\pi )^3}\frac{\mathrm{e}^{\mathrm{i}𝒑(𝒓𝒓^{})}}{p^2+\kappa \overline{m}_0^2},$$
(17)
where we have used the continuum expression since the lattice spacing is much smaller than the correlation length, $`\xi =a/(\kappa \overline{m}_0^2)^{1/2}`$. At large distances, $`G(𝒓)`$ decays exponentially as $`\mathrm{exp}(|𝒓|/\xi )`$, but in a wide range of distances, $`a|𝒓|\xi `$, it is approximately self-similar, $`G(𝒓)1/|𝒓|`$.

The correlations of the density (10), $`\mathrm{\Gamma }(𝒓𝒓^{})`$, have the same long distance behavior as the correlations of the field. To leading order in $`\lambda `$ we have:

$$\mathrm{\Gamma }(𝒓𝒓^{})=\lambda ^2\mathrm{e}^{G(0)}\left[\mathrm{e}^{G(𝒓𝒓^{})}1\right].$$
(18)
For large distances $`\mathrm{\Gamma }(𝒓)G(𝒓)`$, and therefore the density correlations behave as $`1/|𝒓|`$ for $`a|𝒓|\xi `$. Hence, the present approach to self-gravitating systems predicts that density correlations decay as a power law with exponent $`\gamma =1`$ over a vast range of scales. Since $`\xi `$ is proportional to the assumed range of the gravitational interaction, it is well possible that correlations be self-similar at any observable scale.

It is worthwhile stressing that the above analysis refers to the quasi-critical behaviour of each of the two phases as $`\kappa \overline{m}_0^20`$. There is a true critical point at $`\kappa _c\overline{m}_0^2=1/4`$ and $`g=1/\mathrm{e}^2`$ (or $`\rho _c=1/2`$). In this case the parameters are not small and a non-perturbative analysis is required in order to investigate the nature of its universality class. This critical point might be relevant at very high temperatures, when $`\kappa \overline{m}_0^2`$ is of order one.

A similar analysis can be made for small $`\kappa `$ (strong coupling regime). The symmetry is broken at the classical level if $`\kappa \overline{m}_0^2<1/4`$. For small $`\kappa `$ and small $`\overline{m}_0^2`$ short wavelength fluctuations cannot induce tunneling between the two minima since it would imply huge fluctuations of the action, of the order $`1/\kappa \overline{m}_0^4`$. Likewise, long wavelength fluctuations cannot induce tunneling since the barrier is too high, of the order of $`1/(8\kappa \overline{m}_0^2)`$. The two minima describe phases of low and high density. The correlation length on each phase is of the order of $`1/\overline{m}_0^2`$. Hence, the conclusions are the same as in the weak coupling regime: an abrupt first order transition separates the low and high energy phases, and each phase present critical behaviour of mean field type.

## 5 Finite volume effects: the Lane-Emden equation

It might seem paradoxical that the mean field solution is homogeneous, since it is well known that the usual mean field solutions of self-gravitating systems present density profiles that depend on the distance to the center of the system. There is a simple explanation: since we are looking for the mean field solution of a short ranged translational invariant system in the infinite volume limit it is natural that it be homogeneous. Had we looked for the minimum of the action (8) on a finite box of size $`L`$ we would have obtained the equation

$$\underset{n^{}}{}(_L^2)_{nn^{}}\psi _n^{}+m_0^2\psi _n\frac{bge^{b\psi _n}}{1+ge^{b\psi _n}}=\mathrm{\hspace{0.33em}0}.$$
(19)
On a finite volume, the walls break translational invariance and the solution of Eq. (19) will not be homogeneous, due to the boundary conditions imposed to the operator $`(_L^2)_{nn^{}}`$. Since this operator is a discretization of the Laplacian, if $`L1/m_0`$, Eq. (19) is similar to the isothermal Lane-Emden equation  (if $`m_0=0`$, to leading order in $`g`$ we will have exactly the discretized Lane-Emden equation, and the corrections in $`g`$ are due to the cut-off $`a`$), and we will get a solution similar to the known profiles of self-gravitating systems in the mean field approximation. On the other hand, if $`L1/m_0`$, the spatial dependence of the solution will be washed out and we will have the reported constant solutions. However, the two solutions correspond to extremely high and extremely low densities, respectively. The system will only be homogeneous at such extreme densities. For intermediate densities, inhomogeneities will necessarily develop: the system will be in a mixture of high and low density domains.

This instability as the volume increases had been notice long ago by Antonov, who found that the solutions of the Lane-Emden equation ceases to be (local) maxima of the entropy if the size of the box is larger than $`0.335GM^2/(E)`$, where $`M`$ is the total mass and $`E<0`$ the energy, and the system collapses . This was called the gravothermal catastrophe in Ref. .

## 6 Final remarks

In thermal equilibrium, a self-gravitating system is made up of domains of low and high density (voids and clusters). The distribution of domains is a dynamical question that has to do with the way equilibrium is reached. The correlations within each domain are self-similar on a vast range of scales, decaying as $`1/r`$. The transition to homogeneity333The scale at which the density fluctuations cease to be self-similar.  would take place on scales comparable to the range of the gravitational interaction, which may be larger than the deepest observed distances. A similar behavior (first order phase transition and self-similar correlations) was found for the two-dimensional self-gravitating system by using techniques of conformal field theory .

It is straightforward to take into account the cosmic expansion in a simplified way, as in Refs. , by introducing comoving coordinates $`𝒙`$ in (1), such that the physical coordinates are $`𝒓=R(t)𝒙`$, where $`R(t)`$ is the scale factor. As intuitively expected, this is equivalent to rescale the lattice spacing: $`a(t)=R(t)a`$ and, consequently, we have the following rescaling of parameters: $`\kappa R(t)\kappa `$ and $`\overline{m}_0^2R^2(t)\overline{m}_0^2`$. This variation of the parameters implies that the difference between the densities of the two phases decreases. This is not surprising: the expansion of the universe acts as a pressure that competes with the tendency of gravity towards collapse.

It is remarkable that the picture of the self-gravitating system devised in this work is strikingly similar to the observed universe, despite it is not at thermal equilibrium. The scaling exponent of the the correlation function of the density, $`\gamma =1`$, however, is far from the widely accepted exponent of the galaxy correlation function, $`\gamma 1.8`$ for distances between 0.2 and 20 Mpc . As mentioned in the introduction, in the last years there has been a strong debate about the validity of this result. Pietronero and his group analyzed the data from a different perspective and claimed that the exponent of the galaxy correlation function is $`\gamma 1`$, and that the self-similar behavior extends up to the deepest explored distances. The controversy seems not resolved, although most astrophysicist believe the earlier result, $`\gamma 1.8`$. This is theoretically supported by the thermodynamic arguments given in . These arguments, however, are somehow heuristic and their validity may be questioned. Indeed, using Eq. (12) as starting point, the same arguments can be applied step by step to the formulation of the self-gravitating system given in this work, leading to the same predictions of Saslaw’s book. However, the results of this work, based on a rather rigorous treatment of the self-gravitating system, are completely different.

Hence, the conclusion of this paper is clear: either Pietronero and coworkers are right, and $`\gamma 1`$, or the thermodynamical approach is not valid to describe the universe at such scales. The later possibility cannot be discarded, since it is difficult to argue that the universe is at thermal equilibrium, and dynamical effects may play a prominent role in the behavior of the galaxy correlation function. In any case, the results obtained in this paper need not be applied to the universe as a whole, or to the large scale structure of it, but to any self-gravitating system that are close to thermal equilibrium. The interstellar gas, where self-similar behavior has also been observed, might be an instance. In this case the scaling exponent of the density correlations is compatible with $`\gamma =1`$, but with large uncertainties .

Acknowledgments

I thank J.L. Alonso, J.L. Cortés, and A. Galante for useful discussions. This work received financial support from Ministerio de Ciencia y Tecnología (Spain) under the Ramón y Cajal program, and from CICyT (Spain) under the project FPA 2000–1252.

## Appendix

Let us outline in this appendix the derivation of equation (12), which is not completely straightforward. The field $`\varphi _n`$ can be obtained from the derivative of the action as

$$\varphi _n=\kappa ^1\underset{n^{}}{}\overline{𝒟}_{nn^{}}^1\rho _n^{}+\kappa ^1\underset{n^{}}{}\overline{𝒟}_{nn^{}}^1\frac{𝒮}{\varphi _n^{}}\frac{1/2h}{\kappa \overline{m}_0^2},$$
(20)
where we used $`_n^{}\overline{𝒟}_{nn^{}}^1=1/\overline{m}_0^2`$. Now insert the above expression for $`\varphi _n`$ in equation (11)444Not in the term $`\rho _n=\mathrm{e}^{\varphi _n}/(1+\mathrm{e}^{\varphi _n})`$.. The term proportional to $`(1/2h)/\kappa \overline{m}_0^2`$ cancels the $`\mathrm{ln}g`$ term of Eq. (11). Averaging over the thermal fluctuations the resulting expression for $`ϵ`$, we are led to compute $`\rho _n𝒮/\varphi _n^{}`$, which can be written as

$$\rho _n\frac{𝒮}{\varphi _n^{}}=\frac{1}{𝒵}[d\varphi ]\rho _n\frac{}{\varphi _n^{}}\mathrm{e}^𝒮.$$
(21)
Integrating by parts and taking into account that $`\rho _n/\varphi _n^{}=\delta _{nn^{}}\rho _n(1\rho _n)`$ , we get

$$\rho _n\frac{𝒮}{\varphi _n^{}}=\delta _{nn^{}}\rho _n(1\rho _n).$$
(22)
The term linear in $`\rho _n`$ cancels the term proportional to $`\kappa ^1\overline{𝒟}_{nn}^1`$ of the averaged Eq. (11). Hence, we obtain for the averaged energy

$$ϵ=k_BT\left[\frac{3}{2}\rho \frac{1}{2}\frac{1}{V}\underset{nn^{}}{}\kappa ^1\rho _n\overline{𝒟}_{nn^{}}^1\rho _n^{}+\frac{1}{2}\frac{1}{V}\underset{n}{}\kappa ^1\rho _n\overline{𝒟}_{nn}^1\rho _n\right].$$
(23)
It suffices to write $`n^{}=n+r`$ to realize that the above equation is Eq. (12).


format: text

text:
# New Problems for a New Century

## 1 Introduction

This is an optimistic time in X-ray astronomy. Chandra and XMM-Newton are fully operational and producing results at an undigestible rate; Astro-E has been rescheduled for a 2005 launch and will then assume its place as the third essential component of contemporary missions. HETE-II is reportedly performing flawlessly and Swift is on track for a 2003 launch. Two major international missions, Constellation-X and XEUS, as well as a host of smaller projects are at various stages of development.

My task is to provide some overview of what we have learned here and, especially, to suggest some questions that should be addressed in planning new missions. I will have to be somewhat selective and idiosyncratic in what I discuss and will refer to the accompanying articles, using brackets, for more extensive discussions. As the meeting has been organised along astronomical lines, I will attempt an orthogonal organization and try to bring out some connections by concentrating on the underlying physics.

However, before doing this, I would like to recall the scientific context in which Chandra, XMM-Newton and Astro-E were first proposed. The big issues in 1977, when Chandra was proposed were the X-ray background and quasars, the two being thought to be closely related. In addition, there was an understandable envy of optical astronomers and a strong belief that arcsecond imaging would transform our view of the high energy universe. For XMM-Newton, in 1982, the prospect of high dispersion spectroscopy and the ability to perform abundance analyses and velocity measurements were strong motivators. By the time that Astro-E was conceived, around 1991, the lure of the hard X-ray band and, in particular, iron lines was irresistible and it was appreciated that, despite their lower luminosity, the stronger fluxes from Seyfert galaxies and their more rapid variability recommended them as suitable sources with which to study the general properties of AGN. In fact, as we have seen here, these have turned out to be pretty good scientific justifications for the two missions that have been launched and the one whose successor we eagerly await.

X-ray astronomy has become a far richer field than anticipated by these proposals. While some problems, like the X-ray background, supernova remnants, clusters and accretion disks around cokpact objects appear to have been broken, albeit leaving a lot of crucial details to be filled in, others, like GRBs have become extremely interesting and have developed in quite unexpected directions. Completely new phenomena like the Galactic Ridge X-ray Emission and X-ray protostars of all types, which were widely discussed here, have been uncovered. On the technological front there have been tremendous advances culminating in the microcalorimeter arrays which are anticipated to lie at the heart of future advances in X-ray spectroscopy.

If I turn now to the future and look at the Constellation-X \[White\] brochure I see black hole astrophysics, “astro-ecology”, and mapping the distribution of dark matter as the principle scientific challenges. For XEUS \[Jansen\], the emphasis is a little more cosmological and we must add the first groups and quasars, the history of metal production and the astrophysics of the hot intergalactic medium. It will be fascinating to see how these drivers evolve over the coming decade.

## 2 Ten Questions for a New Decade

In order to enter into the spirit of this meeting, I would like to present ten fundamental astrophysical or even physical questions that I believe can be addressed by X-ray astronomy.

### 2.1 How do Cosmic Magnetic Fields really Behave?

Cosmic magnetic field, on all scales is generally described under the MHD approximation. Physically, MHD rests on two complementary principles, that magnetic field is frozen into moving fluid – “go with the flow” – and that the anisotropic, magnetic stress tensor is be responsible for directing the flow – “push-pull”. However, even if these principles suffice, almost always and everywere, they do not tell the whole story. Shocks, reconnection sites and finite Larmor radius effects can all lead to behavior that is not consistent with the precepts with the precepts of MHD. It is only though observations, interpreted using numerical simulation, that we are going to understand how magnetized fluid really behaves under cosmic conditions.

X-ray observations reported here are providing vital clues in a variety of environments. All four types of protostar are found to be X-ray sources along with brown dwarfs \[Kuboi, Linskey\], suggesting that coronae are commonly energised, presumably magnetically, to form high temperature gas. Disk accretion, likewise, shows signs of magnetic activity \[Mukai, Kubota, Mineshiga\]. At the other end of the scale, many putative black holes in galactic nuclei are found to be spectacularly underluminous, with powers as low as $`10^8L_{\mathrm{Edd}}`$ or $`10^5L_{\mathrm{Bondi}}`$ \[Mushotzky\]. This is a very powerful clue as to how magnetic viscosity develops in accreting systems. A quite different consideration is important in rich clusters of galaxies where the Coulomb scattering mean free path is so large that magnetic field must be invoked to limit the heat transport.

### 2.2 How do Cosmic Plasmas really Behave?

Magnetohydrodynamics is only an approximation to the full, kinetic behavior of plasma, dealing, as it does, with the large scale conservation of mass, momentum and energy. Plasma effects may be vital to understanding the Galactic Ridge emission \[Tanaka\] observed from the Galactic bulge. Collective interactions, like transit time damping, rather than pure two body Coulomb interactions may play a role in electron-ion equipartition. These considerations are of vital importance to understanding stellar coronae \[Linsky\] as well as accretion disk coronae and adiabatic accretion in the underluminous sources found in galactic nuclei. They are of no less importance to understanding the thermal balance of hot gas in clusters and the general intergalactic medium and may provide the explanation for why lines from intermediate states of ioniziation, such as would be expected from cooling gas in ionization equilibrium, are generally not seen in cluster spectra.

In addition, I suspect that plasma physics will ultimately be needed to resolve the controversy concerning the X-ray spectra of Seyfert and LINER accretion disks \[Kahn, Iwasawa, Brandt, Boller, Nandra, Zdziarski, Reeves, Yaqoob\]. However, before this happens some purely observational differences must be resolved. How wide and variable are the Fe lines? How much of the observed spectrum is imprinted by the warm absorber? Are broad lines of lower Z elements really present in the spectrum?

### 2.3 What is the Structure of Collisionless Shocks?

Plasma physics is also crucial to the structure of collisionless shocks which are essential to X-ray astronomy as this is almost the only way to heat gas to high temperature in the interplanetary, interstellar and intergalactic media. (Bulk, dissipative heating may be occuring in accretion disks, but even here, this poses problems.) Shocks also accelerate relativistic electrons. They inevitably arise when bulk speeds become supersonic either in expanding flows like stellar winds (Linsky) or as a consequence of gravitational acceleration during structure formation in the expanding universe.

Observationally, we are having a hard time locating shocks. There are very few supernova remnants where we can be sure that we have identified the outer blast wave \[Peter, Aschenbach, Kamae\]. Similarly, with the accretion shocks around clusters \[Bautz, Arnaud\] and even the Galactic center \[Koyama\]. These shocks are responsible for heating gas, accelerating cosmic rays and, possibly, stretching magnetic field lines by large factors. However, we do not understand electron-ion equipartition, whether cosmic rays actually mediate the shock compression and how and whether fields are amplified. This last point is particularly relevant to the mildly relativistic shocks observed in X-ray afterglows from $`\gamma `$-ray bursts, where almost all contemporary models invoke an immediate and magical post shock field enhancement to a fraction of the equipartion value that is just one of many free parameters in the fitted spectra \[Ricker, Nousek, Fiore\]. We really need a more basic understanding of what is happening and the best approach involves careful analysis and simulation of the glorious observations of supernova remnants shown here.

### 2.4 When and how were the Elements Made?

These same observations provide prima facie evidence that heavy elements are synthesized during supernova explosions \[Iyudin, Petre, Aschenbach\]. However, the details remain elusive. We really have no working model of type II supernovae, cannot tell, for example, whether a neutron star or a black hole has been left behind in Cas A and we are in the embarassing position that we cannot even agree upon the progenitors of type Ia supernovae, that are so vital to contemporary cosmography. X-ray observations really provide an excellent handle on the models, probing as they do the mass cut in the explosion of heavy stars, the extent of the inhomogeneity induced by the explosion, the strength and angular distribution of the wind and the centering which also relates to the recoil speed of the remnant star.

There are other, conjectured sites of nucleosynthesis, especially of r-process elements, like hypernovae and neutron star coalescence \[Nomoto\] and where X-ray observations can be vital, such as the Fe lines reported from GRBs \[Fiore\]. On the large scale, opnservations of distant galaxies can explore abundance gradients and provide a quantitative measure of the history of star formation \[Fabbiano, Iwasawa\]. Indeed, claims of abundance enhancements as large as \[Fe/H\] $`5`$ have been made for the centers of galaxies and in rich clusters \[Ohashi, Bohringer\].

### 2.5 How does Matter Behave at High Density?

One of the more important ways through which X-ray astronomy can repay its debt to fundamental physics is to determine the equation of state of cold nuclear matter at supranuclear density \[Tsuruta, Weisskopf\]. We really have little idea from basic physics as to whether the interiors of neutron stars contain free quarks, pion or kaon condensates or just a proton superconductor plus neutron superfluid as conservatively! assumed. X-ray observations are the way to find the answers. This investigation complements the studies of hot nuclear matter being performed at heavy ion colliders and is quite possibly relevant to the early universe.

The most direct method is to measure the mass-radius relation of neutron stars using estimates of the gravitational redshift and the surface gravity. In practice, this is now appearing to be a more difficult task than originally envisaged due to the absence of spectral lines. As feared, the surface composition is largely H and He and there are no strong lines.

A second approach, is less direct but currently more productive. This involves measuring the cooling rate by comparing the measured surface temperature with the ages of the surrounding remants. As neutron star atmospheres and our understanding of remnant dynamics improve, this should become a far more prescriptive as the cooling rates are strongly sensitive to the interior composition.

A third approach is to use the frequencies of neutron star QPOs, though this will require a much more sophisticated understanding of their dynamics than we have at the moment.

### 2.6 Is General Relativity Correct?

In addition to probing nuclear physics, X-ray astronomers have the opportunity to provide quantitative probes of strong gravity, specifically by detecting physical effect that are predicted by the relativistic description of the spacetime around a spinning black hole \[Yaqoob\]. The Kerr metric is the default description and it will take hard evidence to persuade most relativists that it is wrong or incomplete. However, it logically possible that black holes are more hirsute than conventionally thought and that particle (including photon) motion is affected.

To date most attention has centered around attempts to measure the second parameter, the spin or, equivalently, the specific angular momentum. The width of the Fe line has been taken as demonstrating that the orbiting gas gets close to the horizon and that, consequently, holes spin fast. There are some ongoing concerns about the strength and widths of these lines. In addition, some “Devil’s Advocates” have claimed that they can be formed by gas plunging into the horizon. However, it will be hard to make it convincingly quantitative until we understand the radiative transfer. I suspect that there is more promise in the “diskoseismology”, so comprehensively studied most recently using RXTE \[Swank\]. Here, there have also been claims for the detection of spin on the basis of some quite simple and to me, not yet compelling descriptions of what is going on. However, remarkable patterns are starting to emerge among the frequencies and spin must be a major factor in setting the clock. (The observed QPOs are far too hard for the emission to arise from the disk.) It would be wonderful if we could recapitulate the early history of atomic spectroscopy and create a redundant, theoretical description of these modes that verified the detailed form of the Kerr metric.

### 2.7 What is the Nature of Dark Matter?

An even larger physics challenge is to identify the dark matter and chronicle its role in the development of large scale structure. Despite the, perfectly reasonable, focus on neutral supersymmetric particles, and to a lesser extent, axions, we really don’t understand what it is. X-ray astronomy can contribute in a major way by helping trace the shape of the potential well around ellipticals, groups and clusters \[Bautz, Boehringer, Ohahshi, Arnaud\]. In the last case, velocity dispersions, Sunyaev-Zel’dovich dips and weak lensing can all be combined to produce a comprehensive picture. So far, the most important impression left by the X-ray observations is that many, though not all, clusters are surprisingly inhomogenous despite the apparent regularity embodied in the $`LT`$ relation.

### 2.8 How did the Universe Expand?

On the larger scale, CMB and SNIa observations have produced an impressive body of evidence that our 14 Gyr old universe is flat, mostly dark and accelerating. Presuming that this result holds up, it is partly quite unexpected and even more perplexing. (It should not be forgotten that, for a long while, the most compelling evidence that the Einstein-De Sitter universe was not even a close approximation to the truth, came from X-ray observations of clusters \[Henry, Mitsuda\].)

The long term kinematic goal must surely be to measure the variation of the scale factor $`a(t)`$ with cosmic time in a manner that is independent of dynamical assumption. (We do not know if the so-called dark energy is a manifestation of a shortcoming in the left hand side or the right hand side of Einstein’s Field Equations, or both.) It is the eternal hope of observational cosmology that standard rods and candles will be found extending out to large redshifts and X-ray sources, particularly clusters, may yet turn out to be examples \[Hasinger\].

### 2.9 What were the First Structures in the Universe?

Although, it is commonly believed that normal galaxies assembled from the merger of smaller subunits and that they, in turn congregated into clusters, we still do not what were the first self-gravitating structures to form in the expanding universe. What we do know is that powerful quasars were already functioning when the universe was less than $`800`$ million years old. This implies that black holes assembled very early in the life of galaxies and probably played an important role in galaxy formation. Although many of these black holes are likely to be obscured, their hard X-ray emission can be seen redshifted into the intermediate X-ray energy band.

Alternatively, perhaps the first self-gravitating objects were “Population III” stars, that formed inside small dark matter halos, and which subsequently turned into some of the “Intermediate Black Holes” the powerful X-ray sources with luminosities in excess of the Eddington limit for $`100`$ M stars, seen in nearby galaxies or the X-ray binaries found in globular clusters \[Mushotzky\].

### 2.10 How do you do X-ray Astronomy from the Ground?

Much of this meeting has been concerned with recent ideas and progress in mirror and detector development throughout the entire X-ray band \[Elvis\]. The prospect of pushing focusing optics beyond 100 keV \[Strüder, Takahashi\] and of deploying large microcalorimeter arrays is indeed an exciting one \[Stahle\]. Looking further ahead, the time is surely ripe to start doing polarimetry again and work towards making far more sensitive X-ray polarimeters. (I note that it has recent proved possible to measure X-ray circular polarization from laboratory synchrotrons.)

There is even considerable enthusiasm for attempting X-ray interferometry from space \[Cash\]. In its most elemental form this is effectively a Young’s slit experiment which should appeal, at least in principle, to any physicist. How long it takes to develop a system capable of accomplishing the stated goal of imaging a black hole remains to be seen, but there are several important milestones along the way which could lead to more broadly applicable technology development.

Among several other approaches to understanding the extreme physical conditions within the sources are high energy density experiments at laser facilities, numerical simulation, especially of three (and four) dimensional magnetohydrodynamics and laboratory astrophysical studies to augment our library of wavelengths, oscillator strengths, collision integrals and so on that will be necessary for us to interpret X-ray observations of hot plasma\[Porter\].

It will be interesting to see if the questions that have motivated proposals for the next generation of X-ray observatories fare as well as those that motivated Chandra, XMM-Newton and Astro-E.

## Acknowledgements

I thank Professor Kuneida for the invitation to attend this meeting and his hospitality and several other participants for helpful advice.


text:
# A Common Origin of Neutralino Stars and Supermassive Black Holes

## 1 INTRODUCTION

Dark objects with masses of $`(0.01÷1)M_{}`$ were detected in the Galactic halo when the microlensing of stars from the Large Magellanic Cloud was observed. Brown and white dwarfs and Jupiter-like planets were proposed as possible baryonic candidates. According to the latest data , dark objects must account for about 20% of the hidden mass in the Galactic halo. However, the theory of primordial nucleosynthesis imposes much more stringent constraints on the number of such baryonic objects. Gravitationally bound, noncompact objects, which are also called neutralino stars, were considered among nonbaryonic candidates; these objects can explain some microlensing events with evidence of gravitational lenses being noncompact . The hypothetical noncompact objects are the lightest objects in the hierarchical large-scale structure of the Universe. In the model of Gurevich, Zybin, and Sirota , these are formed immediately after the passage of the Universe to the dust stage.

To reconcile the parameters of noncompact objects with data for the observed microlensing events, Gurevich et al.  assumed the cosmological density perturbation spectrum on some small scale to increase sharply with an rms fluctuation of the order of 1 at the peak. For adiabatic perturbations, rms fluctuations of the order of 0.05 correspond to these fluctuations at the radiation-dominated stage on the cosmological horizon scale. As was shown in , fluctuations with such an rms value give rise to primordial black holes (PBHs) at the radiation-dominated evolutionary stage of the Universe. As we will see below, the mass of the forming PBHs exceeds the mass of noncompact objects by several orders of magnitude. The large difference between the masses of noncompact objects and PBHs stems from the fact that the mass of the radiation contained in a fixed comoving volume at the radiation-dominated stage is much larger than the mass of cold dark matter (CDM) in the same volume at the dust stage.

If the power spectrum of primordial cosmological perturbations is a power law with an index $`n>1`$, then PBHs are formed in a wide range of masses. If, however, the spectrum has a peak on some scale, then PBHs are formed mostly in a narrow range of masses, near the mass that corresponds to this peak. A sharp peak emerges in the fluctuation spectrum if the inflationary potential $`V(\varphi )`$ has a flat segment . Indeed, if the derivative at some value of the scalar field is $`V^{}=dV(\varphi )/d\varphi 0`$, then a peak emerges in the perturbation spectrum on the corresponding scale, because the mean density perturbation on the horizon scale is $`\delta _\mathrm{H}M_{\mathrm{Pl}}^3V^{3/2}/V^{}`$, where $`M_{Pl}`$ is the Planck mass. A similar effect is achieved in inflationary models with several scalar fields . In this case, the spectrum outside the peak can have an ordinary shape, for example, it can be a Harrison–Zel’dovich spectrum, and can give rise to galaxies, their clusters and superclusters according to standard scenarios.

Certain evidence for a deviation of the initial perturbation spectrum from a simple power-law shape, more specifically, for a spectral break near large scales, $`k(0.060.6)h`$ Mpc-1, was obtained in the counting of APM galaxies and in the Boomerang and Maxima cosmic microwave background (CMB) anisotropy measurements . Therefore, there is reason to suggest that the spectrum may also exhibit features on small scales.

Here, we show that if the perturbation spectrum has a peak, then there is a clear relationship between the masses of noncompact objects and PBHs; a PBH mass of the order of $`10^5M_{}`$ corresponds to a noncompact object mass of the order of $`0.1M_{}`$. For noncompact objects to be able to serve as gravitational lenses, they must originate from sufficiently large dark matter density fluctuations. At the radiation-dominated stage, these fluctuations logarithmically increase with time and become nonlinear even at this stage. To study the evolution of the fluctuations at the radiation-dominated stage, we use the nonlinear model proposed by Kolb and Tkachev  and specify the initial data for this model according to the linear theory.

Based on standard astrophysical constraints on the PBH number density in the Universe, we obtained constraints on the rms fluctuations at the peak. In turn, the constraints on the fluctuations give constraints on the radius of a noncompact object. As was shown in , a noncompact object can serve as a gravitational lens and can produce observable microlensing events only if its radius exceeds the Einstein radius for this object by no more than a factor of 10. Stringent constraints on the coefficient of nonlinear contraction for noncompact objects follow from this condition.

The inferred relationship between noncompact objects and PBHs holds only for adiabatic cosmological density perturbations. If the density perturbations are entropic, then even if there are large fluctuations in the dust component, the radiation density on the horizon scale is almost uniform and no PBHs are formed.

It should be emphasized that when talking about a common origin of noncompact objects and PBHs, we have in mind not the relationship between individual noncompact objects and PBHs but the fact that the fluctuations from which the entire collection of noncompact objects and the entire collection of PBHs originate have a common spectrum. If a PBH emerged at some point in space, then a noncompact object can no longer emerge at this point. Conversely, if there is a noncompact object, then no PBH was previously (at the radiation-dominated stage) formed at this point.

## 2 THE FORMATION OF PBHs

The PBH formation criterion was analytically derived by Carr  and confirmed by numerical calculations . The region of space with a density $`\rho >\rho _\mathrm{c}=3H^2/8\pi G`$ can be roughly considered to be part of the closed Universe. Gravitational collapse of this region and the formation of a PBH take place if the relative radiation density fluctuation $`\delta _\mathrm{H}=(\rho \rho _\mathrm{c})/\rho _\mathrm{c}`$ at the time it goes under the horizon satisfies the conditions

$$\delta _\mathrm{c}\delta _\mathrm{H}1,$$
(1)
where $`\delta _\mathrm{c}=1/3`$. The left-hand inequality implies that the radius of the perturbed region at the time $`t`$ its expansion stops exceeds the Jeans radius $`ct/\sqrt{3}`$, while the right-hand inequality corresponds to the formation of a PBH rather than an isolated universe. The mass of the forming PBH in this model is

$$M_{\mathrm{BH}}=\frac{M_\mathrm{H}}{3^{3/2}},$$
(2)
where $`M_\mathrm{H}`$ is the mass within the horizon.

In recent years, numerical experiments have revealed the so-called critical gravitational collapse, during which the mass of the forming PBH is

$$M_{\mathrm{BH}}=AM_\mathrm{H}(\delta _\mathrm{H}\delta _\mathrm{c})^\gamma ,$$
(3)
where $`A3`$, $`\gamma 0.36`$, and $`\delta _\mathrm{c}(0.65÷0.7)`$. The mass (3) can be much smaller than $`M_\mathrm{H}`$. However, as shown in , the PBH mass distribution for critical gravitational collapse is concentrated near $`M_{\mathrm{BH}}M_\mathrm{H}`$, and the contribution of low masses to the cosmological PBH density is modest. Here, we consider various cases where $`\delta _\mathrm{c}`$ lies within the range $`1/3\delta _\mathrm{c}0.7`$.

Following \[6–8\], we assume that there is a sharp peak in the fluctuation spectrum on some fixed (in comoving coordinates) scale $`\xi =r/a(t)`$. Since the PBH formation threshold is large, $`\delta _\mathrm{H}>1/3`$, most $`\xi `$-scale fluctuations do not collapse into a PBH but are preserved until the passage to the dust stage if the dark matter particle free streaming length is moderately large . The mass $`M_\mathrm{H}`$ within the horizon as a function of the mass $`M_x`$ of the dust component in fluctuations of the same comoving scale can be calculated by using the standard Friedmann equations. Noncompact objects are formed on time scales $`tt_{\mathrm{eq}}`$, where $`t_{\mathrm{eq}}`$ is the time of equality between the matter and radiation densities. At this epoch, a flat cosmological model serves as a good approximation:

$$\begin{array}{c}a(\eta )=a_{\mathrm{eq}}\left[2\frac{\eta }{\eta _{}}+\left(\frac{\eta }{\eta _{}}\right)^2\right];\hfill \\ \\ ct=a_{\mathrm{eq}}\eta _{}\left[\left(\frac{\eta }{\eta _{}}\right)^2+\frac{1}{3}\left(\frac{\eta }{\eta _{}}\right)^3\right],\hfill \end{array}$$
(4)
where $`\eta _{}^2=2\pi G\rho _{\mathrm{eq}}a_{\mathrm{eq}}^2/3c^2`$, $`a_{\mathrm{eq}}`$ is the scale factor at time $`t_{\mathrm{eq}}`$,

$$\rho _{\mathrm{eq}}=\rho _0(1+z_{\mathrm{eq}})^3=3.210^{20}\left(\frac{h}{0.6}\right)^8\left(\frac{\mathrm{\Omega }_\mathrm{m}}{0.3}\right)^4\text{ g sm}^3,$$
(5)
$`1+z_{\mathrm{eq}}=2.3210^4\mathrm{\Omega }_\mathrm{m}h^2`$, $`\rho _0=1.910^{29}\mathrm{\Omega }_\mathrm{m}h^2\text{ g sm}^3`$ is the present cosmological matter density, and $`h`$ is the Hubble constant (in units of 100 km s -1 Mpc -1). We perform our calculations for two cosmological models: the model with the present matter density parameter $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and the cosmological term $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_\mathrm{m}0.7`$ and the model without the $`\mathrm{\Lambda }`$ term and with $`\mathrm{\Omega }_\mathrm{m}=1`$. The presence of the $`\mathrm{\Lambda }`$ term reduces only to a change in $`\rho _{\mathrm{eq}}`$; it does not affect the formation of noncompact objects, because the $`\mathrm{\Lambda }`$ term contributes negligibly to the total cosmological density on time scales $`tt_{\mathrm{eq}}`$.

For $`M_\mathrm{H}`$ and $`M_x`$, we have the expressions

$$M_\mathrm{H}=\frac{4\pi }{3}\rho _\mathrm{H}(a(\eta _\mathrm{H})\xi )^3,M_x=\frac{4\pi }{3}\rho _0(a(\eta _0)\xi )^3.$$
(6)
On the horizon scale, $`a(\eta _\mathrm{H})\xi =2ct_\mathrm{H}`$ with $`\eta _\mathrm{H}\eta _{}`$ and $`\rho _\mathrm{H}=3/32\pi Gt_\mathrm{H}^2`$. The present density is $`\rho _0=\rho _{\mathrm{eq}}(a_{\mathrm{eq}}/a_0)^3`$.

Using (4), we obtain

$$\begin{array}{c}M_\mathrm{H}=\frac{1}{2^{2/3}}\left(\frac{3}{2\pi }\right)^{1/6}\frac{M_x^{2/3}c}{G^{1/2}\rho _{\mathrm{eq}}^{1/6}}=\hfill \\ \\ =1.9610^5\left(\frac{M_x}{0.1M_{}}\right)^{2/3}\left(\frac{\mathrm{\Omega }_\mathrm{m}}{0.3}\right)^{2/3}\left(\frac{h}{0.6}\right)^{4/3}M_{},\hfill \end{array}$$
(7)
$$t_\mathrm{H}=\frac{GM_\mathrm{H}}{c^3}=0.5\left(\frac{M_\mathrm{H}}{10^5M_{}}\right)\text{ s}.$$
(8)
Denote the rms density fluctuation $`\delta _\mathrm{H}`$ by $`\mathrm{\Delta }_\mathrm{H}\delta _\mathrm{H}^2^{1/2}`$. The fraction of the radiation mass that transformed into PBHs at time $`t_\mathrm{H}`$ is then

$$\begin{array}{c}\beta =\underset{\delta _\mathrm{c}}{\overset{1}{}}\frac{d\delta _\mathrm{H}}{\sqrt{2\pi }\mathrm{\Delta }_\mathrm{H}}\mathrm{exp}(\frac{\delta _\mathrm{H}^2}{2\mathrm{\Delta }_\mathrm{H}^2})=\hfill \\ \\ =\frac{1}{2}\left[\mathrm{erf}\left(\frac{1}{\sqrt{2}\mathrm{\Delta }_\mathrm{H}}\right)\mathrm{erf}\left(\frac{\delta _\mathrm{c}}{\sqrt{2}\mathrm{\Delta }_\mathrm{H}}\right)\right]\frac{\mathrm{\Delta }_\mathrm{H}}{\delta _\mathrm{c}\sqrt{2\pi }}\mathrm{exp}(\frac{\delta _\mathrm{c}^2}{2\mathrm{\Delta }_\mathrm{H}^2}),\hfill \end{array}$$
(9)
where $`\mathrm{erf}(x)`$ is the error integral. Since, according to , the fraction of the collapsing mass of the Universe for critical gravitational collapse is $`0.8\beta `$, Eq. (9) for critical collapse is also a good approximation.

Using (7) and (9), we can calculate the cosmological PBH density parameter $`\mathrm{\Omega }_{\mathrm{BH}}`$ at the present time $`t_0`$:

$$\begin{array}{c}\mathrm{\Omega }_{\mathrm{BH}}=\frac{\beta }{2^{2/3}}\left(\frac{3}{2\pi }\right)^{1/6}\frac{c}{M_x^{1/3}G^{1/2}\rho _{\mathrm{eq}}^{1/6}}\hfill \\ \\ \\ 3.4510^5\left(\frac{M_x}{0.1M_{}}\right)^{1/3}\left(\frac{\mathrm{\Omega }_\mathrm{m}}{0.3}\right)^{2/3}\left(\frac{h}{0.6}\right)^{4/3}\frac{\mathrm{\Delta }_\mathrm{H}}{\delta _\mathrm{c}}\mathrm{exp}\left\{\frac{\delta _\mathrm{c}^2}{2\mathrm{\Delta }_\mathrm{H}^2}\right\}.\hfill \end{array}$$
(10)
To within a factor of order unity, expression (7) can be derived from a simple estimate $`M_\mathrm{H}M_xa(t_{\mathrm{eq}})/a(t_\mathrm{H})M_x(t_{\mathrm{eq}}/t_\mathrm{H})^{1/2}`$, where $`t_\mathrm{H}=GM_\mathrm{H}/c^3`$ and $`t_{\mathrm{eq}}610^{10}`$ s is the completion time of the radiation-dominated stage. In the same way, we can derive (10) from the estimate $`\mathrm{\Omega }_{\mathrm{BH}}\beta a(t_{\mathrm{eq}})/a(t_\mathrm{H})`$.

PBHs are formed in the tail of the Gaussian fluctuation distribution, while most noncompact objects are formed from rms fluctuations. Therefore, we repeat that not each fluctuation, by any means, from which a noncompact object formed could collapse into a PBH at time $`t_\mathrm{H}`$.

## 3 THE EVOLUTION OF PERTURBATIONS

The evolution of a radiation density perturbation at the radiation-dominated stage follows the law

$$\delta _\mathrm{r}=xf(x)+\frac{3x^2}{x^2+6}\frac{d}{dx}f(x),$$
(11)
where for the growing mode $`f(x)=A_{\mathrm{in}}j_1(x/\sqrt{3})`$, $`j_1`$ is the spherical Bessel function, $`A_{\mathrm{in}}`$ is the normalization constant, $`x=k\eta `$, and $`k`$ is the comoving perturbation wave vector. The physical perturbation wavelength satisfies the relations

$$\lambda _{\mathrm{ph}}(\eta )=a(\eta )\left(\frac{2\pi }{k}\right);M_x=\frac{4\pi }{3}\rho _0\left(\frac{\lambda _{\mathrm{ph}}(t_0)}{2}\right)^3.$$
(12)
On the horizon scale, $`\lambda _{\mathrm{ph}}/22ct`$, $`x_\mathrm{H}\pi `$, and we obtain from (11) $`\delta _\mathrm{r}=\delta _\mathrm{H}=A_{\mathrm{in}}\varphi `$, where $`\varphi 0.817`$.

For adiabatic perturbations, the perturbation in nonrelativistic matter at $`x1`$ is $`\delta =3\delta _\mathrm{r}/4`$. In , an analytic solution was found for $`\delta `$ in the entire interval from $`x1`$ to $`x1`$ at the radiation-dominated stage. At $`x1`$, this solution is

$$\delta =\frac{3A_{\mathrm{in}}}{2}\left[\mathrm{ln}\left(\frac{x}{\sqrt{3}}\right)+\gamma _\mathrm{E}\frac{1}{2}\right],$$
(13)
where $`\gamma _\mathrm{E}1/20.077`$ and $`A_{\mathrm{in}}=\delta _\mathrm{H}/\varphi `$ is the same as in (11).

The applicability of (13) is limited to a linear stage, $`\delta 1`$. When passing to a nonlinear stage, we will use the spherical model from . In this model, the evolution of adiabatic perturbations is described by the equation

$$y(y+1)\frac{d^2b}{dy^2}+\left[1+\frac{3}{2}y\right]\frac{db}{dy}+\frac{1}{2}\left[\frac{1}{b^2}b\right]=0,$$
(14)
where $`y=a(\eta )/a_{\mathrm{eq}}`$ and the following parametrization was introduced for the radius of the perturbed region:

$$r=a(\eta )b(\eta )\xi ,$$
(15)
Here, $`\xi `$ is the comoving coordinate and $`b(\eta )`$ allows for the deceleration of cosmological expansion in the region of enhanced density. The quantity $`b`$ in (15) can be expressed in terms of $`\delta `$ as

$$b=(1+\delta )^{1/3}.$$
(16)
This relation means the passage from the Eulerian description of the perturbation evolution (13) to its Lagrangian description (14).

To solve (14) requires specifying $`\delta _\mathrm{i}`$ at some initial $`y_\mathrm{i}`$, $`b_\mathrm{i}=(1+\delta _\mathrm{i})^{1/3}`$ according to (16), and the expansion rate $`db/dy`$. In , Eq.(14) was solved for entropic perturbations when the initial velocity may be disregarded, $`db/dy0`$. In our case of adiabatic perturbations, the initial velocity is large; we specify it according to the solution (13). At $`x1`$ and $`y1`$, we have

$$\begin{array}{c}x=\frac{\pi }{2^{2/3}}\left(\frac{3}{2\pi }\right)^{1/6}\frac{yc}{M_x^{1/3}G^{1/2}\rho _{\mathrm{eq}}^{1/6}}=\hfill \\ \\ =2.8610^6y\left(\frac{\mathrm{\Omega }_\mathrm{m}}{0.3}\right)^{2/3}\left(\frac{h}{0.6}\right)^{4/3}\left(\frac{M_x}{M_{}}\right)^{1/3}.\hfill \end{array}$$
(17)
We obtain from (13), (16), and (17)

$$\frac{db}{dy}|_{y_\mathrm{i}}=\frac{\delta _\mathrm{H}b_\mathrm{i}^4}{2y_\mathrm{i}\varphi }.$$
(18)
We solve Eq. (14) numerically. The time $`y_\mathrm{i}`$ must be chosen in such a way that (13) and (14) are simultaneously valid. An optimum choice is the time when $`\delta _\mathrm{i}=0.2`$. A test shows that the results of our calculations change by no more than 15% for a different choice of $`\delta _\mathrm{i}`$ in the range $`0.10.3`$. Having specified $`\delta _\mathrm{i}`$, we obtain $`x_\mathrm{i}`$ and $`y_\mathrm{i}`$ from (13) and (17). The evolution of $`\delta =b^31`$ is illustrated in Fig. 1.

The cosmological expansion of noncompact objects stops when $`dr/dt=0`$. According to , this is equivalent to

$$\frac{db}{dy}=\frac{b}{y}.$$
(19)
Denote $`b`$ and $`y`$ at the time the expansion stops by $`b_{\mathrm{max}}`$ and $`y_{\mathrm{max}}`$, respectively. The CDM density in noncompact objects at the same time is

$$\rho _{\mathrm{max}}=\rho _{\mathrm{eq}}y_{\mathrm{max}}^3b_{\mathrm{max}}^3$$
(20)
and, consequently, the radius of the noncompact object at the stoppage time is

$$R_{\mathrm{max}}=\left(\frac{3M_x}{4\pi \rho _{\mathrm{max}}}\right)^{1/3}.$$
(21)
The numerically calculated dependence of $`R_{\mathrm{max}}`$ on $`\delta _\mathrm{H}`$ is shown in Fig. 2.

Thus, we have shown how the radius of a noncompact object at the time its cosmological expansion stops can be determined for the specified perturbation $`\delta _\mathrm{H}`$ on the horizon scale.

## 4 THE PARAMETERS OF
GRAVITATIONAL LENSES

After its cosmological expansion stops, a noncompact object begins to contract, with its final radius being

$$R_x=\kappa R_{\mathrm{max}},$$
(22)
where $`\kappa `$ is the coefficient of nonlinear contraction. It is generally believed  that after the cosmological expansion of a noncompact object stops, it is virialized, radially contracting by a factor of 2, i.e., $`\kappa =0.5`$. According to the theory of gravitational instability, $`\kappa 0.3`$ in a multiflow region . At present, the physical processes that could cause a noncompact object to contract to $`\kappa <0.3`$ are unknown. Therefore, we take $`\kappa 0.3`$ as the lower limit.

For a noncompact object to be able to serve as a gravitational lens producing observable microlensing events in the Galactic halo, its radius should not significantly exceed the Einstein radius

$$R_\mathrm{E}=2\sqrt{GM_xd/c^2},$$
(23)
where $`d20`$ kpc for microlensing in the halo. We define $`\epsilon `$ as

$$\epsilon 10\frac{R_\mathrm{E}}{R_x}$$
(24)
According to calculations , the inequality $`\epsilon 1`$ must hold. If $`\epsilon <1`$, then the theory comes into conflict with observational data on the light curves . One of the microlensing events with a lens mass $`M_x0.02M_{}`$ was studied in detail in . A comparison of the observed and calculated light curves showed that if a noncompact object has no baryonic core at its center, then it must have the size $`R_x=1.610^{13}`$ cm and a compactness parameter (in our notation) $`\epsilon 19`$. In the presence of a baryonic core with a mass of $`0.05M_x`$, $`R_x=5.710^{13}`$ cm and $`\epsilon 4.8`$. It was noted in  that the model of a point-like lens for this event is also acceptable because of the large observational errors.

We obtain from (21), (23), and (24)

$$\kappa \epsilon =10\frac{R_\mathrm{E}}{R_{\mathrm{max}}}=\frac{1.910^{15}\text{ cm}}{R_{\mathrm{max}}}\left(\frac{M_x}{M_{}}\right)^{1/2}\left(\frac{d}{20\text{ kpc}}\right)^{1/2}$$
(25)
On the other hand, using our calculations (see the preceding section), we derived a relationship between $`R_{\mathrm{max}}`$ and $`\delta _\mathrm{H}`$; the quantity $`\delta _\mathrm{H}`$ defines the present cosmological PBH density according to Eq. (10). It should be noted that noncompact objects are formed from rms fluctuations. Therefore, $`\delta _\mathrm{H}`$ from the preceding section should be set equal to the rms fluctuation $`\delta _\mathrm{H}=\mathrm{\Delta }_\mathrm{H}`$. Using our numerical calculations, we derive the dependence of $`\mathrm{\Omega }_{\mathrm{BH}}`$ on the product $`\kappa \epsilon `$ (see Figs. 3 and 4).

There are several astrophysical constraints on the mass and number of PBHs. It follows from a constraint on the age of the Universe that $`\mathrm{\Omega }_{\mathrm{BH}}1`$. If PBHs provide the dominant part of dark matter in the Galaxy, then they must tidally interact with globular clusters by disrupting them. The PBH mass was constrained for this case in , $`M_{\mathrm{BH}}10^4M_{}`$. At $`\mathrm{\Omega }_{\mathrm{BH}}1`$, PBHs are capable of distorting the CMB spectrum if they are formed about 1 s after the annihilation of $`e^+e^{}`$\- pairs . Mass accretion by black holes at the pre-galactic and present epochs contributes to the background radiation in different wavelength ranges. However, calculations strongly depend on the model and yield $`\mathrm{\Omega }_{\mathrm{BH}}10^3÷10^1`$ for $`M_{\mathrm{BH}}10^5M_{}`$. In , the constraint $`\mathrm{\Omega }_{\mathrm{BH}}<0.1`$ on intergalactic PBHs was obtained from the condition for the absence of reliable gamma-ray-burst lensing events for $`10^5M_{}<M_{\mathrm{BH}}<10^9M_{}`$. A more stringent lensing constraint, $`\mathrm{\Omega }_{\mathrm{BH}}<0.01`$ for the mass range $`10^6M_{}<M_{\mathrm{BH}}<10^8M_{}`$, was obtained from VLBI observations of compact radio sources .

Let us first consider the microlensing event studied in  by assuming that $`\epsilon =0.3`$. If the noncompact object has no baryonic core (the vertical line $`\kappa \epsilon 5.7`$ in Figs. 3 and 4 corresponds to this case), then the constraint $`\mathrm{\Omega }_{\mathrm{BH}}<0.1`$ rules out the interpretation of this event as microlensing by noncompact objects for all $`\delta _\mathrm{c}=1/3÷0.7`$ in the two cosmological models under consideration. In the presence of a baryonic core (the vertical line $`\kappa \epsilon 1.44`$), the case with $`\delta _\mathrm{c}=0.7`$ remains admissible in the cosmological model without the $`\mathrm{\Lambda }`$ term (Fig. 4).

Consider the less stringent condition $`\kappa \epsilon >0.3`$, which is satisfied at $`\epsilon >1`$ and $`\kappa >0.3`$. At $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$ (Fig. 3), the constraint $`\mathrm{\Omega }_{\mathrm{BH}}<0.1`$ rules out all models with $`\delta _\mathrm{c}<0.5`$. At $`\mathrm{\Omega }_{\mathrm{BH}}>10^6`$, the case with $`\delta _\mathrm{c}0.7`$ remain possible. If $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, then the constraint $`\mathrm{\Omega }_{\mathrm{BH}}<0.1`$ leaves a narrow region (see Fig. 4) with $`M_x>0.1M_{}`$ and $`\delta _\mathrm{c}=1/3`$. The constraint $`\mathrm{\Omega }_{\mathrm{BH}}<10^3`$ completely rules out the models with $`\delta _\mathrm{c}=1/3`$, but allows the cases with $`\delta _\mathrm{c}>1/3`$.

## 5 CONCLUSIONS

We have shown that the formation of noncompact dark-matter objects (neutralino stars) proposed in  to account for the observed microlensing events in the Galactic halo must be preceded by the formation of a PBH with a mass of the order of $`10^5M_{}`$. These PBHs and neutralino stars are formed from the same peak in the primordial fluctuation spectrum.

Astrophysical constraints on the number of PBHs in the Universe allowed us to constrain the coefficient of nonlinear contraction and compactness parameter for a noncompact object $`\kappa \epsilon `$; our constraints are shown in Figs. 3 and 4. The most stringent constraints are obtained in the presence of a cosmological term $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$. In this case, noncompact objects can serve as gravitational lenses only at a large PBH formation threshold, $`\delta _\mathrm{c}>0.5`$, which was calculated in the model of critical gravitational collapse. At smaller $`\delta _\mathrm{c}`$, the model of noncompact objects as microlensing objects is ruled out. The constraint on $`\kappa \epsilon `$ is significantly relaxed if the $`\mathrm{\Lambda }`$ term is small (Fig. 4). However, this possibility is currently considered to be unlikely.

It is important to note that to avoid the situation where $`\mathrm{ln}(\mathrm{\Omega }_{\mathrm{BH}})0`$ requires an accurate adjustment of the parameters for the inflationary model that ensures that $`\mathrm{\Delta }_\mathrm{H}`$ is in a narrow interval, $`0.05÷0.12`$. Therefore, if the observed microlensing events are actually produced by noncompact objects, then these objects most likely originate from entropic density perturbations. In this case, our constraints are removed, as discussed in the Introduction.

The constraints can also be significantly relaxed if for some reason, PBHs are formed in smaller quantities than that given by expression (9). This is possible, for example, in the case of great importance of the nonlinear effects that accompany the generation of metric perturbations at the inflationary stage if these effects suppress the appearance of large fluctuations . Note, however, that the role of nonlinear effects is presently not completely understood, and the results of calculations depend on the specific inflationary model. For example, it was found in  that nonlinear effects cause the PBH formation probability to increase, which is directly opposite to the result from .

Conversely, if the noncompact nature of lenses will be proven in the future, then for adiabatic density perturbations, this leads us to conclude that a large number of PBHs with masses of the order of $`10^5M_{}`$ can be formed at the radiation-dominated sage of the Universe. These PBHs must affect the evolution of galaxies and their nuclei. It may well be that one of such PBHs was found by the Chandra space X-ray observatory in the galaxy M 82 . In a separate paper , we develop a model for the formation of galaxies through multiple mergers of protogalaxies with condensation centers in the form of massive PBHs. The mergers of galaxies and the growth of central massive black holes at the galactic nuclei take place simultaneously with the formation of galaxies. The recently found correlations between the masses of central black holes and the bulge velocity dispersion have been shown to naturally arise in this scenario.

Since noncompact objects and PBHs originate from a common perturbation spectrum, we can in principle reconstruct the shape of the perturbation spectrum and determine the PBH mass function from the distribution of noncompact objects in mass and radius using a Press–Schechter-type theory . Unfortunately, only a few objects were detected by microlensing, and such a calculation will become possible only in the future when the statistics will improve.

We may consider a situation that, in a sense, is reverse to the situation described previously. According to the hypothesis , the dark halo objects responsible for microlensing are PBHs with masses of the order of $`(0.01÷1)M_{}`$. If elementary particles with a free streaming length $`\xi `$ comprise the remaining part of the dark matter, then miniclusters will be formed from these particles at the dust stage. Using formula (7), we can immediately estimate the minicluster mass as $`10^{11}÷10^8M_{}`$. Such masses are possible if the mass of the dark-matter particles exceeds 1 GeV . The formation of PBHs with $`\mathrm{\Omega }_{\mathrm{BH}}1`$ requires rms fluctuations $`\mathrm{\Delta }_\mathrm{H}0.06`$. At time $`t_{\mathrm{eq}}`$, fluctuations in the dust component $`\delta 1`$, which give rise to miniclusters, correspond to them.

ACKNOWLEDGMENTS

We thank K.P. Zybin and V.N. Lukash for helpful discussions. This study was supported by the INTAS (grant no. 99-1065) and the Russian Foundation for Basic Research (project nos. 01-02-17829, 00-15-96697, and 00-15-96632).

FIGURE CAPTIONS

Fig. 1. The evolution of a CDM density perturbation $`\delta `$. The plot corresponds to the parameters $`\delta _\mathrm{H}=0.04`$, $`M_x=0.1M_{}`$, and $`\mathrm{\Omega }_\mathrm{m}=0.3`$. The curve was obtained from formula (13) before the point ($`y_\mathrm{i}=3.9610^6`$, $`\delta _\mathrm{i}=0.2`$), and a numerical solution of Eq. (14) was used at $`y>y_\mathrm{i}`$. The cosmological expansion of a noncompact object stops at the radiation-dominated stage at $`y=a/a_{\mathrm{eq}}0.49`$. The dashed line indicates the evolution of $`\delta `$ according to the linear theory (13).

Fig. 2. Radius $`R_{\mathrm{max}}`$ of noncompact objects at the time their expansion stops versus perturbation magnitude $`\delta _\mathrm{H}`$ on the horizon scale. Solid curves 1, 2, and 3 correspond to the masses of noncompact objects $`M_x=1M_{}`$, $`0.1M_{}`$, and $`0.01M_{}`$ in the cosmological model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$. Dashed lines 4, 5, and 6 were constructed for $`\mathrm{\Omega }_\mathrm{m}=1`$ at the same masses.

Fig. 3. PBH density parameter $`\mathrm{\Omega }_{\mathrm{BH}}`$ versus nonlinear contraction coefficient and compactness parameter for noncompact objects $`\kappa \epsilon `$ in the model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$. Curves (1, 2, 3) correspond to the masses of noncompact objects $`M_x=(0.01,0.1,1)M_{}`$ at $`\delta _\mathrm{c}=1/3`$. Curves (4, 5, 6) correspond to $`\delta _\mathrm{c}=1/2`$ and curves (7, 8, 9) correspond to $`\delta _\mathrm{c}=0.7`$ for the same masses. The horizontal and vertical lines correspond to the upper observational limits $`\mathrm{\Omega }_{\mathrm{BH}}=0.1`$, $`\mathrm{\Omega }_{\mathrm{BH}}=10^3`$ and to $`\kappa \epsilon =0.3`$, $`1.44`$, respectively. The admissible regions lie below the horizontal lines and to the right from the vertical lines.

Fig. 4. Same as Fig. 3 for the cosmological model with $`\mathrm{\Omega }_\mathrm{m}=1`$. The vertical lines correspond to $`\kappa \epsilon =0.3`$, $`1.44`$ and $`5.7`$.


text:
# Domain structure of superconducting ferromagnets

\[

## Abstract

In superconducting ferromagnets the equilibrium domain structure is absent in the Meissner state, but appears in the spontaneous vortex phase (the mixed state in zero external magnetic field), though with a period, which can essentially exceed that in normal ferromagnets. Metastable domain walls are possible even in the Meissner state. The domain walls create magnetostatic fields near the sample surface, which can be used for experimental detection of domain walls.

\]

Recently there has been a growing interest to materials, in which superconductivity and ferromagnetism coexist . A number of unusual phenomena and structures have been predicted and observed, spontaneous vortex phase as an example. But the theory mostly addressed macroscopically uniform structures, whereas ferromagnetic materials, even ideally uniform, inevitably have a domain structure, which is a ground-state property of ferromagnets. So a further progress in studying materials with coexisting ferromagnetism and superconductivity requires an analysis of the domain structure. The present work is the first step in this direction.

An object of the study is a material, in which the magnetic transition occurs earlier, i.e. at a higher temperature, than the superconductivity onset. This was called “superconducting ferromagnet” , in contrast to “ferromagnetic (or magnetic) superconductors” where the superconductivity sets in before the magnetic transition, which have been studied mostly in the past . Competition of ferromagnetism and superconductivity may result in various structures with the magnetic moment rotating in space (spiral structures, cryptoferromagnetism and so on). This also can be considered as a “domain structure”, but with a period determined by intrinsic properties of materials. However, our goal is the domain structure due to magnetostatic fields generated by nonzero average bulk magnetization $`\stackrel{}{M}`$. In this case the domain size depends on a sample size. We shall consider type-II superconductivity, bearing in mind ruthenocuprates , which are type II high-$`T_c`$ superconductors.

Before analyzing the domain structure it is useful to summarize the magnetic properties of a single-domain superconducting ferromagnet. The total free energy of the superconducting ferromagnet can be written as

$$F(\stackrel{}{M},\stackrel{}{B})=f_E+K+\frac{(\stackrel{}{B}4\pi \stackrel{}{M})^2}{8\pi }+\frac{2\pi \lambda ^2}{c^2}j_s^2,$$
(1)
where $`\lambda `$ is the London penetration depth and $`\stackrel{}{B}`$ is the magnetic induction. The energy $`f_E(M,\stackrel{}{M})`$ is the exchange energy, which depends on the absolute value of $`M`$ and on gradients of $`\stackrel{}{M}`$. As a rule , in magnetic materials this is the largest energy, which fixes $`M`$. The anisotropy energy $`K(\stackrel{}{M}/M)`$ is smaller and depends on the direction of $`\stackrel{}{M}`$. We shall consider a stripe magnetic structure, which is possible only if $`K`$ essentially exceeds the magnetostatic energy $`M^2`$ . The latter is determined by the magnetic field $`\stackrel{}{H}=\stackrel{}{B}4\pi \stackrel{}{M}`$ \[the third term in Eq. (1)\]. The expression Eq. (1) includes also the kinetic energy related to the superconducting current

$$\stackrel{}{j}_s=\frac{c\mathrm{\Phi }_0}{8\pi ^2\lambda ^2}\left(\stackrel{}{}\phi \frac{2\pi \stackrel{}{A}}{\mathrm{\Phi }_0}\right),$$
(2)
where $`\mathrm{\Phi }_0`$ is the magnetic-flux quantum, $`\phi `$ is the phase of the superconducting order parameter, and the vector potential $`\stackrel{}{A}`$ determines the magnetic induction $`\stackrel{}{B}=\stackrel{}{}\times \stackrel{}{A}`$. The kinetic energy of superconducting currents is absent in a normal ferromagnet.

Minimization of the energy with respect to the vector potential $`\stackrel{}{A}`$ yields the Maxwell equation

$$\frac{4\pi }{c}\stackrel{}{j}_s=\stackrel{}{}\times (\stackrel{}{B}4\pi \stackrel{}{M})=\stackrel{}{}\times \stackrel{}{B}.$$
(3)
Together with the equation

$$\stackrel{}{}\times \stackrel{}{j}_s=\frac{c}{4\pi \lambda ^2}\stackrel{}{B}$$
(4)
this yields the London equation which determines $`\stackrel{}{B}`$:

$$\lambda ^2\stackrel{}{}\times [\stackrel{}{}\times \stackrel{}{B}]+\stackrel{}{B}=0.$$
(5)
Here we took into account that $`\stackrel{}{}\times \stackrel{}{M}=0`$ inside domains. In contrast to Ref. , we neglect the differential susceptibility ($`M`$ does not depend on a magnetic field), which renormalizes the London penetration depth.

These equations and the boundary conditions at the sample boundary (continuity of the tangential component of $`\stackrel{}{H}`$ and of the normal component of $`\stackrel{}{B}`$) yield the distribution of $`\stackrel{}{B}`$ and $`\stackrel{}{H}`$. This distribution is shown in Fig. 1 for the case of $`\stackrel{}{M}`$ parallel to the sample boundary and for zero external magnetic field. The magnetic induction and the related magnetic flux exist only in the layer of the thickness $`\lambda `$. Meissner currents in this layer screen the internal field $`4\pi \stackrel{}{M}`$, as well as they screen the external magnetic field in a nonmagnetic superconductor.

Let us consider now the mixed state of the superconducting ferromagnetic, in which vortices (magnetic fluxons) are present in the bulk. Since ferromagnetism does not affect the London equation (5), one expect the same magnetic-induction distribution in the mixed state as for nonmagnetic type II superconductors , and the free energy is given by

$$F_m(\stackrel{}{M},\stackrel{}{B})=f_E+K+2\pi M^2\stackrel{}{B}\stackrel{}{M}+F_0(B),$$
(6)
where $`F_0(B)`$ is the free energy of a nonmagnetic type II superconductor, and $`\stackrel{}{B}`$ now is the magnetic induction averaged over the vortex-array cell. The energy $`F_0(B)`$ contains both the magnetic energy $`B^2/8\pi `$ and the kinetic energy of the superconducting currents inside the vortex cell. Determining the magnetic field

$$\stackrel{}{H}=4\pi \frac{F_m}{\stackrel{}{B}}=4\pi \frac{F_0}{\stackrel{}{B}}4\pi \stackrel{}{M},$$
(7)
we see that the magnetization curve of a superconducting ferromagnet is described by $`B=B_0(|\stackrel{}{H}+4\pi \stackrel{}{M}|)`$ where $`B_0(H)`$ is the equilibrium magnetization curve for a nonmagnetic type II superconductor  (Fig. 2a). Note that in this relation the magnetic field $`\stackrel{}{H}`$ has a different physical meaning from that used in the Meissner state. For the Meissner state we introduced $`\stackrel{}{H}=\stackrel{}{B}4\pi \stackrel{}{M}`$, where the moment $`\stackrel{}{M}`$ originates from “molecular” currents responsible for ferromagnetism, the superconducting currents being treated as external currents. In the mixed state, which is considered now, it is more convenient to define the magnetic field as $`\stackrel{}{H}=\stackrel{}{B}4\pi (\stackrel{}{M}+\stackrel{}{M}_s)`$, i.e. the definition includes also the diamagnetic moment $`\stackrel{}{M}_s=(\stackrel{}{B}_0\stackrel{}{H})/4\pi `$ of the superconducting currents circulating around vortex lines in the mixed state. Thus these currents are treated in the same manner as molecular currents responsible for ferromagnetism.

Figure 2b shows that in a superconducting ferromagnet the Meissner state ($`B=0`$) exists until $`H+4\pi M<H_{c1}`$, where $`H_{c1}=(\mathrm{\Phi }_0/\lambda ^2)\mathrm{ln}(\lambda /\xi )`$ is the lower critical field in a nonmagnetic superconductor and $`\xi `$ is the coherence length, which determines the vortex core size. So ferromagnetism decreases the lower critical field $`\stackrel{~}{H}_{c1}=H_{c1}4\pi M`$. If $`4\pi M>H_{c1}`$, the Meissner state is absent (Fig. 2c) and the superconducting ferromagnet is in the mixed state with vortices penetrating into it even in zero external field $`H=0`$. This is spontaneous vortex phase with nonzero magnetic induction $`B=B_0(4\pi M)`$ in the bulk.

Now let us consider formation of the domain structure in the standard geometry : a slab of the thickness $`d`$ along the anisotropy easy axis $`y`$ and infinite in directions of the axes $`x`$ and $`z`$ (Fig. 3). We start from a normal ferromagnet. In the absence of an external magnetic field the average magnetic induction inside the slab must vanish. Therefore, $`B=0`$ in a single-domain structure (Fig. 3a), and there exists an uniform magnetostatic field $`\stackrel{}{H}=4\pi \stackrel{}{M}`$ in the entire sample, an analog of the electrostatic field in a charged plane capacitor. This results in a high magnetostatic energy $`H^2/8\pi M^2`$. However, the domain structure with period $`l`$ ($`ld`$) suppresses this energy in the domain bulk: $`\stackrel{}{H}0`$ and $`\stackrel{}{B}=4\pi \stackrel{}{M}`$, except for the area $`l^2`$ near the sample boundary (Fig. 3b). But the average induction still vanishes, since $`\stackrel{}{M}`$ changes its sign from a domain to a domain. For the stripe structure one can solve the equations of magnetostatics, $`\stackrel{}{}\times \stackrel{}{H}=0`$ and $`\stackrel{}{}\stackrel{}{H}=4\pi \rho _M`$, exactly . Here $`\rho _M=\stackrel{}{}\stackrel{}{M}`$ is the magnetic charge. The magnetostatic energy per unit volume of the slab is

$$E_s=0.852M^2l^2\times \frac{1}{ld}=0.852M^2\frac{l}{d}.$$
(8)
This energy is by a factor $`l/d`$ less than the magnetostatic energy in a single-domain structure. However, the domain walls increase the energy. The energy of one domain wall (per unit length along the slab) is $`\alpha K\delta d`$, where $`\delta `$ is the wall thickness and the numerical factor $`\alpha `$ depends on the detailed definition of $`K`$ and $`\delta `$. Its specification is not essential for the present analysis. The domain-wall energy per unit volume of the sample is

$$E_w=\alpha Kd\delta \times \frac{1}{ld}=\alpha K\frac{\delta }{l}.$$
(9)
The equilibrium value of the period $`l`$ is determined by minimization of the energy $`E_w+E_s`$ :

$$l=\sqrt{\frac{\alpha K}{0.852M^2}\delta d}.$$
(10)
Let us return back to a superconducting ferromagnet. In the Meissner state the magnetic induction must vanish in the bulk, which is compatible only with the single-domain structure. Thus the equilibrium domain structure is impossible in the Meissner state. However, domains with the changing direction of $`\stackrel{}{M}`$ can appear in the spontaneous vortex phase with nonzero $`B=B_0(4\pi M)`$. Like in a normal ferromagnet, the magnetic flux $`B`$ in domains should produce the magnetostatic fields in the area $`l^2`$, but in a superconducting ferromagnet these fields are by the factor $`B_0(4\pi M)/4\pi M`$ smaller. We can take it into account introducing the effective magnetization $`\stackrel{~}{M}=B_0(4\pi M)/4\pi `$. Then the period of the domain structure is given by Eq. (10), where $`M`$ must be replaced with $`\stackrel{~}{M}`$. In the limit of large $`4\pi MH_{c1}`$, one has $`\stackrel{~}{M}M`$ and the effect of superconductivity on the domain structure vanishes. In the opposite limit of small $`M`$, when $`4\pi MH_{c1}`$, $`\stackrel{~}{M}`$ vanishes and the period $`l`$ becomes infinite, as it should be in the Meissner state $`4\pi M<H_{c1}`$. However, this calculation of $`l`$ assumes that the penetration of the magnetostatic field into a superconducting ferromagnet is similar to the penetration into a normal ferromagnet. The assumption is correct if rigidity of the vortex array is negligible and the effective penetration depth is infinite. We can also consider the opposite limit of a very rigid vortex array, when the magnetostatic fields penetrate only into the layer of the thickness $`\lambda `$. If $`\lambda l`$, the penetration of the magnetic flux into a superconductor becomes insignificant. This increases the magnetic fields outside the sample, as well as the total magnetostatic energy, by a factor of 2, while the correspondingly period $`l`$ decreases by a factor of $`\sqrt{2}`$ (Fig. 3c), in analogy with the effect of a superconducting substrate on a domain size in a ferromagnetic slab . Thus avoiding a detailed analysis of the vortex and field pattern in the domains close to the sample border we lose only a numerical factor of not more than $`\sqrt{2}`$. In any case, superconductivity, which coexists with ferromagnetism in the same bulk, always increases the domain size, in contrast to the superconductor-ferromagnet bilayer, where superconductivity shrinks ferromagnetic domains .

The absence of the equilibrium domain structure in the Meissner state does not rule out a possibility of metastable domain walls, as topologically stable planar defects. Domains can appear also because of disorder, or grain structure. The structure of the domain wall should be found by solution of the coupled equations of magnetostatics and the London electrodynamics. We restrict ourselves to the simplest case, when the London penetration depth $`\lambda `$ essentially exceeds the domain wall thickness $`\delta `$. This means that at the spatial scales of order $`\delta `$ the domain-wall structure is governed by large energies \[the exchange energy and the anisotropy energy, see Eq. (1)\] and is not affected by the magnetostatic and kinetic energy. On the other hand, at scales $`\lambda `$ one can find the distribution of $`\stackrel{}{B}`$ and $`\stackrel{}{H}`$ from the London equation at constant $`\stackrel{}{M}`$. This is shown for the Bloch domain wall (the magnetization $`\stackrel{}{M}`$ rotates in the plane of the wall and does not produce the magnetostatic charges) in Fig. 4a. Though our picture corresponds to a 180 wall, a similar picture is expected for any domain wall. The jump of the tangential component of the moment $`\stackrel{}{M}`$ at the wall defines the current sheet, responsible for a jump of the magnetic induction parallel to the wall, whereas a possible jump of the normal component of $`\stackrel{}{M}`$ (a “charged” domain wall) would produce a jump of the normal component of the field $`\stackrel{}{H}`$.

The magnetic flux on the opposite sides from the domain wall creates the magnetostatic fields outside the sample, where the wall meets the sample boundary (Fig. 4b). The magnetic fluxes, which exit from the sample at two sides from the wall, are equal in magnitude ($`4\pi M\lambda `$ per unit length along the wall) but opposite in direction. The magnetostatic field from domain walls could be used for their experimental detection. At distances $`r\lambda `$ from a line, where the wall exits to the sample boundary, this field is a dipole field of the order of $`M\lambda ^2/r^2`$.

In summary, the letter presents the first analysis of the domain structure in superconducting ferromagnets. There is no equilibrium domain structure in a superconducting ferromagnet in the Meissner state. In the spontaneous vortex phase the period of the domain structure may essentially exceed that in the normal ferromagnet. But metastable domain walls can exist even in the Meissner state. They generate the magnetic flux in layers of a thickness $`\lambda `$, which can be revealed by magneto-optical methodic.

The author thanks Yu. Barash, I. Felner, and N. Kopnin for useful discussions. The work has been supported by the grant of the Israel Academy of Sciences and Humanities.