text:
# Gravitational waves from newly born, hot neutron stars

## 1 Introduction

According to the standard theory of star formation, neutron stars are born in the aftermath of successful core–collapse supernova explosions, as the stellar remnant becomes gravitationally decoupled from the expanding ejecta. Following the gravitational collapse, the proto–neutron star (PNS) radiates its binding energy (about 0.1 $`M_{}`$) via neutrino emission in a timescale of tens of seconds. A small fraction of this energy reservoir is radiated in gravitational waves (GWs), and it is interesting to establish for how long, and at which frequencies the gravitational radiation will be emitted. The best currently available 2D-numerical simulations of a core collapse, and the most reliable estimates of gravitational wave emission consider, at best, the first hundred milliseconds from the collapse onset. In a number of axisymmetric simulations \[Zwerger & Müller 1997, Dimmelmeier, Font & Müller 2002, Freyer, Holz & Hughes 2002\], the flux of gravitational waves is evaluated through the quadrupole formalism, by using the time variation of the energy density of the collapsing system obtained either with Newtonian or general relativistic hydrodynamical codes. The emitted energy is estimated to range between $`10^8`$ and $`10^6M_{}`$. Full three dimensional simulations in general relativity are not yet available, but there are some Newtonian studies \[Rampp, Müller & Ruffert 1998, Brown 2001\] that, by assuming small deviations with respect to axial symmetry, show similar results. In any event, the compact object which forms after the collapse possesses a certain amount of mechanical energy, which will partially be released through violent oscillations in the form of gravitational radiation.

The spectrum of the non–radial oscillations of a neutron star (NS) has a very rich structure that depends on how the temperature and the internal composition vary along its history. For instance, the frequency of the fundamental mode has been found to scale approximately as $`\sqrt{M/R^3}`$, where $`M`$ is the gravitational mass of the star and $`R`$ is its radius, whereas that of the acoustic p–modes scales linearly with $`M/R`$. Since during the first tens of seconds of its life the star contracts and its gravitational mass decreases, one expects significant changes in the frequency of these modes. In addition, new families of modes may appear or disappear for several different reasons: i) as the star’s entropy gradient decreases with time, the frequency of the gravity modes (g–modes) that are present at finite temperature migrate toward zero; ii) if no significant diffusion mechanism is present, layers of different chemical composition accumulating on the stellar surface would produce one or more density discontinuities, which would give rise to low frequency (typically $`<200`$ Hz) discontinuity g–modes; iii) if the density of the inner core becomes high enough, a phase transition may occur, involving pion and/or kaon condensation as well as a transition from ordinary nuclear matter to quark matter \[Prakash et al. and references therein 1997\] which may produce a density discontinuity in the inner core of the star. In this case, a new discontinuity g–mode would appear, at a frequency higher than the previous ones - typically $`>700`$ Hz - \[Sotani, Tominaga & Maeda 2002, Miniutti et al. 2002\]; iv) in previous works \[Reisenegger & Goldreich 1992, Lai 1994\] it has been shown that g-modes due to composition gradients in the inner core of cold NSs have typical frequency below $`200`$Hz. In hot PNS, the composition gradients are larger and produce g-modes of higher frequencies, even at constant entropy. v) Finally, it should be mentioned that relativistic stars do also possess a family of modes, called w–modes, that are associated to perturbations of the gravitational field, and couple very little (or do not couple at all, as for the axial w–modes) with the motion of the fluid \[Chandrasekhar & Ferrari 1991, Kokkotas & Schutz 1992\]. The properties of the w–modes depend on the internal composition of the star \[Andersson & Kokkotas 1996, Benhar, Berti & Ferrari 1999\], and therefore they also change as the star evolves.

In this paper we study how the characteristic frequencies and damping times of some relevant classes of quasi–normal modes (QNMs) of a newly born, hot PNS change during its early evolution until it settles as a cold, catalyzed NS. Despite the fact that PNSs are expected to be born with a significant amount of angular momentum, in this work we neglect the effects of rotation. In this way, we can isolate and better understand the effects of the thermal and chemical evolution, and we can also clearly indicate the differences between hot PNS and cold NS oscillation properties. This study gives an interesting indication of the frequency range where the gravitational signals emitted immediately after gravitational collapse should be looked for. The paper is organized as follows. In Section 2 we review the evolutionary stages immediately following the birth of a proto–neutron star and discuss the dissipative precesses and timescales relevant during the Kelvin-Helmholtz epoch. In Section 3, we remind the physical processes that produce the gravity modes, and shortly describe the procedure we use to compute the QNM frequencies. The results of our study are presented in Section 4, where we also discuss the detectability of the gravitational signal associated to the quasi–normal modes by resonant bars or ground based interferometric detectors. The main findings of the paper are summarized in the concluding remarks.

## 2 The first minute of a neutron star life.

A newly born PNS is very different from its cold descendant, i.e. a NS. Indeed, initially, the PNS is hot, lepton rich, and optically thick to neutrinos, which are temporarily trapped within the star. The subsequent evolution is dominated by neutrino diffusion which first deleptonizes and later cools the star \[Burrows & Lattimer 1986\]. The evolution of a PNS proceeds through several distinct stages and has been recently studied in great detail, including the possibility of developing a core of exotic matter such as hyperons, a kaon condensate, or quark matter. The results we present in this paper are based on simulations of PNS evolution described in Pons et al. (1999; 2000; 2001). We consider two models of PNS evolution, that use two different equations of state (EOSs). In both cases the EOS of baryonic matter is a finite-temperature, field-theoretical model solved at the mean field level, where the interactions between baryons are mediated by the exchange of $`\sigma `$, $`\omega `$ and $`\rho `$ mesons. Electrons and muons are included in the model as non interacting particles, being the contribution due to their interactions much smaller than that of the free Fermi gas. We consider two of the models of evolving PNS studied by Pons et al. (1999; 2001) where neutrino transport is treated using the diffusion approximation, that allows to obtain global properties of the neutrino emission, like average energies, integrated fluxes and timescales. Model A has a soft EOS \[Glendenning & Moszkowski 1991\] while model B is stiffer \[Müller & Serot 1996\] and allows for the formation of a quark core at some point of the evolution. The mixed phase of baryons and quarks is constructed to satisfy Gibbs’ phase rules for mechanical, chemical and thermal equilibrium, assuming global charge neutrality. The quark matter EOS is computed using the MIT bag model (with Bag constant B= 150 MeV/fm3). A detailed discussion on the EOS, on the evolution of PNSs and on the microphysical inputs (opacities) can be found in the above references. However, here it is useful to summarize the different stages of the early evolution of a PNS referring, as an illustrative example, to Model A. In Figure 1 we show some relevant physical variables that characterize the PNS status at different instants of time ($`0`$, $`0.1`$, $`1`$, $`5`$, $`10`$, $`40`$ seconds): the entropy per baryon, $`s`$, the temperature $`T`$, the electron neutrino chemical potential $`\mu _\nu `$, and the radius $`R`$. Each quantity is plotted versus the enclosed baryonic mass, $`M_B`$.

The different stages of the early evolution can be summarized as follows.

i) Immediately after the core bounce and the passage of the shock wave through the outer PNS mantle, the star contains an unshocked, low entropy core of mass $`0.7`$ M in which neutrinos are trapped. The core is surrounded by a low density, high entropy mantle (see Figure 1, left–bottom) that is accreting matter falling through the shock from the outer iron shell and it is rapidly losing energy due to electron captures and thermal neutrino emission. The mantle extends up to the shock, which is temporarily stalled at a radius of about $`200`$ km prior to an eventual explosion. If the supernova is successful, after a few hundred milliseconds accretion becomes unimportant and the shock lifts off the stellar envelope. In a few tens of seconds, extensive neutrino losses reduce the lepton pressure provoking the contraction of the mantle. The radius of the PNS is now about $`2030`$ km.

ii) Neutrino diffusion deleptonizes the core on a time scale of $`1015`$ s. The diffusion of degenerate, high-energy ($`200`$ MeV) neutrinos from the core to the surface, from where they escape as low-energy neutrinos ($`1020`$ MeV), generates a large amount of heat within the star, resulting in temperatures up to $`50`$ MeV, while the core entropy approximately doubles (see Figure 1, left panels). The inverse temperature gradient which arises in the inner region produces a heat wave that propagates inward. At the end of the deleptonization epoch, the star reaches its maximum central temperature.

iii) After approximately $`15`$ s, the PNS has become lepton-poor but it is still hot. The net number of neutrinos in its interior is low, but thermally produced neutrino pairs of all flavors are abundant and dominate the emission. Thus, neutrino diffusion continues to cool the star, while the average neutrino energy decreases, and the neutrino mean free path increases. After approximately $`50`$ seconds, the mean free path becomes comparable to the stellar radius, and the star finally becomes transparent to neutrinos. By this time, the temperature has dropped to $`15`$ MeV ($`10^{10}`$ K) and the star has radiated almost all of its binding energy. A neutron star is born.

### 2.1 Dissipation time scales

Since the focus of this paper is on the emission of gravitational radiation produced by the violent pulsations of a newly born PNS, it is crucial to understand whether these oscillations may be damped by other dissipative processes acting on timescales shorter than the typical gravitational timescales.

We have seen that in a PNS the energy is transferred and redistributed dominantly by neutrino diffusion. Hence, when one considers hydrodynamical processes, the relevant dissipation time-scales are those corresponding to the neutrino viscosity, diffusivity, thermal conductivity, or thermodiffusion, rather than those associated to the matter kinetic coefficients. Transport properties of degenerate neutrinos in dense matter were studied in detail by Goodwin & Pethick (1982) and by van den Horn & van Weert (1984). They found that for the thermodynamical conditions typical in PNSs and for strongly degenerate neutrinos ($`\mu _\nu k_BT`$) the dominant dissipative process is neutrino chemical diffusion, followed by thermal conduction and viscosity. Thus, the hierarchy of dissipation timescales is $`\tau _D<\tau _C<\tau _{visc},`$ where $`\tau _D`$, $`\tau _C`$, and $`\tau _{visc}`$ are the diffusion, thermal conduction, and viscous timescales, respectively. $`\tau _D`$ can be estimated by the following approximate relations \[van den Horn & van Weert 1984\]

$$\tau _D(\mathrm{s})24sn_B^{1/3}\left(\frac{0.01}{Y_\nu }\right)\left(\frac{T}{10\mathrm{M}\mathrm{e}\mathrm{V}}\right)\left(\frac{\mu _\nu }{100\mathrm{M}\mathrm{e}\mathrm{V}}\right),$$
(1)
where $`s`$ is the entropy per baryon in units of the Boltzmann constant, $`n_B`$ is the baryon density (in fm-3), and $`Y_\nu `$ is the neutrino fraction.

Conversely, if neutrinos are non-degenerate ($`\mu _\nu k_BT`$), the hierarchy between neutrino diffusion and thermal conduction is reversed and $`\tau _C<\tau _D<\tau _{visc}`$. In this case $`\tau _C`$ can be expressed as \[van den Horn & van Weert 1984\]

$$\tau _C(\mathrm{s})1.5sn_B^{1/3}\left(\frac{0.01}{Y_\nu }\right)\left(\frac{T}{10\mathrm{M}\mathrm{e}\mathrm{V}}\right)^2,$$
(2)
In both cases, the timescale $`\tau _{visc}`$ is a factor $`1050`$ longer.

In summary, during the first $`1020`$ seconds of the star life, diffusion is the dominant dissipative process in the inner regions (where neutrinos are degenerate), while thermal conduction dominates in the outer region. Afterwards, thermal conduction becomes the more efficient mechanism of energy dissipation until the matter composing the star becomes transparent to neutrinos. The order of magnitude of both $`\tau _D`$ and $`\tau _C`$ varies from a few seconds to tens of seconds, and it is comparable to the overall evolutionary timescale. This is shown in Figure 2, where the diffusion (solid line) and the conductive (dashed line) timescales are plotted as a function of the radius for model A, and $`t=0.25`$ seconds after the birth of the PNS. The dotted line is the degeneracy parameter of neutrinos ($`\eta =\mu _\nu /T`$), and it is shown to indicate in which region $`\tau _D`$ and $`\tau _C`$ dominate. As a rule of thumb, $`\tau _D`$ must be used when $`\eta >1,`$ while $`\tau _C`$ is the relevant timescale when $`\eta <1`$. In conclusion, a conservative estimate of the dissipative timescales before neutrino transparency is $`\tau _D\tau _C1020`$ seconds. In the following we will denote by $`\tau _{diss}`$ the minimum of the typical dissipative timescales, $`\tau _{diss}=\mathrm{min}(\tau _D,\tau _C,\tau _{visc})`$.

## 3 Quasi-normal modes of relativistic stars

In order to understand what are the differences between the oscillation properties of PNSs and those of old, cold NSs, we shall briefly recollect some results obtained in the literature for cold NSs. A systematic analysis of the pulsation properties of cold NSs has been done by Andersson and Kokkotas (1998) for a large range of EOSs; they found a set of empirical relations that connect the mode frequencies and damping times to the mass and radius of the star. As an example, from their relations one finds that for a cold NS with $`M=1.4M_{}`$ and $`R=10`$ km the inferred values of the frequency and damping times of the f–mode, and of the first p– and w– modes for $`l=2`$ are

$`\nu _f`$ $``$ $`2\mathrm{k}Hz,\nu _{p_1}7\mathrm{k}Hz,\nu _{w_1}11.8\mathrm{k}Hz,`$
$`\tau _f`$ $``$ $`10^1\mathrm{s},\tau _{p_1}5\mathrm{s},\tau _{w_1}2.4\times 10^5\mathrm{s}.`$ (3)
A further family of modes which may appear in NSs are the so-called gravity or g–modes, whose restoring force is buoyancy. It is known \[e.g. Thorne 1966\] that the radial acceleration of a fluid element displaced by its equilibrium position can be written as

$$a=\frac{e^{\lambda /2}}{(ϵ+p)\gamma p}\left(\frac{dp}{dr}\right)S(r)\mathrm{\Delta }r,$$
(4)
where $`ϵ`$ and $`p`$ are the energy density and the pressure, $`\mathrm{\Delta }r`$ is the radial displacement of the fluid element, $`\gamma =[(ϵ+p)/p][p/ϵ]_{s,Y_L}`$ is the adiabatic index, with $`Y_L`$ being the lepton fraction, $`e^\lambda `$ is the $`g_{rr}`$ component of the unperturbed metric tensor, and

$$S(r)=\frac{dp}{dr}\frac{\gamma p}{ϵ+p}\frac{dϵ}{dr}$$
(5)
is the relativistic Schwarzschild discriminant. If $`S(r)>0`$, the displaced fluid element feels a force that tends to restore the initial equilibrium position and the star (or the region inside the star) is stable against convection. In this case, a non–zero frequency spectrum of stable g–modes appear. Conversely, if $`S(r)<0`$, the force tends to move the fluid element away from equilibrium and the star becomes unstable against convection. Finally, if $`S(r)=0`$, all g–modes are degenerate at zero frequency.

By introducing the adiabatic sound velocity, $`c_s^2`$, which is defined at constant entropy $`s`$ and chemical composition, and the so–called equilibrium velocity, $`c_0^2`$,

$$c_s^2=\left(\frac{p}{ϵ}\right)_{\mathrm{s},\mathrm{Y}_\mathrm{L}}=\frac{p}{ϵ+p}\gamma ,c_0^2=\frac{dp/dr}{dϵ/dr},$$
(6)
the Schwarzschild discriminant (5) can be cast into a form which will be useful in following discussions

$$S(r)=\frac{dp}{dr}\left(\mathrm{\hspace{0.17em}1}\frac{c_s^2}{c_0^2}\right).$$
(7)
Thus, a stable g–mode arises only if $`c_s^2>c_0^2`$, and the buoyancy which provides the restoring force can be induced either by entropy and/or composition gradients \[McDermott, van Horn & Hansen 1988, Reisenegger & Goldreich 1992, Lai 1994\], by density discontinuities in the outer envelopes of NSs \[Finn 1987, Strohmayer 1993\], or by density discontinuities in the inner core as a consequence of phase transitions at high density \[Sotani, Tominaga & Maeda 2002, Miniutti et al. 2002\]. The frequency of the g–modes is always lower than that of the f–mode, it is roughly proportional to the Schwarzschild discriminant, and depends on the physical process which is source of the buoyancy. In general, the g–modes of cold NSs have a typical frequency of less than 200-300 Hz and very long damping times (typically, $`\tau _g>10^6`$ s). We shall see that the situation is different during the first second of life of a PNS.

### 3.1 Computation of the modes

In this subsection we clarify some technical aspects of the computational procedure we use to find the frequencies and damping times of the quasi–normal modes. We integrate the equations describing the polar, non radial perturbations of a non rotating star in general relativity. We use the same formulation of the perturbed equations as in Miniutti et al. (2002), based on the work of Detweiler & Lindblom (1985), who derived a fourth order system of equations for the metric and fluid perturbations inside the star. However, we do not use their original method, because we integrate the system in the complex frequency domain directly, rather than working in the real frequency domain. The system of perturbed equations is closed by the equation of state. In particular, the equilibrium profiles of the energy density and of the sound velocity are sufficient to completely solve the perturbation equations in the stellar interior. The profiles we use are the same as those obtained in Pons et al. (1999; 2001).

The set of equations, that couple the perturbations of the fluid with those of the gravitational field, are Fourier expanded and separated by expanding the perturbed tensors in tensorial spherical harmonics. For each assigned value of the harmonic index $`l`$ and of the (complex) frequency $`\nu `$, the system of equations admits only two linearly independent solutions regular at $`r=0`$. The general solution can be found as a linear combination of the two, such that the Lagrangian perturbation of the pressure, $`\mathrm{\Delta }p`$, vanishes at the surface. Outside the star, the fluid perturbations vanish and the system reduces to the Zerilli equation \[Zerilli 1970\]. A quasi normal mode of the star is defined to be a solution of the perturbed equations belonging to a complex eigenfrequency, which is regular at the center, continuous at the surface, and which behaves as a pure outgoing wave at infinity. To compute the QNM complex eigenfrequencies we follow the method developed by Leins, Nollert and Soffel (1993) who derived a procedure (the continued fractions method) to evaluate a function which is proportional to the ingoing part of the wave, and therefore vanishes when $`\nu `$ is an eigenfrequency of a QNM. A very clear account of the method can be found in Sotani, Tominaga & Maeda (2002).

## 4 The gravitational wave spectrum of a proto-neutron star.

Using the method described in Section 3, we have computed the quasi–normal mode frequencies and damping times of the two models of proto–neutron star, at different instant of time from $`t=0.2`$s to $`t=50`$s after their formation. The configurations we consider as initial refers to $`t=0.2`$s for both models, because at earlier times the PNS evolution could be affected by dynamical processes like accretion, and cannot be accurately approximated by a quasi–stationary series of solutions of the hydrostatic equations. Thus we exclude from our study the collapse, the bounce and the first two hundred milliseconds after formation, for which more detailed hydrodynamical simulations are needed.

At $`t=0.2`$s, model A has $`M=1.58M_{}`$ and $`R=34.2`$km, while model B has the same mass and $`R=35.1`$km. Both models are chosen such that the final stage of the PNS evolution, at $`t=50`$s in our simulations, represents a “canonical” NS with a mass of about $`M=1.46M_{}`$ and a radius of $`R1213`$ km. The gravitational mass difference between the initial and final configurations ($`0.12M_{}`$) has been taken away by neutrinos during the PNS evolution. In the following, we shall discuss the behavior of the frequencies of the fundamental mode, and the first g–, p– and w–modes for $`l=2`$, as the star evolves and cools down. To identify the different classes of fluid modes we study the behavior of the Lagrangian radial displacement. The eigenfunction of the f–mode has no radial nodes, while the g1– and p1–modes have one radial node. The w1–mode is easily identified, because it produces negligible motion in the fluid, and because its damping time is much shorter than that of the fluid modes.

The QNM frequencies and damping times for model A and B are tabulated in Tables 1 and 2 at different values of the evolution time. The oscillation properties within the first $`5`$ seconds of the PNS evolution are summarized in Figure 3 for both models.

We find that during the first seconds of life of the PNS the frequencies of the f–, p1– and w1– modes are much lower than those of the cold neutron star which results at the end of the evolution. Initially, they cluster in a narrow range within $`[9001500]`$ Hz, and if we compare the value that a given mode frequency has at $`t=0.2`$s and at $`t=50`$s, we find a difference of about $`40\%`$ for the f–mode, $`75\%`$ for the p1–mode and of $`65\%`$ for the w1–mode. This significant difference can be attributed partly to the fact that at earlier times the PNS is less compact, and partly to the very reach thermodynamical structure that characterizes the first few seconds of the PNS life (see Figure 1). In this respect, the behavior of the the f–mode frequency is of particular interest, since it does not scale with $`\sqrt{M/R^3}`$ as it does for cold NSs: indeed, Tables 1 and 2 show that, during the first second, the mass of the star remains approximately constant while the radius rapidly decreases, and consequently $`\nu _f`$ should increase; on the contrary, it initially decreases, at about $`t0.5`$s approaches the value of the g–mode frequency, and only after about one second, it starts to increase linearly with $`\sqrt{M/R^3}`$. At the same time, the g1–mode frequency increases when the star is younger and hotter and large thermal and composition gradients dominate in the interior, it reaches a maximum value of about $`800900`$ Hz within the first second of life, and afterwards, as the star cools down and entropy gradients are smeared out, it migrates towards lower values. Even though there is a point of the evolution when $`\nu _f`$ and $`\nu _{g_1}`$ are very close (see Figure 3), the two modes are always distinguishable because of the different number of nodes in their eigenfunctions.

From Tables 1 and 2, we see that after a few seconds from birth, the damping times of the f, p1 and w1–modes of the PNS are quite close to those of its cold descendant. At very early times ($`0.2<t<2`$ s) the situation is different; for instance, $`\tau _f`$ decreases rapidly from $`70`$ s to $`23`$ s, and it is interesting to compare its values with $`\tau _{diss}`$, which is about 10-20 seconds, as explained in Section 2.1. For $`t\stackrel{<}{}0.5`$ seconds, $`\tau _{diss}<\tau _f`$, therefore during this time, a small fraction of the amount of energy initially stored in the f–mode, will be dissipated by neutrino processes. However, for $`t\stackrel{>}{}0.5`$ seconds, $`\tau _f`$ becomes smaller than $`\tau _{diss}`$, and gravitational wave emission becomes the dominant dissipation mechanism. In addition, apart from the first few tens of milliseconds, the p1–mode damping time is always comparable to $`\tau _{diss}`$, while $`\tau _{w_1}`$ is basically unaffected by the evolution and it is by far the shortest one. Therefore, if some energy is initially stored into these modes, it will be emitted in gravitational waves.

Furthermore, we find that during the first second, the g1–mode damping time is much shorter than that expected from previous calculations for cold NSs. In particular, for both PNS models $`\tau _{g_1}`$ is comparable to $`\tau _f,`$ to $`\tau _{p_1}`$ and also to the neutrino dissipation timescale during the first half a second. Thus, the common belief, based on cold NSs dynamical behavior, that the damping time of the f–mode is an order of magnitude shorter than that of the first p–mode and many orders of magnitude smaller than the g–modes damping times, is no longer correct during the first second of life, which is likely to be the most important for the emission of GWs.

For g–modes of order higher than one, we find damping times much larger than that of the g1–mode; For instance, the g2–mode has a damping time of the order of $`10^4`$s at $`t=0.4`$s. Thus, g–mode pulsations of order higher than 1 are unlikely to be significant sources of gravitational radiation and will be neglected in our analysis.

From Tables 1 and 2 we see that, for $`t\stackrel{>}{}1020`$seconds, the g1–mode frequency of both models increases. As explained in Section 3, g–modes can be associated to several physical processes, and during the early evolution of a PNS they arise mainly because there are large entropy and composition gradients. While during the deleptonization epoch entropy gradients dominate the dynamical behavior of the star, after some time composition gradients become important, and this change produces a modulation in the frequency of the g–modes. In addition, at $`t=50`$ s we see that $`\nu _{g_1}`$ of model A and B are different. The reason can be understood by looking at Figure 4, where we plot $`c_s^2`$ and $`c_0^2`$ as a function of the radial distance at $`t=50`$ s. For model A the difference between the two quantities is small compared to model B, which develops a quark core and an associated larger composition gradient. As discussed in Section 3, this explains why the frequency of the g1–mode for model B is larger than that of model A.

Finally, it should be mentioned that the frequency of the first w–mode ranges within $`[1.41.8]`$ kHz, which is about one order of magnitude smaller than the typical value of a cold NS, while the damping time is similar, and very short ($`10^4`$ s), much shorter than $`\tau _{diss}`$.

### 4.1 Gravitational wave detectability

The fact that during the early evolution of a PNS the frequencies of the quasi–normal modes are much lower than those of the final NS is of particular interest, since the region where the interferometric antennas of first generation are more sensitive does not extend much beyond 1 kHz. This means that waves emitted at the frequencies of the QNMs of a cold NS will not be detectable (see Eqs. 3).

Unfortunately we do not know how the mechanical energy is distributed among the various modes, but we cannot exclude a priori that at very early times, besides the f–mode, the first g–mode or the first p–mode could give a contribution to the emitted wave, therefore we will consider each mode separately and we will proceed as follows. Having calculated how the frequencies and damping times evolve in time, we can model the gravitational signal emitted by the star oscillating in a given mode as

$$h(t)=h_0\mathrm{e}^{(tt_0)/\tau (t)}\mathrm{sin}[2\pi \nu (t)(tt_0)],$$
(8)
where $`h_0`$ is the initial amplitude at the detector site, $`t_0`$ is the arrival time, $`\nu (t)`$ is the frequency of the oscillation, and $`\tau (t)`$ is the gravitational damping time of the oscillation. The signal–to–noise ratio (SNR) of a given detector is evaluated by using the standard formula for optimal matched filtering

$$SNR=2\left[_0^{\mathrm{}}𝑑\nu \frac{|\stackrel{~}{h}(\nu )|^2}{S_n(\nu )}\right]^{1/2},$$
(9)
where $`S_n(\nu )`$ is the noise power spectral density of the detector. We shall use the noise curves of VIRGO I, LIGO I, GEO600 \[Grishchuk et al. 2001, and references therein\], and that of EURO (http://www.astro.cf.ac.uk/geo/euro/), an interferometric detector with extremely high sensitivity in the kHz region whose performances have been devised by a joint panel of European experts.

We discuss the results of our calculations for model B, but similar results are obtained for model A. Since the emission of gravitational waves is efficient only if the gravitational damping time $`\tau (t)`$ is comparable or shorter than the dissipative timescale $`\tau _{diss}`$, the gravitational signal produced by the various modes has been considered different from zero only if the condition $`\tau (t)\stackrel{<}{}\tau _{diss}`$ is satisfied.

In Table 3 we show the values of the initial amplitude $`h_0`$ that a given signal with the form (8) should have to be detectable with a SNR=5. In addition, we can estimate the amount of energy emitted by a source by integrating over the frequency and the surface the expression of the energy flux

$$\frac{dE_{\mathrm{GW}}}{dSd\nu }=\frac{\pi }{2}\nu ^2|\stackrel{~}{h}(\nu )|^2,$$
(10)
In Table 3 we also give the energy $`E_{GW}`$ that should be emitted by a galactic source oscillating in a given mode to reach SNR=5. Notice that, since the energy scales as $`D^2`$, the energy required to detect the same signal emitted by a source in the Virgo cluster (15 Mpc) with SNR=5 is simply given by $`2.25\times 10^6E_{GW}`$.

From Table 3, we see that a signal detectable by the first generation of ground–based interferometers with SNR=5 must have an initial amplitude larger than $`10^{22}`$. For instance, $`h_0=2.310^{22}`$ for VIRGO I, if emitted by an f–mode oscillation. This amplitude corresponds to an emission of energy of $`3.5\times 10^8M_{}c^2`$, in the event of a galactic Supernova. This is consistent with recent studies of the axisymmetric collapse of the core of a massive star, where the energy emitted in gravitational waves, averaged over all the considered models, is estimated to range within $`[3.610^88.210^8]M_{}c^2`$ \[Dimmelmeier, Font & Müller 2002\]. In the above reference, the corresponding average amplitude for a source at $`10`$ kpc is about $`1.7\times 10^{20}`$, and the average spectrum of gravitational emission is peaked at 930 Hz.

The advanced interferometer EURO would be able to detect signals with initial amplitude $`h_0\stackrel{>}{}810^{25}`$ with SNR=5 (see Table 3). For sources located at the distance of the Virgo cluster this corresponds to an emission energy of the order of $`10^8M_{}c^2`$, with a multiplicative factor depending on the particular fluid mode. Thus, a detector like EURO would increase the detection rate for Supernova events from a few events per century to a few events per month.

The role of the spacetime w1–mode deserves a separate discussion. It is interesting to note that for the w–modes the amount of energy required for the GW to be detected is similar to that of fluid modes (see Table 3) despite the wave amplitude has to be larger. This is essentially due to the smaller gravitational damping time, $`\tau _{w_1}10^4`$ s. However, these modes are very weakly coupled to the fluid, and it remains unclear whether this small amount of energy can be effectively conveyed into spacetime oscillations.

We have also analyzed the response of a resonant mass detector to the gravitational signals computed in this paper. In Figure 5 we plot the experimental noise strain amplitude of EXPLORER obtained in December 2001 (http://www.roma1.infn.it/rog/explorer/), the strain noise curve that the EXPLORER team expects to obtain in 2003 reducing the noise and increasing the quality factor \[Astone 2002\], and the strain amplitude produced by the f– and g1–mode, corresponding to a SNR$`=5`$. From the figure we see that only the f–mode signal falls in the region of best sensitivity. By using the 2003 noise curve, we find that it would be detected with SNR$`=5`$ if the amplitude on Earth is $`h_0310^{20}`$, i.e. $`E_{GW}610^4M_{}c^2`$ for a galactic source. However, this is only an indication; indeed, rotation or different details of the dynamical evolution of the star may shift the two contributions in a way that either the g–modes, or both modes may fall inside the sensitivity band. In addition, better sensitivities will be reached in the next future by the bars further refining the experimental setting. With the expected improvements (M. Bassan, private communication) we find that $`E_{GW}`$ can be lowered by at least an order of magnitude.

## 5  Concluding Remarks

In this paper we show that the frequency and the damping time of the QNMs of a newborn proto–neutron star change significantly as the star cools down and deleptonizes. In particular, at early times, when most of the gravitational energy is emitted, the typical frequencies are considerably lower than those of cold NSs. In order to understand whether the dissipative processes acting in a young PNS could damp the oscillations of the star more efficiently than gravitational waves, we have compared the gravitational damping times to the timescales associated to neutrino viscosity, diffusivity and thermal conductivity ($`\tau _{diss}1020`$ s). We find that, during the first half a second, neutrino processes dominate with respect to the gravitational emission due to the f–mode. However, this time interval is very short compared to $`\tau _{diss}`$, and therefore no much of the energy initially stored into the mode will be dissipated by neutrino diffusion. After half a second, the gravitational damping time of the f–mode becomes shorter than $`\tau _{diss}`$ and, consequently, the remaining mechanical energy can be freely released in gravitational waves. The p–mode damping time is always comparable to $`\tau _{diss}`$, while the gravitational emission from g–modes could be expected only during the first half a second of the evolution, when the g–mode damping time is comparable to $`\tau _{diss}`$. The w–mode damping time is always orders of magnitude smaller than $`\tau _{diss}`$.

We have restricted our analysis to two models of PNS evolution, but we do not expect qualitative changes if different EOSs are used. For example, by changing the symmetry energy of the equation of state one can vary the final composition gradient. For changes in a physically reasonable range, the g1–mode frequency of the final, cold neutron star may vary between 100-300 Hz for the model without a phase transition, and higher values would be reached if a phase transition to quark matter occurs. From our data (Tables 1 and 2) we see that during the first second of evolution $`\nu _{g_1}`$ varies within $`[600800]`$ Hz, $`\nu _f`$ within $`[9001100]`$ Hz, and $`\nu _{g_1}`$ reaches a maximum approximately when $`\nu _f`$ has a minimum; it is reasonable to expect that these features will not change significantly for different initial models.

As shown in Section 4.3, a gravitational signal, emitted at the frequency of the quasi–normal modes of a newborn proto–neutron star in the galaxy, could be detected with a signal to noise ratio of 5 by the first generation interferometers if the energy stored in the modes is greater than $`10^8M_{}c^2`$ (or by a resonant antenna if it is greater than $`10^4M_{}c^2`$). It is interesting to stress that the f–, p– and w–modes have frequencies much lower than those of the cold neutron star which forms at the end of the evolution, in a region where the detectors sensitivity is higher and the energy threshold for detection is lower.

Unfortunately, a galactic Supernova is a very rare event ($`3`$ per century). However, there is another astrophysical scenario to which the results of our study could be applied. Recent numerical simulations \[Shibata & Uryu 2000, Shibata & Uryu 2002\] have shown that short-lived supra-massive neutron stars can also be formed following the merger of two neutron stars with nearly equal masses and low compactness. These highly non axisymmetric objects are supported by differential rotation and would emit quasi-periodic gravitational waves with typical frequencies of $`23`$ kHz. The results of these simulations show that up to 0.01 $`M_{}`$ could be radiated in gravitational waves, as opposed to the 10-8–10-6 $`M_{}`$ usually obtained for stellar core collapse. In these works, the coalescing stars are assumed to have a polytropic EOS and thermal effects are not included, thus the peak frequency corresponds to the excitation of the f-mode frequency of a cold neutron star. The fastly rotating PNS born as a consequence of the coalescence of low compactness neutron stars, should have thermodynamical properties similar to those of the PNSs we study in this paper. Therefore, we think that the effects of thermal and composition gradients we have found should produce a similar shift in frequencies and damping times with respect to the cold star case, but with the remarkable difference that these metastable PNSs would be much stronger emitters of GWs. The large amount of emitted energy would make the ringing down of the PNS observable from the Virgo cluster, increasing the event rate to a few per year.

Acknowledgements

We thank P. Astone for kindly providing the experimental curves of the noise of EXPLORER and for interesting discussions.

This work has been supported by the EU Programme ‘Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137) and the Spanish Ministerio de Ciencia y Tecnologia grant AYA 2001-3490-C02-01. JAP was supported by the Marie Curie Fellowship No. HPMF-CT-2001-01217.


format: text

text:
# X-Ray Evidence for Flare Density Variations and Continual Chromospheric Evaporation in Proxima Centauri

## 1 Introduction

The identification of the physical mechanisms underlying the heating of coronae of magnetically active stars to some 10–100 Million degrees (MK) remains as one of the fundamental problems of high-energy stellar astrophysics. Among the non-steady heating mechanisms, explosive release of energy by flares has received considerable attention in solar physics (e.g., Porter et al. 1995; Krucker & Benz 1998; Parnell & Jupp 2000; Aschwanden et al. 2000), in particular since the statistical ensemble of flares, if extrapolated to small “micro-” or “nano-”flares, may provide enough energy to explain the total coronal energy losses.

There is mounting evidence that flares also play a fundamental role in coronal heating of magnetically active stars. The quiescent X-ray luminosity $`L_\mathrm{X}`$ is correlated with the flare rate in X-rays (Audard et al., 2000) and in the optical (Doyle & Butler, 1985; Skumanich, 1985). Statistically, the hotter X-ray emitting plasma component is more variable than the cooler one, indicating that flares may be involved (Giampapa et al., 1996). Recent studies of the flare frequency in magnetically active stars as a function of the flare energy indicate power-law distributions that may be sufficient to explain all coronal radiative losses (Audard et al., 2000; Kashyap et al., 2002; Güdel et al., 2003).

The physical mechanisms that transport the hot plasma into the corona are, however, less clear. A favored model in solar physics predicts that an initial release of high-energy particles (detected by their radio or prompt optical emission, see below) deposits energy in the cool chromosphere, inducing an overpressure that drives hot material into the corona (“chromospheric evaporation”). As for any coronal heating mechanism, a large density increase is predicted which should be measurable in high-resolution X-ray spectra, in particular in density-sensitive He-like triplets of O VII and Ne IX.

We present novel observations of Proxima Centauri, the nearest star to the Sun (distance = 1.3 pc), with the XMM-Newton satellite, providing new important evidence for frequent chromospheric evaporation. Detailed models will be presented in a forthcoming paper.

## 2 Target and Observations

Proxima Centauri is a dM5.5e dwarf revealing considerable coronal activity. Its “quiescent” X-ray luminosity $`L_\mathrm{X}510\times 10^{26}`$ erg s-1 is similar to the Sun’s despite its $``$50 times smaller surface area. It has attracted the attention of most previous X-ray observatories (Haisch et al., 1980, 1981, 1983, 1990, 1995, 1998; Wargelin & Drake, 2002).

XMM-Newton (Jansen et al., 2001) observed the star on 2001 August 12 during 65 ks of exposure time. For our light curve analysis, we used the most sensitive soft X-ray detector system presently available, namely the three combined European Photon Imaging Cameras (EPIC; Turner et al. 2001; Strüder et al. 2001) with a total effective area of $`2000`$ cm2 at 1 keV and $`2500`$ cm2 at 1.5 keV (with the medium filters inserted). The flares thus recorded are the weakest yet observed on an active star. To minimize CCD pile-up, the small window modes were used. Background light curves were extracted from CCD fields outside the source region. The Reflection Grating Spectrometers (RGS; den Herder et al. 2001) simultaneously recorded high-resolution X-ray spectra between 0.35–2.5 keV, with a resolving power of $`E/\mathrm{\Delta }E300`$ (FWHM) around the O VII lines at $`22`$ Å. The Optical Monitor (OM; Mason et al. 2001) observed the star through its U band filter in high-time resolution mode. All data were analyzed using the XMM-Newton Science Analysis System (SAS versions 5.3 and 5.3.3).

Relevant extracts of the observations are shown in Figures 1 and 2. The first 45 ks of the observation reveal low-level emission around $`L_\mathrm{X}=6\times 10^{26}`$ erg s-1 (0.15-10 keV). Thanks to the high sensitivity, much of the emission is resolved into a succession of weak X-ray flares (Figure 1). The smallest discernible events show $`L_{\mathrm{X},0.1510}=2\times 10^{26}`$ erg s-1 and an integrated X-ray energy loss of $`E_{\mathrm{X},0.1510}1.5\times 10^{28}`$ erg, corresponding to modest solar flares. Most importantly, many (although not all) of these flares are preceded by a short pulse in the optical U band (Figure 1). A large flare monitored almost in its entirety governed the final 20 ks of the observation, with a peak $`L_{\mathrm{X},0.1510}3.7\times 10^{28}`$ erg s-1 and $`E_{\mathrm{X},0.1510}1.5\times 10^{32}`$ erg (Figure 2). Due to pile-up in the MOS detectors, we used only the PN data for the light curve analysis of this flare. Again, a large optical burst accompanies the X-ray flare during its increase. The figure also shows spectral extracts of the He-like triplets of O VII and Ne IX referring to different phases of the large flare. The line flux ratios are clearly variable, an aspect that will be discussed in detail below.

## 3 Results and Discussion

We now seek an interpretation within the evaporation scenario. During a flare, high-energy electrons are thought to be accelerated in the corona from where they travel toward the magnetic footpoints in the chromosphere. They induce a prompt signal in the optical U band (3000-4000 Å), probably owing to Balmer emission or blackbody radiation, or both (Hudson et al., 1992; Neidig & Kane, 1993; Hawley et al., 1995). Their bulk kinetic energy heats the chromospheric gas to several MK, building up a pressure gradient that drives the hot plasma into the corona, thus increasing the coronal density and the soft X-ray emission. Since the optical signal $`L_\mathrm{O}`$ approximately traces the rate of energy deposited by the electrons (Hudson et al., 1992), and the X-ray losses $`L_\mathrm{X}`$ are a measure of the accumulated coronal thermal energy, the time integral of the optical signal should resemble the X-ray light curve,

$$L_\mathrm{X}(t)_0^tL_\mathrm{O}(t^{})𝑑t^{}$$
(1)
where the integration starts at a point in time preceding the flare ($`t^{}=0`$). In other words, the time derivative of the X-ray light curve should mimic the optical signal:

$$\frac{dL_\mathrm{X}}{dt}L_\mathrm{O}.$$
(2)
Equation 1 is appropriate for testing faint signals, whereas equation 2 requires strong signals but is sensitive to details in the light curves. This important diagnostic has been observed in the solar corona, known as the Neupert Effect (using radio or hard X-ray emissions instead of U band radiation; Neupert 1968; Dennis & Zarro 1993), but has rarely been detected in stellar observations, especially in low-level emissions (Doyle et al., 1986, 1988), presumably owing to the limited sensitivity of previous X-ray devices. Qualitative evidence is seen in some strong flares (e.g., de Jager et al. 1986, 1989; Schmitt et al. 1993; Hawley et al. 1995). We suggest that each of our optical signals is the equivalent of an impulsive burst, while the associated X-ray flares are the signatures of the subsequent chromospheric evaporation. A surprising aspect is the high repetition rate of such processes in Proxima Centauri.

While hot plasma is accumulated in the corona, the X-ray emission of any plasma packet decays nearly exponentially due to cooling. Equation 1 in this form does not consider cooling; we therefore include a decay term by convolving the optical light curve with a cut-off exponential decay profile ($`e^{t/\tau }`$ for $`t0`$ and zero for $`t<0`$). A decay time of $`\tau =200`$ s is found to be appropriate, and the convolution is shown by the thick solid curves in Figures 1a-f. After the convolution, a high level of correlation with the X-rays reveals itself in many coincident peaks and similar light curve shapes (Figures 1a,b). Suggestive examples are marked by vertical dotted lines. The optical light curve often predicts the subsequent X-ray behavior, revealing an intimate connection between the two emissions. The correlation holds even for some flare groups (Figures 1c,d). In some cases, the detailed correlation is poor although optical flare groups cluster during X-ray flux increases (Figures 1e,f). The large flare allowed us to compute its time derivative (Figure 2b). It traces the optical signal closely until after the optical flare peak, when uncorrelated cooling mechanisms take over. Here, $`\tau `$ was not considered as it appears that ongoing heating on time scales longer than $`\tau `$ dominates over cooling.

The X-ray flux $`F`$ from a collisional plasma is proportional to the volume emission measure (EM), $`F=\mathrm{const}n_en_\mathrm{H}𝑑V0.84\mathrm{const}n_e^2𝑑V`$ where $`n_e`$ is the electron density, $`n_\mathrm{H}`$ is the hydrogen (or proton) density, and $`dV`$ is the source volume element; the proportionality constant is the line emissivity (Mewe et al., 1985). The determination of $`V`$ and $`n_e`$ from fluxes alone is therefore a degenerate problem. It can be circumvented by using the He-like spectral line triplets as a diagnostic for $`n_e`$. The relevant triplets of the O VII and Ne IX ions are each formed by radiative decays to the ground state, namely the resonance transition $`r`$ ($`1s^2{}_{}{}^{1}S_{0}^{}1s2p{}_{}{}^{1}P_{1}^{}`$), the intercombination transition $`i`$ ($`1s^2{}_{}{}^{1}S_{0}^{}1s2p{}_{}{}^{3}P_{1,2}^{}`$), and the forbidden transition $`f`$ ($`1s^2{}_{}{}^{1}S_{0}^{}1s2s{}_{}{}^{3}S_{1}^{}`$). In cool stars with negligible ultraviolet continuum radiation, the flux ratio $`R=f/i`$ is sensitive to $`n_e`$ because electron collisions can excite the $`{}_{}{}^{3}S_{1}^{}`$ to the $`{}_{}{}^{3}P`$ state before the former decays radiatively (Gabriel & Jordan, 1969; Porquet et al., 2001). The data reveal strong variations of $`R`$ along the flare (examples in Figure 2c). The characteristic $`n_e`$ derived from O VII reaches $`4\times 10^{11}`$ cm-3 during the flux maxima, but is much lower ($`2\times 10^{10}`$ cm-3) during the decays (crosses in Figure 2a). Previous attempts to estimate stellar coronal flare densities spectroscopically used lines in the extreme ultraviolet range, but density variations were difficult to constrain (e.g., Monsignori-Fossi et al. 1996; for a marginal X-ray result, see Stelzer et al. 2002). Qualitatively similar variations are seen in the Ne IX triplet (Figure 2c), although blends of the $`i`$ line thwart a precise quantitative interpretation. The strongly variable $`f`$ line alone suggests that plasma around 4–5 MK (the maximum formation temperature of Ne IX) attained densities across the density sensitivity range of this triplet, $`10^{11}2\times 10^{12}`$ cm-3 (Porquet et al., 2001).

To derive the associated EM, we measured the total flux $`F_{\mathrm{OVII}}`$ of the O VII lines used above. We found that $`F_{\mathrm{OVII}}`$ is nearly proportional to the 0.15–10 keV EPIC count rate. This is important because O VII lines are efficiently formed only in a temperature range of $`1T4`$ MK. We conclude that variations of the EM-weighted formation temperature of these lines are not the dominant factor for the variability of $`F_{\mathrm{OVII}}`$. The ratio $`G=(f+i)/r`$ provides an independent test since it is a temperature indicator for the O VII line-forming plasma (Gabriel & Jordan, 1969; Porquet et al., 2001). We infer characteristic temperatures, $`T_G=1.53.5`$ MK, that cluster around the maximum formation temperature ($`T_{\mathrm{OVII}}2`$ MK) but do not show a clear correlation with $`F_{\mathrm{OVII}}`$. We will therefore approximate the O VII emitting plasma as being isothermal at $`T_{\mathrm{OVII}}`$; systematic errors thus introduced for the EM are no larger than a factor of $``$2 (see tables in Mewe et al. 1985). For a homogeneous source at $`T_{\mathrm{OVII}}`$, $`F_{\mathrm{OVII}}n_e^2V`$. Therefore, $`VF_{\mathrm{OVII}}/n_e^2`$, and the total mass $`M`$ involved is $`MF_{\mathrm{OVII}}/n_e`$. We used emissivities from Mewe et al. (1985) referring to solar photospheric composition. From spectral fits in XSPEC (Arnaud, 1996), we found an abundance of oxygen of $`0.52\pm 0.08`$ times the solar photospheric value (Anders & Grevesse, 1989). For the primary peak, its decay, the secondary peak, and its decay, the masses of O VII emitting plasma thus amount to $`M4.8\times 10^{14}`$ g, $`2.1\times 10^{15}`$ g, $`3.2\times 10^{14}`$ g, and $`2.5\times 10^{15}`$ g, respectively, while $`V7.0\times 10^{26}`$ cm3, $`2.3\times 10^{28}`$ cm3, $`4.3\times 10^{26}`$ cm3, and $`5.3\times 10^{28}`$ cm3, respectively. The potential energy of a $`10^{15}`$ g plasma packet at a height of $`9\times 10^3`$ km (the pressure scale height at $`2\times 10^6`$ K for a radius of $`10^{10}`$ cm and a stellar mass of 0.2$`M_{}`$, see Haisch et al. 1983 and references therein) is $`2.5\times 10^{29}`$ erg. This compares with the total thermal energy content of the same plasma, $`E1.65\times (3/2)MkT/m_\mathrm{H}4\times 10^{29}`$ erg, while the total radiated X-ray energy from the cooler plasma component at 0.2 keV (determined from a spectral fit) is approximately $`1.7\times 10^{31}`$ erg during the complete flare (using a peak EM of $`1\times 10^{50}`$ cm-3, a decay time of 4300 s, and a cooling loss rate of $`4\times 10^{23}`$ erg s-1cm3), necessitating considerable replenishment of the cool material during the flare.

Even after consideration of the uncertainties due to variations of $`T_G`$, it is clear that densities, masses and volumes of the plasma accessible by O VII increase during the flare. Coronal density increases are a consequence of any heating mechanism in closed loops. It is the preceding optical impulsive emission that strongly supports chromospheric evaporation induced by electron beams. Solar white-light flares are closely correlated in time with hard X-ray bursts, suggesting that this phase relates to electron bombardment of the chromosphere (Hudson et al., 1992; Neidig & Kane, 1993), while soft X-rays and higher densities are evolving more gradually. Since plasma cooling becomes important early in the flare development, the relatively cool O VII emitting plasma component is not only augmented due to heating of chromospheric material, but also due to cooling of hotter material that has initially reached temperatures beyond the formation range of O VII lines, hence the continuous increase of masses and volumes during the decay.

The low flux of the weaker flares and their strong mutual overlap prohibit individual density measurements, but integration of larger segments shows appreciable changes in $`R`$ as well, with the highest-density episode related to the segment in Figure 1b (an increase of at most a factor of ten over the average). Furthermore, the time behavior of several optical and X-ray flares that expresses itself as a Neupert Effect is similar to the behavior seen in the large flare, indicating that the same processes repeat at a much higher cadence at lower energy levels such that they govern the X-ray emission most of the time. We do not find a one-to-one correspondence between all optical and X-ray flares, but this is no different on the Sun (Dennis & Zarro, 1993). Some groups of small optical flares occur, however, during a general increase of the soft X-ray light curve.

The present observations suggest an important role of chromospheric evaporation in active stellar coronae. The resulting high plasma densities appear to be a critical element for the high X-ray luminosities of magnetically active stars.

Helpful comments by the referee are acknowledged. Research at PSI has been supported by the Swiss National Science Foundation (grant 2000-058827). The present project is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).


text:
# Sand stirred by chaotic advection

## Abstract

We study the spatial structure of a granular material, $`N`$ particles subject to inelastic mutual collisions, when it is stirred by a bidimensional smooth chaotic flow. A simple dynamical model is introduced where four different time scales are explicitly considered: i) the Stokes time, accounting for the inertia of the particles, ii) the mean collision time among the grains, iii) the typical time scale of the flow, and iv) the inverse of the Lyapunov exponent of the chaotic flow, which gives a typical time for the separation of two initially close parcels of fluid. Depending on the relative values of these different times a complex scenario appears for the long-time steady spatial distribution of particles, where clusters of particles may or not appear.

Many different physical processes can be studied in the general framework of transport of finite-size particles by an external flow. Sand suspended on the surface of a river or driven by the wind, chemical pollutants advected by atmospheric flows, and impurities driven by oceanic currents, are just but a few of the many important real examples that fit into this.

Recently, much effort coming from the Chaotic Advection gen:adv  and Turbulence gen:tur  communities has been devoted to the study of advection of non-reacting inertial particles by chaotic or turbulent flows. Though some works considered the effect of elastic collisions, see for example collins , the obvious fact that in real systems collisions among particles dissipate kinetic energy into heat (inelastic collisions) has not been, to our knowledge, much subject of study. From another perspective, the so-called Granular Media community has also widely considered granular  the properties of particles colliding inelastically, and externally driven, for example, by a shear flow or by vibrating periodically the container of the particles campbell . However, the influence of an external chaotic or turbulent flow seems to be overcome in the literature (see, however, shinbrot1  for the case of densely packed granular systems).

In this paper, we try to investigate numerically the full problem. Equivalently, we study the effect of collisions on inertial particles immersed in a flow showing chaotic advection aref , or rather, the properties of a dilute granular system (in the literature the term granular gas gas  is widely used) under chaotic stirring. It is by now clear that transport behavior in the so-called Batcherlor’s regime, that is, in the range of scales between the smallest typical length scale of the velocity field and the characteristic diffusion length scale, is equivalent in a turbulent flow and a chaotic advection flow.

We introduce a simple model of granular material consisting of many particles subject to inelastic mutual collisions, and immersed in a smooth two-dimensional chaotic flow without gravity. Granular particles are assumed to be much heavier than fluid. In particular, we focus on the steady spatial structure of the system. It is known, in fact, that the presence of inertia can induce preferential concentration gen:adv ; gen:tur . One may ask what happens if the diffusive role of collisions (in the elastic limit) among particles is taken into account. At the same time is also well known that dissipative collisions originate particle clustering puglio . Therefore, in spite of its simplicity, we will see that the model shows many transitions from clustering to unclustering depending on the interplay between the inertia of the particles, the chaotic flow and the collisions.

We consider $`N`$ identical particles of mass $`m=1`$ on a two-dimensional domain $`L\times L`$ (with periodic boundary conditions) driven by an external velocity field $`𝐮(𝐱,t)`$

$`{\displaystyle \frac{d𝐯_i(t)}{dt}}`$ $`={\displaystyle \frac{1}{\tau }}\left(𝐯_i(t)𝐮(𝐱_i(t),t)\right),`$ (1a)
$`{\displaystyle \frac{d𝐱_i(t)}{dt}}`$ $`=𝐯_i(t),`$ (1b)
where $`i=1,\mathrm{},N`$, $`𝐯_i`$ is the velocity of the particle $`i`$, $`\tau `$ is the standard Stokes time, that depends on the particles diameter, viscosity, and densities of fluid and particles. Hence the term in the r.h.s of Eq. (1a) is simply a viscous Stokes drag. In addition the particles are subject to binary instantaneous inelastic collisions according to the rule

$$𝐯_{i(j)}^{}=𝐯_{i(j)}\frac{1+r}{2}((𝐯_{i(j)}𝐯_{j(i)})\widehat{𝐧})\widehat{𝐧}$$
(2)
where $`i`$ and $`j`$ are the indexes of colliding particles, the primes denote the velocities after collision and $`r[0,1]`$ is the restitution coefficient ($`r=1`$ for the elastic case). This rule ensures the conservation of momentum, while the kinetic energy is dissipated if $`r<1`$.

Some further clarifications on the model are worth mentioning. In the absence of collisions the equations in (1) are simply the equations of motion of a rigid sphere in a flow where the Faxen corrections, the added mass term, and the Bernoulli term are neglected maxey . This is a consistent approximation when the particles of granular material are much heavier than that of the fluid, which is the case we are considering. Moreover, we adopt the Direct Simulation Monte Carlo (DSMC) scheme to integrate the dynamics of the system bird . The DSMC algorithm, widely used in numerical hydrodynamics, consists of a time discretization with fixed time step: in every time step there is a free flow-step and a collision step. In the free-flow step the motion of particles is integrated according to equations (1) without taking into account the possible collisions. In the collision step every particle has a probability of colliding with neighboring particles (in a disk with a radius smaller than the mean free path): the probability of particle $`i`$ to collide with particle $`j`$ is proportional to the relative velocity $`|𝐯_i𝐯_j|`$. This scheme has been proved to converge to the solution of the Boltzmann equation of the corresponding hard disk gas and its key feature is (through the randomization of collisions) the assumption of Molecular Chaos, i.e. lack of correlations for colliding particles, $`P(𝐯_i,𝐯_j,t)=P(𝐯_i,t)P(𝐯_j,t)`$. It is important to note that this assumption does not rule out correlations at scales larger than the scale of a diameter, and the DSMC is often used to study vortices vortices  and other instabilities. In the context of inelastic gases it has been previously used also for the investigation of the clustering phenomenon puglio .

A statistically stationary state is reached when the dissipation of energy because of the inelastic collisions is balanced with the continuous energy injection coming from the external flow note . We determine that the statistically stationary state is reached when the typical fluctuations of the total energy of the system, $`E=v_i^2`$ (the so-called granular temperature), are rather small, typically some small percentage of $`E`$. As already mentioned, the properties of this state are the central objective of this work, with special emphasis to the clusterization of particles.

The formal introduction of the collision operator and adimensionalization of the dynamical system Eqs. (1a,1b) will help in order to see the relevance of the distinct terms. If we denote $`T`$ as the typical time of the flow, $`L`$ the typical length and $`𝐅_{\mathrm{𝐜𝐨𝐥}}/\tau _𝐜`$ the collision operator, with $`\tau _c`$ the mean collision time, we can rewrite Eq. (1a) as follows

$$\frac{\tau }{T}\frac{d\widehat{𝐯}_i}{d\widehat{t}}=(\widehat{𝐯}_i\widehat{𝐮})+\frac{\tau }{\tau _c}𝐅_{\mathrm{𝐜𝐨𝐥}},$$
(3)
where $`\widehat{t}=t/T`$, $`\widehat{𝐯}_i=𝐯_iT/L`$, and $`\widehat{𝐮}=𝐮T/L`$.

For our numerical studies we use the time periodic flow $`𝐮(𝐱,t)=(u_x(x,y,t),u_y(x,y,t))`$ given by:

$`u_x`$ $`=`$ $`2\pi L\alpha /T\mathrm{sin}(2\pi y/L),`$
$`u_y`$ $`=`$ $`0,`$ (4)
if $`t`$ mod $`T<T/2`$, and

$`u_x`$ $`=`$ $`0,`$
$`u_y`$ $`=`$ $`2\pi L\alpha /T\mathrm{sin}(2\pi x/L),`$ (5)
if $`t`$ mod $`T>T/2`$.

Periodic boundary conditions are assumed, and for $`\alpha >0.4`$ a unique chaotic region without KAM tori (numerically distinguishable) is obtained berthier . The Lyapunov exponent for the flow is given by $`\lambda =1.96\mathrm{ln}(3.35\alpha )/T`$ and introduces a new time scale, $`\tau _f=1/\lambda `$, for the exponential separation rate of close fluid parcels. Our results has been checked to be independent of the kind of flow. In particular, we have alternatively used the cellular flow derived from the stream function $`\mathrm{\Psi }(x,y,t)=U\mathrm{sin}(\frac{n\pi x+B\mathrm{cos}(\omega t)}{L})\mathrm{sin}(\frac{n\pi y}{L})`$, where $`U`$ is the maximal velocity of the flow, $`n`$ the number of cells, $`B`$ the strength of the time dependent forcing and $`\omega `$ its frequency gollub .

Depending on the relative values of the above introduced time scales: $`\tau `$, $`T`$, $`\tau _f`$ and $`\tau _c`$, a very complex scenario appears for the spatial distribution of the particles in the stationary state. Some preliminary general results can be anticipated. For example, it seems obvious that when $`\tau _f<\tau _c`$ the exponential separation produced by the chaotic driving of the fluid may avoid any clusterization. The opposite may happen if $`\tau _c<\tau _f`$ as the mechanism of aggregation due to inelastic collisions resists the dispersing effect of the chaotic flow. However, as we will see it is also very important to take into account the compaction due to the inertia. In any case, discussing all the different possibilities (taking also into account some values of the parameter $`r`$) goes beyond the scope of this Paper. Thus, we just restrict ourselves to study some of the most interesting cases.

First we consider the case of irrelevant inertia, $`\tau <<T`$. Obviously, in the absence of collisions, particles follow strictly the fluid and no aggregation appears (see Fig. 1a)). In fact, $`\tau `$ is the relaxation time of the particles to the driving flow; without collisions the particles are continuously flowing through the entire spatial domain since no KAM tori are present in the flow. However, everything changes when the collisions are taken into account. Thus, if $`\tau _c<\tau `$, the last term in Eq. (3) is the most relevant one. If we have also that $`\tau _c<\tau _f`$ an aggregation of the particles is always observed for any value of $`r`$ (even, and most noticeably, at $`r=1`$). This can be understood from Eq. (3), since now the velocity of the particles can be written as:

$$\widehat{𝐯}\widehat{𝐮}+\frac{\tau }{\tau _c}𝐅_{\mathrm{𝐜𝐨𝐥}}.$$
(6)
Thus, the effect of collisions is similar to that of inertia when considered alone: the velocity of the particles is modified respect to that of the fluid parcels, then it is no longer compressible ($`\widehat{𝐯}\mathrm{𝟎}`$), and there are regions where particles aggregate. Figs. 1b) and 1c) show a pattern of the distribution of particles for $`r=0.99`$ and $`r=0.6`$, respectively. It is also very important to note here a specific behavior of elastic or quasi-elastic collisions ($`r1`$). The effect of multiple elastic collisions is equivalent to macroscopic diffusion, which, if large enough, may induce a dispersing mechanism that prevents clustering. Summing up, one can say that elastic collisions induce a clustering inertia-like effect, that disappears when diffusion is strong enough. Quantitatively this is controlled by the adimensional parameter $`\tau /\tau _c`$ since the macroscopic diffusivity $`D1/\tau _c`$. Increasing $`\tau /\tau _c`$ clustering like that of Fig. 1b) disappears and a homogeneous distribution is attained (see discussion below on the $`P(m)`$ function and Fig. 2a)).

On the contrary, when $`\tau _f<\tau _c`$ (and also $`\tau _c<<\tau `$) the strong chaoticity of the flow dominates and avoids the formation of any clustering in the system, see Fig. 1d). Finally, it is obvious from Eq.(3) that the limit $`\tau <<\tau _c`$ is equivalent to the absence of collisions, thus no clusterization appears.

A more quantitative analysis of clustering follows. This is performed via a particle-in-cell histogram, i.e., after dividing the system in $`M`$ small boxes (we use $`N/M=5`$) an histogram is made with the number of boxes containing $`m`$ particles, denoted the corresponding function as $`P(m)`$. For an homogeneous system of particles $`P(m)`$ is a Poisson distribution $`exp(\lambda )\lambda ^m/m!`$, where $`\lambda =N/M`$. As the clustering is stronger developed, the deviations from a Poissonian are more evident. Fig. 2a) calculates $`P(m)`$ for the patterns shown in Fig. 1, and one more case. We observe the appearance of clustering for $`r<1`$, and even for $`r1`$. Also, it is shown how for the elastic case the Poisson distribution is approached when $`\tau /\tau _c<<1`$.

It is very relevant to know the effect of collisions in the context of transport of inertial particles in chaotic flows. Therefore, we next consider the general case where the inertial term is relevant, $`\tau >T`$ (but otherwise $`\tau `$ small enough such that the l.h.s is not negligible with respect to first term in the r.h.s). It is well-known gen:tur  that, in the absence of collisions, heavy particles accumulate in regions of low strain and high vorticity, the so-called preferential concentration phenomenon. In our case, the clustering steady pattern is shown in Fig. 3a). After including collisions one has: if $`\tau <<\tau _c`$ no effect is observed as the last term in the r.h.s of Eq. (3) can be neglected. A different behavior is observed when $`\tau _c<<\tau `$. Regardless of the value of $`\tau _f`$, in the limit of elastic collisions ($`r1`$) no aggregation of the particles is observed, Fig. 3b). As explained before, when $`\tau /\tau _c`$ is large, the continuous elastic colliding processes of the particles induce a diffusive effect that prevents clustering. On the contrary, as we go into the inelastic regime, $`r<1`$, the chaotic time scale, $`\tau _f`$, now plays a very important role. Thus, when $`\tau _c<\tau _f`$ particles aggregate although in a very different pattern to that obtained in the absence of collisions (Fig. 3c)). In particular, at variance with clustering in the absence of collisions (see Fig. 3a)), clusters are created and destroyed continuously, and moving with the flow. In this case, the particles aggregate thanks to the inelastic collisions and because the flow is not chaotic enough to disperse completely the clusters. The opposite is observed, Fig. 3d), when $`\tau _f<\tau _c`$; the clusters do not resist the dispersion of the flow and the density of particles is homogeneous. The corresponding $`P(m)`$ function is plotted in Fig. 2b).

For completeness, we briefly mention the limit $`\tau >>T`$. This is completely equivalent to a freely cooling granular gas since the external flow is irrelevant. Thus, clustering is usually observed for inelastic collisions goldhirsch .

In this letter we have studied numerically the stationary spatial structure of a granular material advected by a chaotic flow. We have seen that the relative importance of collisions, inertia and chaoticity of the flow gives rise to transitions from homogeneous to inhomogeneous density of grains. We think that our simple model is relevant for real physical phenomena that involve the transport of a large amount of heavy particles, like for example those appearing in many geophysical situations.

###### Acknowledgements.

We acknowledge fruitful discussions with Andrea Baldassarri, Fabio Cecconi, Umberto Marini Bettolo Marconi, and Angelo Vulpiani, and Emilio Hernández-García for a critical reading of the manuscript. C.L. acknowledges support from the Spanish MECD. A.P. acknowledges support from the INFM Center for Statistical Mechanics and Complexity (SMC).


text:
# Real space statistical properties of standard cosmological models

## 1 Introduction

In standard theories of structure formation in cosmology the density field in the early Universe is described as a perfectly uniform and isotropic matter distribution, with superimposed tiny fluctuations characterized by some particular correlation properties (e.g. pee93 ). These fluctuations are believed to be the initial seeds from which, through a complex dynamical evolution, galaxies and galaxy structures have emerged. In particular, in all the standard models, the initial fluctuations are taken to be Gaussian and with a power spectrum (PS) $`P(k)`$ satisfying the so-called Harrison-Zeldovich (HZ) condition of being proportional to $`k`$ at small $`k`$ har .

The present paper has two main porpuses prb . Firstly, to clarify the real space statistical properties of the mass density fluctuations common to all the standard cosmological models (which have been almost completely overlooked in the literature on the subject). And secondly, through this discussion, to relate and compare these models of the primordial Universe to correlated systems encountered in statistical physics.

In particular we find that all these standard cosmological models are characterized by a “superhomogeneous” (or superuniform) matter distribution. This means that mass fluctuations over sufficiently large spatial scales are “sub-poissonian”, i.e. they increase with the spatial scale more slowly than in a random poissonian matter distribution. Moreover, it is shown that the scaling of mass fluctuations satisfying the HZ condition approaches the slowest possible for any stochastic matter distribution. We will see that in the context of usual statistical physics this kind of behavior is recovered in the case of lattice-like or glass-like particle distributions or in the so called One Component Plasma (OCP) ocp  at thermal equilibrium. In this context we present some simple recipe to build particle distributions satisfying the HZ condition at small $`k`$. This last point can be useful first of all in the context of $`N`$-body simulations for the study of structure formation from primordial fluctuations through a gravitational dynamics thierry .

## 2 Basic concepts of correlation analysis

Inhomogeneities of the mass density field in cosmology are described using the general framework of stationary stochastic processes (hereafter SSP). Let us consider in general the description of a continuous or a discrete homogeneous mass distribution $`\rho (\stackrel{}{r})`$ in terms of such a process. A stochastic process is completely characterized by its “probability density functional” $`𝒫[\rho (\stackrel{}{r})]`$ which gives the probability that the result of the stochastic process is the density field $`\rho (\stackrel{}{r})`$ (e.g. see Gaussian functional distributions Vanmarcke ). For a discrete mass distribution, the space (e.g. the infinite three dimensional space) is divided into sufficiently small cells and the stochastic process consists in occupying or not any cell with a point-particle following certain correlated probabilities, and $`\rho (\stackrel{}{r})`$ can be written in general as:

$$\rho (\stackrel{}{r})=\underset{i=1}{\overset{\mathrm{}}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i),$$

where $`\stackrel{}{r}_i`$ is the position vector of the particle $`i`$ of the distribution .

The word “stationarity” refers in the present context to the spatial stationarity of the process, and means that the functional $`𝒫[\rho (\stackrel{}{r})]`$ is invariant under spatial translation. We suppose also that the distribution is statistically isotropic (invariance of $`𝒫[\rho (\stackrel{}{r})]`$ under spatial rotation), and has a well defined positive average value:

$$\rho (\stackrel{}{r})=\rho _0>0,$$

where $`\mathrm{}`$ is the ensemble average over all the possible realizations of the stochastic process, i.e. the average over the functional $`𝒫[\rho (\stackrel{}{r})]`$. Moreover it is usually assumed that $`𝒫[\rho (\stackrel{}{r})]`$ is ergodic. This means that spatial averages in single infinite realization coincide with the ensemble averages.

The quantity $`\rho (\stackrel{}{r_1})\rho (\stackrel{}{r_2})\mathrm{}\rho (\stackrel{}{r_l})`$ is called the complete $`l`$-point correlation function. In the discrete case $`\rho (\stackrel{}{r_1})\rho (\stackrel{}{r_2})\mathrm{}\rho (\stackrel{}{r_l})dV_1,dV_2,\mathrm{},dV_l`$ gives the a priori probability of finding $`l`$ particles, in a single realization, placed in the infinitesimal volumes $`dV_1,dV_2,\mathrm{},dV_l`$ respectively around $`\stackrel{}{r_1},\stackrel{}{r_2},\mathrm{},\stackrel{}{r_l}`$. In both cases the statistical stationarity and isotropy (hereafter SSI) imply that the $`l`$-point correlation functions, for any $`l`$, depend only on the scalar relative distances among the $`l`$ points gsl01 .

Let us analyze in further detail the auto-correlation properties of these systems. As aforementioned, in our hypothesis, $`\rho (\stackrel{}{r_1})\rho (\stackrel{}{r_2})`$ depends only on $`r_{12}=|\stackrel{}{r_1}\stackrel{}{r_2}|`$.

The reduced two-point correlation function $`\stackrel{~}{\xi }(r)`$ is defined by:

$$\rho (\stackrel{}{r_1})\rho (\stackrel{}{r_2})\rho _0^2\left[1+\stackrel{~}{\xi }(r_{12})\right].$$

The correlation function $`\stackrel{~}{\xi }(r)`$ is one way to measure the ”persitence of memory” of spatial variations in the mass density huang . For a discrete distribution of particles this means that $`\rho _0\stackrel{~}{\xi }(r)d^3r`$ measures the excess of probability of finding a particle in a small volume $`d^3r`$ at a distance $`r`$ from another fixed particle, with respect to the a priori probability $`\rho _0d^3r`$. Obviously, since $`\rho _0`$ must be the average density in any single realization in the infinite volume limit, $`\stackrel{~}{\xi }(r)`$ satisfies the condition of vanishing in the limit $`r\mathrm{}`$.

### 2.1 The Poisson particle distribution

Let us consider the paradigm of a stochastic homogeneous point-mass distribution: the Poisson case (see gsl01 ). It can be defined as follows: let us partition the space in cubic cells of volume $`\mathrm{\Delta }V`$, and then occupy each cell independently of the others, with a point-particle with a probability $`n_0\mathrm{\Delta }V`$, where $`n_0>0`$ and $`\mathrm{\Delta }Vn_0^1`$. It is simple to show that in the limit $`\mathrm{\Delta }V0`$, one can write:

$`\rho _0=n_0`$ (1)
$`\stackrel{~}{\xi }(r)={\displaystyle \frac{\delta (\stackrel{}{r})}{\rho _0}}`$ (2)
where $`\delta (\stackrel{}{r})`$ is the usual Dirac delta function. Equation (2) is a direct consequence of the fact that there is no correlation between different spatial points in the definition of the stochastic process. That is, the reduced correlation function $`\stackrel{~}{\xi }`$ has only the so called diagonal part. It is simply shown that this diagonal part is present in the reduced two-point correlation functions of any statistically homogeneous discrete distribution of particles with correlations. Therefore in general Landau  for a SSI distribution of particles the reduced correlation function can be written as

$$\stackrel{~}{\xi }(r)=\frac{\delta (\stackrel{}{r})}{\rho _0}+\xi (r)$$

where $`\xi `$ is the non-diagonal parts which is meaningful only for $`r>0`$.

## 3 The mass variance in a sphere

In this section we consider the amplitude of the mass fluctuations in a generic sphere of radius $`R`$ with respect to the average mass. First let $`M(R)=_{C(R)}\rho (\stackrel{}{r})d^3r`$ be the mass (for a discrete distribution the number of particles) inside the sphere $`C(R)`$ of radius $`R`$ and then volume $`C(R)=\frac{4\pi }{3}R^3`$. The normalised mass variance is defined as

$$\sigma _M^2(R)=\frac{M(R)^2M(R)^2}{M(R)^2}.$$

It is simple to show that $`\sigma _M^2(R)`$ in a stationary stochastic mass density can be rewritten as:

$$\sigma _M^2(R)=\frac{1}{C(R)^2}_{C(R)}d^3r_1_{C(R)}d^3r_2\stackrel{~}{\xi }(|\stackrel{}{r_1}\stackrel{}{r_2}|).$$
(3)
This formula will be useful in the following for a complete classification of the mass density fields with respect to the large scale behavior of mass fluctuations.

Note that since $`\stackrel{~}{\xi }(r)0`$ for $`r\mathrm{}`$, then

$$\underset{R\mathrm{}}{lim}\sigma _M^2(R)=0.$$

This is nothing but the condition of existence of a well defined average density $`\rho _0`$.

### 3.1 A classification of mass fluctuations

If we apply Eq. (3) to the case of a Poisson point process as above defined, we find

$$\sigma _M^2(R)=\frac{1}{\rho _0C(R)}\frac{1}{M(R)},$$

that is $`\mathrm{\Delta }M^2(R)C(R)`$. The same widespread behavior is found in systems characterized by mainly positive short range correlations (e.g. a perfect gas at high temperature).

We can characterize an arbitrary mass distribution with respect to the behavior of mass fluctuations as follows. Let us suppose $`\sigma _M^2(R)R^\alpha `$ at large $`R`$, then:

* If $`0<\alpha <3`$ the system has critical fluctuations, typical of the order parameter of a thermodynamical system at the critical point of a second order phase transition. As we show below, this happens in the case of mainly positive long range correlations such that $`d^3r\stackrel{~}{\xi }(r)=+\mathrm{}`$;
* If $`\alpha =3`$ we have a substantially poissonian system, characterized by mainly positive short range correlations, such that $`d^3r\stackrel{~}{\xi }(r)=c`$ where $`0<c<+\mathrm{}`$;
* If $`\alpha >3`$ the system has sub-poissonian fluctuations typical of lattice or glass-like systems with an almost ordered arrangement of mass perturbations (of particles in the discrete case), characterized by a balance between positive and negative correlations, such that $`d^3r\stackrel{~}{\xi }(r)=0`$. For this reason we call these systems superhomogeneous. We show also below that $`\alpha <4`$ ($`\alpha <d+1`$ in $`d`$ dimensions) for any stochastic mass distribution. As shown below all the standard cosmological models for primordial mass fluctuations belong to this class with $`\alpha =4`$ (for a difference between a Poisson and a superhomogeneous distribution see Fig. 1).

## 4 The power spectrum $`P(k)`$

Much more used in cosmology than $`\stackrel{~}{\xi }(r)`$ is the equivalent $`k`$-space quantity, the power spectrum $`P(\stackrel{}{k})`$ which is defined as

$$P(\stackrel{}{k})=\underset{V\mathrm{}}{lim}\frac{|\delta _\rho (\stackrel{}{k})|^2}{V}$$

where $`\delta _\rho (\stackrel{}{k})=_Vd^3re^{i\stackrel{}{k}\stackrel{}{r}}(\rho (\stackrel{}{r})\rho _0)/\rho _0`$ is the Fourier integral in the volume $`V`$ of the normalized fluctuation field $`(\rho (\stackrel{}{r})\rho _0)/\rho _0`$. In a SSI matter distribution this depends only on $`k=|\stackrel{}{k}|`$ and can be written as

$$P(\stackrel{}{k})P(k)=d^3re^{i\stackrel{}{k}\stackrel{}{r}}\stackrel{~}{\xi }(r)=4\pi _0^{\mathrm{}}𝑑rr^2\frac{\mathrm{sin}(kr)}{kr}\stackrel{~}{\xi }(r).$$
(4)
It follows from its definition that $`P(k)0`$ (this property and Eq. (4) constitute the so-called Kh$`\widehat{\text{ı}}`$ncin theorem of SSP’s). In particular for a Poisson particle distribution one has $`P(k)=\rho _0^1`$. At last note the important relation $`P(0)=d^3r\stackrel{~}{\xi }(r)`$ for any SSP.

### 4.1 The cosmological power spectrum, and the HZ condition

In current cosmological models it is generally assumed that the power spectrum $`P(k)k`$ at small $`k`$. This is the famous Harrison-Zeldovich (HZ) condition. It is believed to describe the “primordial” fluctuations at very early times, the putative remnants of a period of “inflation” prior to the ordinary Big Bang phase pee93 ; padm . The HZ behaviour is appropriately cut-off at a short distance scale (i.e. for $`k`$ larger than some cut-off $`k_c`$). The type of the cut-off depends on the single cosmological model (e.g. Cold Dark Matter, Hot Dark Matter, etc.). The reason for the HZ condition is tied to considerations about the consistency of the Friedmann-Robertson-Walker (FRW) metrics under perturbations (see pee93 ; padm ): any other spectrum will give mass fluctuations which dominate over the homogeneous background either at some time in the future or past.

One can see prb  that in all the standard cosmological models the corresponding $`\stackrel{~}{\xi }(r)`$ is positive at small $`r`$ and has a negative tail $`r^4`$ at large $`r`$ with only one zero.

### 4.2 The power spectrum vs. the mass variance

Quite generally, by studying Eq. (3), it is not difficult to show prb  that for a spectrum $`P(k)k^n`$ for $`k0`$ (note that the Poisson-like system are given by $`n=0`$) and appropriately cut-off at large $`k`$, one has at large $`R`$

$$\sigma _M^2(R)\{\begin{array}{ccc}1/R^{3+n}\text{if}n<1\hfill & & \\ \mathrm{log}(R)/R^4\text{if}n=1\hfill & & \\ 1/R^4\text{if}n>1\hfill & & \end{array}$$
(5)
In terms of the non-normalized quantity $`\mathrm{\Delta }M^2(R)\sigma _M^2(R)R^6`$, Eq. (5) says that for $`n>0`$ (i.e. $`d^3r\stackrel{~}{\xi }(r)=0`$) we have a scaling behavior of $`\mathrm{\Delta }M^2`$ as a function of $`R`$ slower than for Poisson fluctuations, which correspond to $`n=0`$ (i.e. $`d^3r\stackrel{~}{\xi }(r)=c>0`$), with the limiting behavior for $`n1`$ corresponding to quadratic mass fluctuations which are proportional to the surface area of the sphere. These systems are thus characterised by surface fluctuations, ordered (or homogeneous enough, one could say) to give this very particular behaviour. The case $`n<0`$ (i.e. $`d^3r\stackrel{~}{\xi }(r)=+\mathrm{}`$) is typical of critical phenomena in which long range mainly positive correlations determine a “super-poissonian” behavior of quadratic fluctuations.

Equation (5) shows clearly that cosmological models, characterized by the HZ condition ($`n=1`$), are not only superhomogeneous, but are at the limit of the slowest possible behavior of quadratic mass fluctuations with the spatial scale.

## 5 Two “superhomogeneous” examples

In this section we present two examples of superhomogeneous point-particle distribution: the shuffled lattice prb  and the One Component Plasma (OCP) ocp .

1) In the case of particles placed on the sites of a regular cubic lattice (in any dimension $`d`$), one has prb  $`\rho _0=a^d`$ where $`a`$ is the lattice spacing and:

$$P_l(\stackrel{}{k})=(2\pi )^d\underset{\stackrel{}{H}0}{}\delta (\stackrel{}{k}\stackrel{}{H}),$$

where the sum is extended to all the vector of the reciprocal lattice but the origin $`\stackrel{}{H}=\stackrel{}{0}`$. Therefore we can say that around $`k=0`$ $`P_l(k)k^+\mathrm{}`$, and then it is the most superhomogenehous particle distribution with $`\sigma _M^2(R)R^{(d+1)}`$.

At this point let us introduce a stochastic displacement field, in which each particle is displaced from its lattice site $`\stackrel{}{R}`$ of a random and statistically isotropic vector $`\stackrel{}{u}_\stackrel{}{R}`$, each particle independently of the others. One can show preprint  that, if the probability density of the displacements $`p(u)`$ has a finite variance, then the resulting particle distribution has a power spectrum $`P_{sl}(k)k^2`$ at small $`k`$. If instead, the variance is infinite, and $`p(u)`$ has a tail going as $`u^a`$ (with $`d<a<d+2`$ for normalizability, but with infinite variance), then $`P_{sl}(k)k^{ad}`$ at small $`k`$. Therefore, choosing appropriately $`a`$ we can obtain any kind of superhomogeneous particle distribution with the exponent of $`P_{sl}(k)`$ at small $`k`$ ranging from $`0`$ to $`2`$. In particular for $`a=d+1`$ ($`a=4`$ in three dimension) we obtain the HZ condition.

2) Let us now present briefly another model taken from equilibrium statistical physics displaying superhomogeneous mass fluctuations: the OCP. It is a system of interacting point particles carrying a unit positive charge and mass repulsing each other via a Coulomb potential $`V(r)=1/r`$, in a continuous uniform background negatively charged giving overall charge neutrality ocp . One can show quite generally that, at thermodynamic equilibrium the system has only one phase, and at small $`k`$ the power spectrum of the particle distribution satisfies:

$$P(k)\frac{k^2}{4\pi \rho _0^2\beta },$$

where $`\rho _0`$ is the average density of particles and $`\beta `$ is the usual inverse temperature. It is very interesting to note that lebo  if, instead of having the Coulomb potential $`V(r)r^1`$, one has a potential $`V(r)r^2`$ one would obtain

$$P(k)\frac{k}{2\pi ^2n^2\beta },$$

that is the HZ condition. These results can be very important if applied to $`N`$-body simulations of gravitational dynamics to study structure formation starting from the primordial mass fluctuations thierry . In fact, in this context a great importance is given to the accurate preparation of the initial conditions for the simulation. For instance, this can be set through a modified OCP with $`V(r)r^2`$ at large scales and an appropriate modified behavior at smaller scales lebo .

Both 1) and 2) provide particle distributions similar to the right side picture of Fig. 1.