Yves Pomeau

Martine Le Berre

Pierre-Henri Chavanis

Bruno Denet

Supernovae: An example of complexity in the physics of compressible fluids

Keywords: 97, 60, Bw Supernovae -47, 27, ed Dynamical systems approaches

published or not. The documents may come

Introduction

It is a great pleasure to write this contribution in honor of Paul Manneville. We present below work belonging to the general field where he contributed so eminently, nonlinear effects in fluid mechanics. However, our topic is perhaps slightly unusual in this respect because it has to do with fluid mechanics on a grand scale, namely the scale of the Universe.

We all know that Astrophysics has to tackle a huge variety of phenomena, mixing widely scales of space and time. Our contribution below is perhaps the closest one can imagine of a problem of nonlinear and highly non trivial fluid mechanics in Astrophysics, the explosion of supernovae. In this fascinating field, many basic questions remain to be answered. The most basic one can be formulated as follows: stars evolve on very long time scales, in the billions years range, so why is it that some stars abruptly die (as far as their luminous output is concerned) in a matter of days or even of seconds? This was the point which motivated our work. Actually astrophysicists have worked since almost a century in the aim of finding plausible mechanisms in agreement with what has been observed, which is not the death process itself but the radiation emitted after the star death. Yet no definite answer is given, and it seems that the understanding of supernovae is "in an unsatisfactory state of affairs" to quote a recent review by Burrows [1]. The state of the art distinguishes two kinds of supernovae [2], corresponding to two kinds of death: explosion and collapse. The first kind concerns accreting white dwarfs (Type 1a) whose mass increases (and radius decreases because of gravity) to about the Chandrasekhar limit (the limit mass of a sphere of ultrarelativistic electrons when its radius tends to zero). The explosion which lasts about 10 seconds is thought to have a thermonuclear origin. It is followed by the expansion of a cloud of ejecta glowing brightly for many weeks as radioactive matter produced in the explosion decays, and leads to the total disintegration of the initial star without any remnants. The second kind concerns all the rest (young massive stars) which are supposed to collapse first and then emit matter and radiation, leading to the formation of remnants besides a neutron star (or a black hole).

We focus here on the second kind of star death, the collapse process which was understood and modeled via more and more sophisticated Newtonian and relativistic hydrodynamical simulations including weak interactions, nuclear equations of state, 3D effects, and new insights into the nucleosynthesis occurring before and during the collapse. While this is not our approach, let us present some general results of these studies. As a result of their evolution massive stars are supposed to pass through successive stages of hydrogen, helium, carbon, neon, oxygen and silicon fusion in their center, getting a layered structure like onions as progressively heavier atomic nuclei build up at the center. In this picture, although the fusion of hydrogen and helium takes millions of years, the last burning phase (silicon) lasts only two weeks, leading to a core of iron-group elements which undergoes catastrophic collapse as soon as its mass exceeds the Chandrasekhar value because the degeneracy pressure cannot support it. Theories predict emission of a huge quantity of neutrinos during the collapse which were supposed to play an important role in turning the collapse into the explosion seen with optical telescopes. Neutrinos were indeed observed once, three hours before the visible light from SN 1987A reached the Earth, at three separate neutrino observatories.

Our approach of the phenomenon of supernova explosion is not to try to describe quantitatively this immensely complex phenomenon, something which could well be beyond reach because it depends on so many uncontrolled and poorly known physical phenomena, like equations of state of matter in conditions not realizable in laboratory experiments, the definition of the initial conditions for the star collapse, the distribution of various nuclei in the star, etc. We argue that the huge difference of time scales involved should lead to a description of the supernovae in the light of catastrophe theory [START_REF] Thom | Stabilité structurelle et morphogénese, Benjamin[END_REF], the basic mechanism for star collapse being the loss of equilibrium between pressure and self-gravity. We consider the case where the star is in equilibrium during a long period, as the time dependent control parameter changes, then the series of equilibria presents a bifurcation such that no equilibrium exists beyond the critical point. The bifurcation should be a saddle-center because our system is Hamiltonian (it is described by the Euler-Poisson equations). In the slow regime we assume that an adiabatic approximation is possible, neglecting the diffusion and/or mechanical arrange-ment processes which could be responsible of some discontinuities in the route to equilibrium, because these processes have a short time scale compared to the slow evolution time.

While the theory of this equilibrium with the relevant equations is well-known, we did not find in the immense literature devoted to supernovae any interpretation of the transition from the slow evolution before the collapse to the fast collapse itself in the light of bifurcation theory. Such a description is however powerful because it provides the order of magnitude of the short time scale, as soon as the transition results from the sweeping of the bifurcation via a slow time-dependent parameter. In our description the slow growth of pressure and density is continuous, and the collapse involving the whole star does not happen because there is a drop of pressure, but because the equilibrium state no more exists. This is not what happens in the onion models where each time one fuel runs out, the star contracts because there is a drop of pressure, heats up, and burns the next one (the burning phase of hydrogen and helium taking millions of years, but the silicon burning in two weeks). Our aim was to follow the solution from the initial quasi-equilibrium state of the star, up to the post-collapse state, by sweeping the bifurcation, without changing the equation of state, in order to make clear that the collapse is an intrinsic process to any system which displays such a global loss of equilibrium. We study the solution close to the saddle-center, before the collapse, and beyond the singularity.

We introduce in section 2 a kind of universal equation valid for Hamiltonian systems (no dissipation) undergoing a saddle-center bifurcation (merging of a center and a saddle at the fold point) which is of the Painlevé I form. We set up the hypothesis required for this equation to describe the early stage of the loss of equilibrium in supernovae.

Using a particular equation of state, we show in section 3 that by a slow decrease of a given parameter (here the temperature), the series of equilibria do show a saddlecenter bifurcation.

In section 4 we study the approach towards the saddlecenter. We show that the full Euler-Poisson equations can be reduced to a normal form of the Painlevé I form valid at the first stage of the catastrophe, then we compare the numerical solution of the full Euler-Poisson equations with the solution of this universal equation.

Section 5 is devoted to the final stage of the collapse, just before the appearance of the singularity (divergence of the density and velocity). We show that the existence of a self-similar collapsing solution which agrees with the numerical simulations imposes that the gravity forces are stronger than the pressure ones, a situation which was not understood before. Usually the self-similar collapse, also called "homologous" collapse, is treated by assuming that pressure and gravity forces are of the same order that leads to scaling laws such as ρ ∼ r -α for the density with parameter α equal to 2. This corresponds to the Penston-Larson solution [START_REF] Penston | [END_REF]5]. Assuming that the gravity forces are larger than the pressure ones inside the core, we show first that a collapsing solution with α larger than 2 displays relevant asymptotic behavior in the outer part of the core, then we prove that it requires that α takes the value 24/11, which is larger than 2. We show that this result is actually in agreement with the numerical works of Penston (see Fig. 1 in [START_REF] Penston | [END_REF]) and Larson (see Fig. 1 in [5]) and many others (see Figs. [START_REF] Penston | [END_REF], 9, 10 and the first stage of Fig. 8, in [6]) and that this small discrepancy between α = 2 and α = 24/11 leads to non negligible consequences for the collapse characteristics. Contrary to the α = 2 case for which the velocity remains finite close to the center and tends to a constant supersonic value at large distances, our self-similar solution (in the sense of Zel'dovich) displays a velocity diverging at the center, and slowly vanishing as the boundary of the star is approached (an important point, see below).

Finally, in section 6, we describe the post-collapse dynamics without introducing any new ingredient in the physics. We point out that just at the collapse time, there is no mass in the center of the star, as in the case of the Bose-Einstein condensation [7,8]; the mass begins to accumulate in the inner core just after the singularity. Within the same frame as before (gravity forces dominant with respect to pressure ones), we derive the self-similar equations for the post-collapse regime and compare the solutions with a generalized version of the parametric free-fall solution proposed by Penston [START_REF] Penston | [END_REF].

In the present study we try to solve a simple model in a, what we believe, completely correct way. The interest of our model and analysis is that we fully explain the transition from the slow evolution before the collapse to the fast collapse itself. Continuing the evolution we observe and explain the occurrence of a new finite time singularity at the center, a singularity where the velocity field diverges. This singularity is not the standard homologous Penston-Larson collapse where all terms in the fluid equations are of the same order of magnitude. This singularity is of free-fall dynamics because the pressure force becomes (locally) negligible compared to the gravitational attraction 1 . The novelty of this solution is more than a mathematical nicety because the laws for this collapse, contrary to the ones of the homologous Penston-Larson collapse, are such that the velocity of infall tends to zero far from the center instead of tending to a constant supersonic value. This makes possible that the shock wave generated by the collapse escapes the center without the additional help of neutrinos as needed in models where the post-collapse initial conditions are a homologous Penston-Larson collapse far from the center. 1 Of course, the free-fall solution of a self-gravitating gas is well-known [START_REF] Penston | [END_REF]. However, it has been studied assuming either a purely homogeneous distribution of matter or an inhomogeneous distribution of matter behaving as ρ(r, t) -ρ(0, t) ∼ r 2 for r → 0, leading to a large distance decay ρ ∼ r -α with an exponent α = 12/7. We show that these assumptions are not relevant to our problem, and we consider for the first time a behavior ρ(r, t) -ρ(0, t) ∼ r 4 for r → 0, leading to the large distance decay with the exponent α = 24/11.

The Painlevé equation and the scaling laws

A supernova explosion lasts about ten seconds, when measured by the duration of the neutrino burst in SN1987A, and this follows a "slow" evolution over billions of years, giving an impressive 10 13 to 10 14 ratio of the slow to fast time scale. Such hugely different time scales make it a priori impossible to have the same numerical method for the slow and the fast dynamics. More generally it is a challenge to put in the same mathematical picture a dynamics with so widely different time scales. On the other hand the existence of such huge dimensionless numbers in a problem is an incentive to analyze it by using asymptotic methods.

Recently it has been shown [9] that such a slow-to-fast transition can be described as resulting from a slow sweeping across a saddle-node bifurcation. Let us recall that in such a bifurcation two fixed points, one stable the other unstable, merge and disappear when a parameter (usually supposed as constant) is changed. In [9] the transition is shown to be dynamical, because it occurs when a parameter changes slowly as a function of time. It means that the relevant parameter drifts in time until it crosses a critical value at the time of the catastrophe, this critical time being at the onset of saddle-node bifurcation for the dynamical system. Such a slow-to-fast transition is well known to show up in the van der Pol equation in the relaxation limit [START_REF] Dorodnicyn | [END_REF]. Interestingly, the analysis shows that this slowto-fast transition occurs on a time scale intermediate between the slow and long time scale, and that it is described by a universal equation solvable by the Riccati method. This concerns dynamical systems with dissipation, where the "universal equation" is first order in time. The supernovae likely belong to the class of dynamical catastrophes in our sense, because of the huge difference of time scales, but one must turn to a model of non dissipative dynamics if one assumes that the early post-bifurcation dynamics is described by inviscid fluid dynamics.

Such a dynamical model of catastrophes without dissipation and with time dependent sweeping across a bifurcation is developed below and applied to supernovae. We deal mostly with the early stage of the collapse, which we assume to be described by compressible fluid mechanics, without viscosity. Indeed the slow evolution of a star before the transition is a highly complex process not modeled in this approach because of the large difference in time scales: it is enough to assume that this slow evolution makes a parameter cross a critical value where a pair of equilibria (a saddle and a center) are merging. The universal equation describing the saddle-center bifurcation is the Painlevé I equation. We explain how to derive it from the fluid mechanical equations in the inviscid case, assumed to be valid for the interior of the star. Although applications of the ideas developed below could be found in more earthly situations like in subcritical bifurcation of Euler's Elastica with broken symmetry or the venerable Archimedes problem of (loss of) stability of floating bodies in an inviscid fluid [11], we shall refer below explicitly to the supernova case only. Our starting point is the following equation of Newtonian dynamics,

d 2 r 0 dt 2 = - ∂V ∂r 0 , ( 1
)
where r 0 can be seen as the radius of the star and V (r 0 , t) a time dependent potential. No mass multiplies the acceleration, which is always possible by rescaling the potential V (.). We shall derive later this equation for an inviscid compressible fluid with gravitation and an equation of state changing slowly as a function of time, for a radially symmetric geometry and with a finite mass. Contrary to the case studied in [9], this equation is second order in time because one neglects dissipation compared to inertia. The potential V (.) on the right-hand side represents the potential energy of the star, with the contributions of gravity and of internal energy [START_REF] Landau | Statistical physics[END_REF]. At equilibrium the right-hand side is zero. Given the potential V (.) this depends on two parameters (linked to the total mass and energy), r 0 and another physical parameter which may be seen as the temperature. Because of the long term evolution of the star interior by nuclear reactions and radiation to the outside, its temperature changes slowly. We shall assume that this slow change of parameter makes the equilibrium solution disappear by a saddle-center bifurcation when the temperature T crosses a critical value. We note that a saddle-node bifurcation is sometimes called turning, or tipping point instability, whereas the word "saddle-node" (noeud-col in french) was coined by H. Poincaré in his Ph.D. thesis. Such a bifurcation is a fairly standard problem treated by Emden [START_REF] Emden | Gaskugeln anwendungen der mechanischen wärmetheorie auf kosmologische und meteorologie probleme[END_REF] for a selfgravitating gas at finite (and changing, but not as function of time) temperature in a spherical box. It was also discussed by Ebert [START_REF] Ebert | [END_REF], Bonnor [15], and McCrea [16] by varying the pressure, and by Antonov [17] and Lynden-Bell and Wood [18] by varying the energy. See Chavanis [19,20] for recent studies. A saddle node is also present in the caloric curve of self-gravitating fermions at finite temperature which has the form of a "dinosaur's neck" [25]. A saddle-center may be worked out in the mass-radius relation of neutron stars determined by Oppenheimer and Volkoff [21] when the mass crosses a critical value M OV (see also section 109 of [START_REF] Landau | Statistical physics[END_REF], figure 52) and in the massradius relation of boson stars [22][23][24].

As we do not solve the energy equation, the parameter T could be any parameter describing the smooth changes of the star interior prior to the fast transition. Following the ideas of reference [9] we look for a finite change in the system on a time scale much shorter than the time scale of the control parameter (here the temperature T ). Two time scales are involved: the long time scale of evolution of T , denoted as θ below, and the short time scale τ which is the fundamental period of a pressure oscillation in the star. Our approach will show that the early stage of the collapse is on a time scale intermediate between the fast and slow scale and give a precise definition of the initial conditions for the fast process.

Let us expand the potential V (.) in Poincaré normal form near the saddle-center bifurcation: In the expression above, R, a relative displacement, can be seen as the difference between r c , the value of the radius of the star at the saddle-center bifurcation and its actual value, R = (r 0 -r c )/r c , a quantity which decreases as time increases, because we describe the collapse of the star. Actually the quantity R will be seen later as the Lagrangian radial coordinate, a function depending on r, the radial distance. The saddle-center bifurcation is when the -now time dependent -coefficient a of equation ( 2) crosses 0. Setting to zero the time of this crossing, one writes a = -ct, where c, a constant, is small because the evolution of V is slow. This linear time dependence is an approximation because a(t) is, in general, a more complex function of t than a simple ramp. However, near the transition, one can limit oneself to this first term in the Taylor expansion of a(t) with respect to t, because the transition one is interested in takes place on time scales much shorter than the typical time of change of a(t). Limiting oneself to displacements small compared to r c , one can keep in V (R) terms which are linear and cubic (the coefficient b is assumed positive) with respect to R because the quadratic term vanishes at the saddle-center transition (the formal statement equivalent to this lack of quadratic term in this Taylor expansion of V (R) is the existence of a non trivial solution of the linearized equation at the bifurcation). Moreover higher order terms in the Taylor expansion of V (.) near R = 0 are neglected in this analysis because they are negligible with the scaling law to be found for the magnitude of R near the transition. This is true at least until a well defined time where the solution has to be matched with the one of another dynamical problem, valid for finite R. At t = 0, the potential V (.) is a cubic function of R, exactly the local shape of a potential in a metastable state. For a and b positive, the potential has two extrema, one corresponding to a stable equilibrium point at R = √ a/b and one unstable at R = -√ a/b. In the time dependent case, the potential evolves as shown in Fig. 1 and the equations ( 1)-(2) become

V = -aR + b 3 R 3 + ..., (2)
d 2 R dt 2 = R = -ct -bR 2 , ( 3
)
where the parameter c is supposed to be positive, so that the solution at large negative time is close to equilibrium and positive, crosses zero at a time close to zero and diverges at finite positive time.

To show that the time scale for the dynamical saddlecenter bifurcation is intermediate between the long time scale of the evolution of the potential V (.) and the short time scale of the pressure wave in the star, let us derive explicitly these two relevant short and long time scales. For large negative time the solution of equation ( 3) is assumed to evolve very slowly such that the left-hand side can be set to zero. It gives

R(t) ≃ √ c b (-t) (4)
which defines the long time scale as θ = b/c (recall that R, a relative displacement scaled to the star radius r c , has no physical scale).

As for the short time scale, it appears close to the time t = t * where the solution of equation (3) tends to minus infinity. In this domain the first term in the right-hand side is negligible with respect to the second one, the equation reduces to R = -bR 2 , which has the characteristic time τ = 1/ √ b. Let us scale out the two parameters b, c of equation [START_REF] Thom | Stabilité structurelle et morphogénese, Benjamin[END_REF]. Defining R = R/r s and t = t/t 0 the original equation takes the scaled form

d 2 R d t2 = -t -R2 , (5)
when setting c = r s /t 3 0 and b = 1/(r s t 2 0 ). Inversely, t 0 = 1/(bc) 1/5 and r s = c 2/5 /b 3/5 . The solution of equation ( 5) is called the first Painlevé transcendent, and cannot be reduced to elementary functions [26].

A non trivial question is the sign of b and c. which cannot be decided in a simple way from the general features of the model. To give sense to this theory, the bifurcation must occur in the forward time direction, namely b and c must have the same sign for this theory to make sense. This sign common to b and c determines if the radius R grows or decays just after the bifurcation. It is commonly assumed that, in the early stage of supernova explosion, the star collapses onto itself, and then that this inward motion is reversed to yield the observed explosion, the mechanism of reversal being not well understood yet. In the model to be studied later, the radius actually decreases after the bifurcation, as believed to occur in the early stage of supernovae explosions. But it is well possible that in other models of equation of state with a saddlecenter bifurcation, the early stage and perhaps even all the dynamics triggered by the bifurcation is an outward motion, without any phase of collapse.

The writing of the Painlevé equation in its parameter free form yields the characteristic time scale t 0 of equation (3) in terms of the short and long times,

t 0 = (θτ 4 ) 1/5 . ( 6
)
This intermediate time is such that τ ≪ t 0 ≪ θ; it could be of the order of several hours when taking θ ∼ one billion years, τ ∼ 10 sec. The corresponding spatial extension R is of order

r s = ( τ θ ) 2/5 , (7)
much smaller than unity. The one-fifth power in equations ( 6) and ( 7) is "typical" of the Painlevé I equation, which has a symmetry expressed in terms of the complex fifth root of unity.

To solve equation (5) we have to define the initial conditions. Choosing the initial conditions at large negative time t i , we may assume that the asymptotic relation ( 4) is fulfilled at this time, that gives,

   R( ti ) = √ -ti , Ṙ( ti ) = -1 2 √ -ti . ( 8
)
The numerical solution of equation ( 5) is drawn in Fig. 2 leading to a finite time singularity. With the initial conditions (8) the solution is a non oscillating function (blue curve) diverging at a finite time t * ≃ 3.4 (note that the divergence is not yet reached in Fig. 2). But we may assume that, at very large negative time, the initial conditions slightly differ from the asymptotic quasi-equilibrium value (8). In that case the solution displays oscillations of increasing amplitude and period as time increases, in agreement with a WKB solution of the linearized problem. Let us put R( t) ≈

√

-t+δ R, δ R small which satisfies the linear equation

δ R = -2 √ -tδ R. (9)
A WKB solution, valid for (-t) very large is

δ R = ∑ ± c ± (-t) -1/4 e ±i 4 √ 2 5 (-t) 5/4 . ( 10
)
It represents oscillations in the bottom of the potential V ( R, t) = t R + R3 /3 near R = √ -t. The two complex conjugate coefficients c ± defining the amplitudes are arbitrary and depend on two real numbers. Therefore, the cancelation of the oscillations defines uniquely a solution of the Painlevé I equation. This is illustrated in Fig. 2 where the blue curve has no oscillation (see above) while the red curve displays oscillations of increasing period and a shift of the divergence time.

Near the singularity, namely just before time t = t * , the dominant term on the right-hand side of equation ( 5) is R2 so that R becomes approximately R(t) ≃ -6/( t * -t) 2 , or in terms of the original variables R and t,

R(t) ≃ -6r s ( t 0 t * -t ) 2 . ( 11
)
This behavior will be compared later to the full Euler-Poisson model (see Fig. 16 and relative discussion). Note that this divergence is completely due to the nonlinearity, and has little to do with a linear instability. The applicability of this theory requires R ≪ 1, because it relies on the Taylor expansion of V (.) in equation ( 1) near

r 0 = r c . It is valid if |t -t * | ≫ τ .
Therefore the collapse (we mean by collapse the very fast dynamics following the saddlecenter bifurcation) can be defined within a time interval of order τ , the center of this interval being the time where the solution of equation (5) diverges, not the time where the linear term in the same equation changes sign. Moreover the duration of the early stage of the collapse is, physically, of order (θτ 4 ) 1/5 , much shorter than the time scale of evolution of the temperature, but much longer than the elastic reaction of the star interior.

The blow-up of the solution of equation (3) at finite time does not imply a physical singularity at this instant. It only shows that, when t approaches t * by negative values, R(t) grows enough to reach an order of magnitude, here the radius of the star, such that the approximation of V by the first two terms (linear and cubic with respect to R) of its Taylor expansion is no longer valid, imposing to switch to a theory valid for finite displacements. In this case, it means that one has to solve, one way or another, the full equations of inviscid hydrodynamics, something considered in section 3. A warning at this stage is necessary: we have to consider more than one type of finite time singularity in this problem. Here we have met first a singularity of the solution of the Painlevé I equation, a singularity due to various approximations made for the full equations which disappear when the full system of Euler-Poisson equations is considered. But, as we shall see, the solution of this Euler-Poisson set of dynamical equations shows a finite time singularity also, which is studied in section 5 and which is related directly to the supernova explosion.

Below we assume exact spherical symmetry, although non spherical stars could be quite different. A given star being likely not exactly spherically symmetric, the exact time t * is not so well defined at the accuracy of the short time scale τ because it depends on small oscillations of the star interior prior to the singularity (the amplitude of those oscillations depends on the constants c ± in the WKB part of the solution, and the time t * of the singularity depends on this amplitude). One can expect those oscillations to have some randomness in space and so not to be purely radial. The induced loss of sphericity at the time of the collapse could explain the observed expulsion of the central core of supernovae with large velocities, up to 500 km per second [START_REF] Hoflich | Cosmic Explosions in Three Dimensions: Asymmetries in Supernovae and Gamma Ray Bursts[END_REF] a very large speed which requires large deviations to sphericity. However there is an argument against a too large loss of sphericity: the time scale t 0 for the part of the collapse described by the Painlevé equation is much longer than τ , the typical time scale for the evolution of the inside of the star. Therefore one may expect that during a time of order t 0 , the azimuthal heterogeneities are averaged, restoring spherical symmetry on average on the longer time scale t 0 . However this does not apply if the star is intrinsically non spherically symmetric because of its rotation.

Within this assumption of given slow dependence with respect to a parameter called T , we shall derive the dynamical equation (3) from the fluid equations with a general pressure-density relation and the gravity included. To streamline equations and explanations, we shall not consider the constraint of conservation of energy (relevant on the fast time scale).

3 Euler-Poisson system for a barotropic star presenting a saddle-center

Barotropic Euler-Poisson system

We shall assume that the star can be described as a compressible inviscid fluid with a barotropic equation of state p = p(ρ). The relevant set of hydrodynamic equations are the barotropic Euler-Poisson system. These are dynamical equations for a compressible inviscid fluid with a pressuredensity relation, including the gravitational interaction via Poisson equation. Note that there is no dynamical equation for the transport of energy. They read

∂ρ ∂t + ∇ • (ρu) = 0, (12)
ρ

[ ∂u ∂t + (u • ∇)u ] = -∇p -ρ∇Φ, ( 13
)
∆Φ = 4πGρ, ( 14
)
where u is the fluid velocity vector, ρ the mass density, and G Newton's constant. Using the equation of continuity [START_REF] Landau | Statistical physics[END_REF], the momentum equation ( 13) may be rewritten as

∂ ∂t (ρu) + ∇(ρu ⊗ u) = -∇p -ρ∇Φ. ( 15
)
The potential energy of this self-gravitating fluid is V = U + W where

U = ∫ ρ ∫ ρ p(ρ ′ ) ρ ′ 2 dρ ′ dr, ( 16
)
is the internal energy and

W = 1 2 ∫ ρΦ dr, ( 17
)
is the gravitational energy. The internal energy can be written as

U = ∫ [ρh(ρ) -p(ρ)] dr = ∫ H(ρ)
dr where we have introduced the enthalpy h(ρ), satisfying dh(ρ) = dp(ρ)/ρ, and its primitive H(ρ) = ∫ ρ 0 h(ρ)dρ.

Hydrostatic equilibrium and neutral mode

In this section we briefly recall different formulations of the equilibrium state of a self-gravitating gas. From equation ( 13), the condition of hydrostatic equilibrium writes

∇p + ρ∇Φ = 0. ( 18
)
Dividing this equation by ρ, taking the divergence of the resulting expression, using Poisson equation ( 14), and recalling that p = p(ρ) for a barotropic gas, we obtain a differential equation for ρ that is

∇ • [ p ′ (ρ) ρ ∇ρ ] + 4πGρ = 0. ( 19
)
For a barotropic equation of state by definition p = p(ρ).

The condition of hydrostatic equilibrium ( 18) implies ρ = ρ(Φ). Substituting this relation in Poisson equation ( 14), we obtain a differential equation for Φ that is

∆Φ = 4πGρ(Φ). ( 20
)
Introducing the enthalpy, satisfying ∇h = ∇p/ρ, the condition of hydrostatic equilibrium (18) can be rewritten as

∇h + ∇Φ = 0. (21)
Therefore, at equilibrium, h(r) = -Φ(r) + C where C is a constant. Since the gas is barotropic, we also have ρ = ρ(h). Taking the divergence of equation ( 21) and using Poisson equation ( 14), we obtain a differential equation for h that is

∆h + 4πGρ(h) = 0. ( 22
)
These different formulations are equivalent. In the following, we will solve the differential equation (22).

To determine the dynamical stability of a steady state of the Euler-Poisson system ( 12)-( 14), we consider a small perturbation about that state and write f (r, t) = f (r) + δf (r, t) for f = (ρ, u, Φ) with δf (r, t) ≪ f (r). Linearizing the Euler-Poisson system about that state, and writing the perturbation as δf (r, t) ∝ e λt , we obtain the eigenvalue equation

λ 2 δρ = ∇ • [ρ(∇δh + ∇δΦ)] . ( 23
)
The neutral mode (λ = 0) which usually signals the change of stability is the solution of the differential equation

∇δh + ∇δΦ = 0. ( 24
)
Taking the divergence of this equation and using Poisson equation ( 14), it can be rewritten as

∆δh + 4πGρ ′ (h)δh = 0. ( 25
)
This equation may also be written in terms of δρ by using δh = p ′ (ρ)δρ/ρ. We get

∆ ( p ′ (ρ) ρ δρ ) + 4πGδρ = 0. ( 26
)
In the following, we will solve the differential equation ( 25).

An isothermal equation of state with a polytropic envelope implying a saddle node

The series of equilibria of an isothermal self-gravitating gas with p = ρT is known to present a saddle node [START_REF] Emden | Gaskugeln anwendungen der mechanischen wärmetheorie auf kosmologische und meteorologie probleme[END_REF]19]. Therefore a self-gravitating isothermal gas is a good candidate for our investigation. However, it has the undesirable feature to possess an infinite mass because its density decreases too slowly (as r -2 ) at large distances. Therefore, to have a finite mass, it must be confined artificially into a "box". In order to skip this difficulty, we propose to use here an equation of state that is isothermal at high densities and polytropic at low densities, the polytropic equation of state serving as an envelope that confines the system in a finite region of space without artificial container. Specifically, we consider the equation of state

2 p(ρ) = ρ * T ( √ 1 + ρ/ρ * -1 ) 2 . ( 27
)
For ρ → +∞, it reduces to the isothermal equation of state p = ρT . For ρ → 0, it reduces to the polytropic equation of state p = Kρ 2 with polytropic index γ = 2 and polytropic constant K = T /(4ρ * ). The enthalpy function h(ρ) defined by dh = dp/ρ is explicitly given by h(ρ) = 2T ln

( 1 + √ 1 + ρ/ρ * ) -2T ln(2), (28)
where the constant of integration has been determined such that h(ρ = 0) = 0. With this choice, the enthalpy vanishes at the edge of the star. The inverse relation writes

ρ(h) = 4ρ * ( e h/T -e h/2T ) . ( 29
)
In the following, it will be convenient to use dimensionless variables. The parameters regarded as fixed are ρ * , M , and G. From ρ * and M we can construct a length L * = (M/ρ * ) 1/3 . Then, we introduce the dimensionless quantities

ρ = ρ ρ * , r = r L , Φ = Φ Gρ * L 2 . ( 30
)
and

T = T Gρ * L 2 , p = p GL 2 ρ 2 * , t = t √ Gρ * . ( 31
)
Working with the dimensionless variables with tildes amounts to taking G = ρ * = M = 1 in the initial equations, a choice that we shall make in the following.

Equilibrium solution and temperature-radius relation

The equilibrium solution is obtained by solving equation (22) with equation ( 29). Using the dimensionless variables defined in Sec. 3.3, assuming spherical symmetry, and setting r = r/ √ T , ĥ = h/T , Φ = Φ/T , ρ = ρ, and M = M/T 3/2 , we obtain ĥ,r 2 + 2 r ĥ,r + 4π ρ( ĥ) = 0,

where ρ( ĥ) = 4

( e ĥ -e ĥ/2

) . (

Using Gauss theorem Φ ,r = M (r)/r 2 , where

M (r) = ∫ r 0 ρ(r ′ )4πr ′ 2 dr ′ , ( 34
)
is the mass profile, and the equilibrium relation Φ ,r = -h ,r , we obtain Φ,r = -ĥ,r = M (r)/r 2 that allows us to determine the mass profile from the enthalpy profile using 3 M (r) = -r 2 ĥ,r . (

The boundary conditions of equation ( 32) at r = 0 are ĥ(0) = ĥ0 and ĥ,r (0) = 0. For a given value of ĥ0 , the smallest root of ĥ(r), which is also the one of ρ(r), see Figs. 3 and4, defines the normalized radius r0 of the star. The radius r 0 of the star is therefore r 0 = √ T r0 . On the other hand, Gauss theorem applied at the surface of the star where M = 1 (i.e. M0 = 1/T 3/2 ) leads to ĥ,r (r 0 ) = -1/( √ T r 2 0 ). From these equations, we obtain 4 r 0 = ( r0 -ĥ,r (r 0 )

) 1/3 , T = 1 ( -r 2 0 ĥ,r (r 0 ) ) 2/3 . ( 36
)
3 Equation ( 35) may also be obtained by multiplying equation (32) by r2 and integrating between 0 and r. 4 We can come back to the original (dimensional) variables by making the substitution R → R/L = Rρ 32) and (37), radial profile of the enthalpy ĥ(r) (solid red curve) and neutral mode j(r) (dashed blue curve) for ĥ0 = 2.296 corresponding to the saddle-center, point A of Fig. 5.

1/3 * /M 1/3 and T → T /(Gρ * L 2 ) = T /(Gρ 1/3 * M 2/3 ).
The solution of equation (32), drawn in Fig. 4 solid line, has a single free parameter ĥ0 since its Taylor expansion near r = 0 is like ĥ = ĥ0 + h 2 r2 + ... with ĥ0 free, h 2 = -2π 3 ρ( ĥ0 ), and so on for the higher order coefficients. By varying ĥ0 from 0 to +∞ we can obtain the whole series of equilibria r 0 (T ) giving the radius of the star as a function of the temperature, using the quantities ĥ0 (or r0 ) as a parameter. The result is a spiralling curve shown in Fig. 5 where only the upper part is stable, the solution loosing its stability at the saddle-center (turning point A), as studied in the next subsection 5 . The saddle-center is found nu-merically to occur at ĥ0 = 2.296.., or ρ0 = 27.1299.., that leads to the following critical values for the mass, temperature and radius respectively, Mc = 0.52, T c = 1.546... and rc = 0.385.. (hence r c = √ T c rc = 0.479...). The center of the spiral is obtained for ĥ0 → ∞. 32)-( 36

) (increasing the input parameter ĥ0).

There is a saddle-center bifurcation when equation ( 32) linearized about the profile ĥ(r) determined previously has a non trivial solution. This corresponds to the neutral mode δh defined by the unscaled equation (25). In terms of the scaled variables this linearized equation reads

Ω[j(r)] = j ,r 2 + 2 r j ,r + 4π dρ d ĥ j(r) = 0, ( 37
)
where Ω is a linear operator acting on function j of r. Let us precise that we have the following boundary conditions: arbitrary j(0) and j ′ (0) = 0. Furthermore, we automatically have j ′ (r 0 ) = 0 since δM (r 0 ) = 0. The neutral mode j(r), valid at the critical temperature T c , is pictured in Fig. 4, dashed blue line. We consider below the dynamics of the function M (r, t) which is the mass contained inside the sphere of radius r in the star.

Dynamics close to the saddle-center: derivation of Painlevé I equation

In this section we show that the dynamics close to the saddle-center reduces to Painlevé I equation. This property will be proved first by showing that the normal form of the full Euler-Poisson system ( 12)-( 14) is of Painlevé I form, secondly by comparing the normal form solutions to the full Euler-Poisson ones derived by using a numerical package for high-resolution central schemes [29].

polytrope n = 1 that is known to be stable, we conclude that the upper branch in Fig. 5 is stable up to the turning point A.

Then, the series of equilibria loses a mode of stability at each turning point of temperature T and becomes more and more unstable.

Simplification of the hydrodynamic equations close to the saddle-center

We now consider the dynamical evolution of the star, in particular its gravitational collapse when the temperature falls below T c . In this section and in the following one we use a simplified model where advection has been neglected, an approximation valid in the first stage of the collapse only. In the following we restrict ourselves to spherically symmetric cases, likely an approximation in all cases, and certainly not a good starting point if rotation is present. However this allows a rather detailed analysis without, hopefully, forgetting anything essential. Defining u as the radial component of the velocity, let us estimate the order of magnitude of the various terms in Euler's equations during the early stage of the collapse, namely when equation (3) is valid (this assuming that it can be derived from the fluid equations, as done below). The order of magnitude of u ,t is the one of R, that is Ṙ/t 0 , with t 0 the characteristic time defined by equation ( 6). The order of magnitude of the advection term uu ,r is Ṙ2 /r 0 (here R is dimensional), because one assumes (and will show) that the perturbation during this early stage extends all over the star. Therefore uu ,r ∼ u ,t (R/r 0 ) is smaller than u ,t by a factor R/r 0 , which is the small a-dimensional characteristic length scale defined by the relation (7). Neglecting the advection term in equations ( 13) and (15) gives

∂ ∂t (ρu) = ρ ∂ ∂t u = -∇p -ρ∇Φ. ( 38
)
In the spherically symmetric case it becomes

u ,t = - 1 ρ p ,r - 4πG r 2 ∫ r 0 dr ′ r ′2 ρ(r ′ , t), ( 39
)
where we used Gauss theorem

Φ ,r = 4πG r 2 ∫ r 0 dr ′ r ′2 ρ(r ′ , t), ( 40
)
derived from Poisson equation [START_REF] Ebert | [END_REF]. Taking the divergence of the integro-differential dynamical equation (39) allowing to get rid of the integral term, we obtain

( 2 r u + u ,r ) ,t = - ( h ,r 2 + 2 r h ,r + 4πGρ(h) ) , ( 41
)
which is the dynamical equation for the velocity field. This equation has been derived from the Euler-Poisson system (12)-( 14) where the advection has been neglected, that is valid during the time interval of order t 0 before the critical time. To derive the Painlevé I equation from the dynamical equation (41) we consider its right-hand side as a function of ρ with an equation of state of the form p(ρ) = ρ * T f (ρ/ρ * ) depending on a slow parameter T , and we expand the solution near a saddle-center bifurcation which exists when there is more than one steady solution of equation (41) for a given total mass M = 4π ∫ ∞ 0 dr ′ r ′2 ρ(r ′ ) and temperature T , two solutions merging and disappearing as the temperature crosses a critical value T c . This occurs for the equation of state defined by equation [START_REF] Hoflich | Cosmic Explosions in Three Dimensions: Asymmetries in Supernovae and Gamma Ray Bursts[END_REF], see Fig. 5 where a saddle-center exists at point A. Although this formulation in terms of the velocity field u(r, t) is closely related to the heuristic description developed in Sec. 2, in the following we find it more convenient to work in terms of the mass profile M (r, t). Obviously the two formulations are equivalent.

The equation for the mass profile M (r, t)

In view of studying the dynamics of the solution close to the saddle-center, let us assume a slow decrease of the temperature versus time, of the form T = T c (1 -γ ′ t) with positive γ ′ in order to start at negative time from an equilibrium state. Taking the time derivative of the equation of continuity [START_REF] Landau | Statistical physics[END_REF] and using equation (38), we get the two coupled equations 6

∂ 2 ρ ∂t 2 = ∇ • (∇p + ρ∇Φ), ( 42
)
∆Φ = 4πGρ. ( 43
)
According to the arguments given in Sec. 4.1, these equations are valid close to the saddle-center during the early stage of the collapse 7 . By contrast, when we are deep in the collapse regime (see Secs. 5 and 6) the advection term is important and we must come back to the full Euler-Poisson system ( 12)-( 14).

In the following, we use the dimensionless variables of Sec. 3.3. In the spherically symmetric case, using Gauss theorem (40), the system (42)-(43) writes

∂ 2 ρ ∂t 2 = 1 r 2 [ r 2 p ,r + ρ ∫ r 0 dr ′ 4πr ′2 ρ(r ′ ) ] ,r . ( 44
)
It has to be completed by the boundary conditions imposing zero mass at the center of the star, and a constant total mass

∫ r0 0 dr ′ 4πr ′2 ρ(r ′ , t) = 1, (45)
where r 0 is the star radius (practically the smallest root of ρ(r) = 0). Let us define the variable

M (r, t) = ∫ r 0 dr ′ 4πr ′2 ρ(r ′ , t) (46)
which represents the mass of fluid contained inside a sphere of radius r at time t. Multiplying the two sides of equation (44) by 4πr 2 , and integrating them with respect to 6 These equations are similar to the Smoluchowski-Poisson system (describing self-gravitating Brownian particles in the strong friction limit) studied in [START_REF] Chavanis | [END_REF] except that it is second order in time instead of first order in time. 7 These equations are also valid for small perturbations about an equilibrium state since we can neglect the advection term u • ∇u at linear order.

the radius, we obtain the dynamical equation for the mass profile M (r, t),

∂ 2 M (r, t) ∂t 2 = 4πr 2 p ,r + 1 r 2 M ,r M, ( 47
)
where the term p ,r = p ′ (ρ)ρ ,r has to be expressed as a function of ρ(r, t) = 1 4πr 2 M ,r and ρ ,r (r, t)

= 1 4πr 2 (M ,r 2 - 2 r M ,r
). Using the relation ( 27), one has

p ′ (ρ) = T ( 1 - 1 √ 1 + ρ ) . ( 48
)
The first term of equation (47) becomes

4πr 2 p ,r = T L(M )g(M ,r ) (49) with { L(M ) = M ,r 2 -2 r M ,r g(M ,r ) = 1 - 1 √ 1+ 1 4πr 2 M,r . ( 50
)
Introducing this expression into equation ( 47), the dynamical equation for M (r, t) writes

∂ 2 M (r, t) ∂t 2 = T L(M )g(M ,r ) + 1 r 2 M ,r M. ( 51
)
The boundary conditions to be satisfied are { M (0, t) = 0

M (r 0 (t), t) = 1 = 4π ∫ r0(t) 0 dr ′ r ′2 ρ(r ′ , t). ( 52
)
In the latter relation the radius of the star r 0 (t) depends on time. However this dependance will be neglected below, see equation ( 68), because we ultimately find that the star collapses, therefore its radius will decrease, leading to r 0 (t) < r c , or M (r 0 (t), t) = M (r c ) as time goes on.

Equilibrium state and neutral mode

A steady solution of equation ( 51) is determined by

T L(M )g(M ,r ) + 1 r 2 M ,r M = 0. ( 53
)
Using Gauss theorem Φ ,r = M (r)/r 2 , and the equilibrium relation Φ ,r = -h ,r , we can easily check that equation (53) is equivalent to equation (32). We now consider a small perturbation about a steady state and write M (r, t) = M (r) + δM (r, t) with δM (r, t) ≪ M (r). Linearizing equation (51) about this steady state and writing the perturbation as δM (r, t) ∝ e λt , we obtain the eigenvalue equation

λ 2 δM = T [L(δM )g(M ,r ) + L(M )g ′ (M ,r )δM ,r ] + 1 r 2 (M δM ) ,r . ( 54
)
The neutral mode, corresponding to λ = 0, is determined by the differential equation

T [L(δM )g(M ,r ) + L(M )g ′ (M ,r )δM ,r ] + 1 r 2 (M δM ) ,r = 0.
(55) Using Gauss theorem δΦ ,r = δM (r)/r 2 , and the relation δΦ ,r = -δh ,r satisfied at the neutral point (see Sec. 3.2), we can check that equation ( 55) is equivalent to equation (37). This implies that the neutral mass profile is given by δM (r) = -r 2 j ,r .

(56)

Normal form of the mass profile M (r, t)

The derivation of the normal form close to the saddlecenter proceeds mainly along the lines of [START_REF] Chavanis | [END_REF] 8 . The mass profile is expanded as

M (r, t) = M (c) (r) + ϵM (1) (r, t) + ϵ 2 M (2) (r, t) + ... ( 57
)
where M (c) (r) is the equilibrium profile at T = T c (see above) drawn in solid line in Fig. 6, and ϵ is a small parameter which characterizes a variation of the temperature with respect to its value at the collapse. We set

T = T c (1 -ϵ 2 T (2) ), (58)
which amounts to defining ϵ 2 T (2) = γ ′ t, and rescaling the time as t = t ′ /ϵ 1/2 (this implies that γ ′ ∼ ϵ 5/2 is a small quantity). Substituting the expansion (57) into equation (51), we get at leading order the equilibrium relation

T c L (c) g (c) + 1 r 2 M (c) ,r M (c) = 0, (59)
which has to satisfy the boundary conditions

M (c) (0) = M (c) ,r (0) = 0; M (c) (r c ) = 1. ( 60
)
To order 1 we have

T c ( L (1) g (c) + L (c) g ′(c) M (1)
,r

) + 1 r 2 (M (1) M (c)
) ,r = 0, (61) and to order 2

∂ 2 M (1) (r, t ′ ) ∂t ′2 = T c F (2) + 1 r 2 [ (M (2) M (c) ) ,r + M (1) M (1)
,r

] ,

(62) where

F (2) = F (2) 1 + F (2) 2 + F (2) 3 (63)
with

F (2) 1 = ( L (2) -T (2) L (c) ) g (c) , ( 64
)
F (2) 2 = L (1) g ′(c) M (1) ,r , (65)

M ζ

Fig. 6. Mass M (solid blue line) inside the star versus the radial variable r at the saddle-center, solution of equations ( 32)-( 34) for T = Tc, i.e. ĥ0 = 2.296. The dashed red line is for ζ(r), solution of equation ( 83) with appropriate initial conditions for solving the adjoint problem (in this caption, we have restored the "hat" on the variables).

F (2) 3 = L (c) [ g ′(c) M (2) ,r + g ′′(c) 2 (M (1) ,r ) 2 ] , ( 66
)
where

L (c) = L(M (c) ), L (n) = L(M (n) ), g (c) = g(M (c) ,r ), g ′(c) = ( dg dM,r ) (c) and g ′′(c) = ( d 2 g dM 2 ,r
) (c) . The r-dependent quantities can be written in terms of the equilibrium density function ρ (c) (r) as

           L (c) = 4πr 2 ρ (c) ,r , g (c) = 1 - 1 √ 1+ρ (c) , g ′(c) = 1 8πr 2 (1+ρ (c) ) 3/2 , g ′′(c) = - 3 4(4πr 2 ) 2 (1+ρ (c) ) 5/2 . ( 67
)
The boundary conditions are

{ M (n) (0, t ′ ) = 0; M (n) ,r (0, t ′ ) = 0; M (n) (r c , t ′ ) = 0. ( 68
)
Let us rescale the quantities in equations ( 44)-(68) by using the critical value T c for the temperature in the rescaled variables. We thus define T = T /T c , r = r/ √ T c , t = t, M = M/T 3/2 c , ĥ = h/T c , and ρ = ρ. This rescaling leads to the same expressions as the unscaled ones in equations (44)-( 68), except that T c is set to one. Furthermore, at the critical point, the rescaled variables coincide with those introduced in Sec. 3.4. In the following, we drop the superscripts to simplify the notations.

The foregoing equations have a clear interpretation. At zeroth order, equation (59) corresponds to the equilibrium state (53), equivalent to equation (32), at the critical point T c . The critical mass profile is drawn in Fig. 6 solid line. At order 1, equation (61) has the same form as the differential equation (55), equivalent to equation (37), determining the neutral mode (corresponding to the critical point). Because equation ( 61) is linear, its solution is

M (1) (r, t ′ ) = A (1) (t ′ )F (r), (69)
where

F (r) = δM (r) = -r 2 j ,r , (70)
according to equation (56). This solution, drawn in Fig. 7-(a), thick black line, fulfills the boundary conditions (68).

The corresponding density profile ρ (1) (r, t ′ ) = A (1) (t ′ )δρ(r) is drawn in Fig. 7-(b), where

δρ(r) = F ,r 4πr 2 = j(r) ( dρ dh ) (c) . ( 71
)
At order 2, equation (62) becomes

F (r) Ä(1) (t ′ ) = -T (2) L (c) g (c) + D(F )A (1)2 + C(M (2) ), (72)
where

D(F ) = 1 r 2 F F ,r + 1 2 L (c) g ′′(c) F 2 ,r + g ′(c) L(F )F ,r , ( 73
)
and

C(M (2) ) = L (2) g (c) + 1 r 2 (M (2) M (c) ) ,r + L (c) g ′(c) M (2) ,r .
(74) To write the dynamical equation for A(t) (1) in a normal form, we multiply equation (72) by a function ζ(r) and integrate over r for 0 < r < r c , where r c is the radius of the star at T = T c . We are going to derive the function ζ(r) so that the term C(M (2) ) disappears after integration (see Appendix A for details about the boundary conditions). Introducing the slow decrease of the temperature versus time, T (2) ∼ γ ′ t/ϵ 2 , and making the rescaling A = ϵA (1) to eliminate ϵ (we note that A(t) is the true amplitude of the mass profile δM (r, t)), the result writes

Ä(t) = γt + KA 2 , ( 75
)
where

γ = -γ ′ ∫ rc 0 drL (c) (r)g (c) (r)ζ(r) ∫ rc 0 drF (r)ζ(r) (76)
is found equal to γ = 120.2...γ ′ and

K = ∫ rc 0 drG(r)ζ(r) ∫ rc 0 drF (r)ζ(r) , ( 77
)
with

G(r) = 1 2 L (c) (r)g ′′(c) (r)F 2 ,r + g ′(c) (r)F ,r (F ,r 2 - 2 r F ,r ) + 1 r 2 F (r)F ,r (78)
is found to have the numerical value K = 12.32.... We have therefore established that the amplitude A(t) of the mass profile δM (r, t) satisfies Painlevé I equation. By definition the function ζ must satisfy, for any function M (2) (r), the integral relation

∫ rc 0 dr C(M (2) )(r)ζ(r) = 0. ( 79
)
Let us expand C as

C(M (2) ) = g (c) M ( 2
)
,r 2 + bM (2) ,r + cM (2) (80

) with b(r) = -2g (c) /r + M (c) /r 2 + L (c) g ′(c) and c(r) = M (c)
,r /r 2 , or in terms of the equilibrium values of the density and potential functions at the saddle-center

       g (c) (r) = 1 - 1 √ 1+ρ (c) , b(r) = -2g (c) r -h (c) ,r + ρ (c) ,r 2(1+ρ (c) ) 3/2 , c(r) = 4πρ (c) .
(81) Comparison between theory (thick black curves) and the numerics (thin colored curves) for the first order terms: (a) mass M (1) (r, t) (b) density ρ (1) (r, t), in scaled variables. The numerical curves correspond to the times t = 0.2 to 0.6 in Fig. 8.

Integrating equation (79) by parts, and using M (2) = 0 on the boundaries r = 0 and r = r c (see Appendix A), we find that ζ(r) must be a solution of the second order differential equation

(g (c) ζ) ,r 2 -(bζ) ,r + cζ = 0, ( 82
)
with the initial condition ζ(0) = 0 (the radial derivative ζ ,r (0) is a free parameter since the differential equation is of second order). At the edge of the star we do not have ζ(r c ) = 0, see below, but rather ζ ,r (r c ) = 0: the radial derivative of ζ vanishes because the second order differential equation (83) becomes a first order one (since g (c) (r c ) = 0, see equation ( 67)). This does not happen in the case studied in [START_REF] Chavanis | [END_REF] where the pressure-density relation was p = ρT , that leads to similar relations as here, but g (c) (r c ) = 1. The differential equation for the unknown function ζ(r) writes

g (c) (r)ζ ,r 2 + a 1 (r)ζ ,r + a 0 (r)ζ = 0, (83)
where the coefficients

{ a 1 (r) = 2g (c) ,r -b(r), a 0 (r) = c(r) + g (c)
,r 2 (r) -b ,r (r), (84) may be expressed in terms of the radial density using equations (67) and (81). It turns out that for r = r c we have g (c) = a 0 = 0, but a 1 (r c ) ̸ = 0, that gives the boundary relation ζ ,r (r c ) = 0.

T (t) = 1 -γ ′ t with γ ′ = 0.1.
The solution of equation ( 83) with the condition ζ(0) = 0 is shown in Fig. 6, red dashed line, where ζ ,r (r c ) = 0. Figure 8 shows the evolution of the maximum value M

(1) max (t) of the profile M (1) (r, t) with time (solid line). This quantity is proportional to the function A(t) that is the solution of Painlevé equation (75). It is compared with the numerical solution of the full Euler-Poisson equations (dots). We see that the results agree for small amplitudes but that the agreement ceases to be correct at large amplitudes where our perturbative approach loses its validity. It particular, the real amplitude increases more rapidly, and the singularity occurs sooner, than what is predicted by Painlevé equation.

Remark: According to the results of Sec. 2, and coming back to the original (but still dimensionless) variables, we find that the collapse time in the framework of Painlevé equation is t * = t * /(K γ) 1/5 with t * ≃ 3.4, i.e.

t * = 0.79... T c Ṫ 1/5 . ( 85
)
On the other hand, close to the collapse time, the amplitude of the mass profile diverges as

A(t) ∼ (6/K)(t * -t) -2 i.e. A(t) ∼ 0.487 1 (t * -t) 2 .
(86)

Discussion

This section was devoted to an explicit derivation of the "universal" Painlevé I equation for the beginning of the collapse following the slow crossing of the saddle-center bifurcation for the equilibrium problem. We have chosen to expose this detailed derivation in a simple model of equation of state and without taking into account exchange of energy in the fluid equations. Of course this makes our analysis qualitatively correct (hopefully!) but surely not quantitatively so for real supernovae, an elusive project anyway. We have shown that the Painlevé I equation represents the actual solution of the full Euler-Poisson system until the changes out of the solution at the saddle-center equilibrium are too large to maintain the validity of a perturbative approach. Our analysis explains well that the collapse of the star can be a very fast process following a very long evolution toward a saddle-center bifurcation.

As we shall explain in the next section, after the crossing of the saddle-center bifurcation, the solution of the Euler-Poisson equations have a finite time singularity at the center. We point out that this happens when the radius of the star has the order of magnitude it had at the time of the saddle-center bifurcation. Therefore the size of the core should remain orders of magnitude smaller than the star radius, as found for the Penston-Larson solution which predicts a core containing a very small portion of the total star mass. If the saddle-center bifurcation is the key of the implosion mechanism, this result should not depend on the equation of state. However the question of how massive is the self-collapsing core has received various answers. For supernovae in massive stars, starting from the hypothesis that pressure and gravity forces are of the same order during the collapse, Yahil [31] considered equations of state of the form p = Kρ Γ with adiabatic indices in the range 6/5 < Γ ≤ 4/3. He found that the ratio of the mass inside the core and the Chandrasekhar mass is almost constant, between 1.1 and unity in this range of Γ . Moreover he found that the core moves at less than the sound speed, that was considered as essential for all its parts to move in unison [2]. In the next section we show that the hypothesis that pressure and gravity forces are of the same order is not relevant to describe the collapse. Our derivation leads to a drastically different velocity field, which is supersonic in the core and subsonic outside, tending to zero at the edge of the star.

Finite time singularity of solutions of Euler-Poisson equations: pre-collapse

The perturbation analysis presented so far can deal only with perturbations of small amplitude, that is corresponding to a displacement small compared to the radius of the star. We have seen that, at least up to moderate values of the amplitude of perturbations to the equilibrium solution, the analysis derived from Painlevé equation yields correct results, not only for the exponents, but also for all the numerical prefactors. This defines somehow completely the starting point of the "explosion of the star". But there is still a long way toward the understanding of supernovae. As a next step forward, we shall look at the dynamics of the solution of the Euler-Poisson equations with radial symmetry, starting with a quasi-equilibrium numerical solution of the equations of motion. We emphasize the importance of the initial conditions for solving the dynamics, a delicate problem which could lead to various solutions as discussed and illustrated in [6] for instance.

The most noticeable feature of our numerical study is the occurrence of a singularity at the center after a finite time.

To describe the numerical results, we must invoke a singularity of the second kind, in the sense of Zel'dovich [32]. Contrary to the singularity of the first kind where the various exponents occurring in the self-similar solution are derived by a simple balance of all terms present in the equations, a singularity of the second kind has to be derived from relevant asymptotic matching, that may require to neglect some terms, as described in the present section. The occurrence of a finite time singularity in the collapse of a self-gravitating sphere has long been a topic of investigations. An early reference is the paper by Mestel [33] who found the exact trajectory of a particle during the free-fall 9 of a molecular cloud (neglecting the pressure forces), assuming spherically symmetry. The exact Mestel solution displays a self-similar solution of the pressureless Euler-Poisson system as shown later on by Penston [START_REF] Penston | [END_REF], that leads to a finite time singularity with an asymptotic density as ρ(r) ∼ r -α with α = 12/7, smaller than 2 (an important remark, as will be shown in the next subsection). Taking account of the pressure forces, another self-similar solution was found independently by Penston [START_REF] Penston | [END_REF] and Larson [5] which is usually called the Penston-Larson solution. It is characterized by α = 2. This solution was proposed to describe the gravitational collapse of an isothermal gas assuming that pressure and gravitational forces scale the same way. This corresponds to a self-similarity of the first kind (the exponent being defined simply by balancing all the terms in the original equations) by contrast to self-similarity of the second kind, or in the sense of Zel'dovich, that we are considering below. In the Penston-Larson solution, the magnitude of the velocity remains finite, something in contradiction with our numerical findings. Moreover this solution has a rather unpleasant feature, noticed by Shu [34]: it implies a finite constant inward supersonic velocity far from the center, although one would expect a solution tending to zero far from the center, as observed numerically. We present below another class of singular solution which better fits the numerical observations than the one of Penston [START_REF] Penston | [END_REF] and Larson [5]. In the numerics we start from a physically relevant situation which consists in approaching slowly the saddle-center bifurcation in a quasi-equilibrium state. As time approaches the collapse, we observe that the numerical velocity tends to infinity in the core of the singularity and decays to zero far from the center, in agreement with the theoretical solution proposed, equations (99)-(100) below with α larger than 2. The equations we start from are the Euler-Poisson equations for the mass density ρ(r, t) 9 By free-fall, we mean a situation where the collapse is due only to the gravitational attraction, i.e. in which pressure forces are neglected. This corresponds to the Euler-Poisson system ( 12)-( 14) with p = 0. and radial speed u(r, t),

ρ ,t + 1 r 2 ( r 2 ρu ) ,r = 0. (87) ρ (u ,t + uu ,r ) = -T ρ ,r - GM (r, t)ρ r 2 , ( 88
)
with

M (r, t) = 4π ∫ r 0 dr ′ r ′2 ρ(r ′ , t). ( 89
)
In the equations above, we consider the case of an isothermal equation of state, p = ρT , which amounts to considering the equation of state [START_REF] Hoflich | Cosmic Explosions in Three Dimensions: Asymmetries in Supernovae and Gamma Ray Bursts[END_REF] in the limit of large density, that is the case in the central part of the star. The temperature T has the physical dimension of a square velocity, as noticed first by Newton, and G is Newton's constant. The formal derivation of self-similar solutions for the above set of equations is fairly standard. Below we focus on the matching of the local singularity with the outside and on its behavior at r = 0. A solution blowing-up locally can do it only if its asymptotic behavior can be matched with a solution behaving smoothly outside of the core. More precisely, one expects that outside of the singular domain (in the outer part of the core) the solution continues its slow and smooth evolution during the blow-up, characterized in particular by the fact that the velocity should decrease to zero at the edge of the star meanwhile the local solution (near r = 0) evolves infinitely fast to become singular. In summary, contrary the Penston-Larson derivation which imposes the value α = 2 by balancing the terms in the equations and leads to a free parameter value R(0), our derivation starts with an unknown α value (larger than 2), but leads to a given value of R(0). In our case the unknown α value is found after expanding the solution in the vicinity of the center of the star. This yields a nonlinear eigenvalue problem of the second kind in the sense of Zel'dovich [32], as was found, for instance, in the case of the Bose-Einstein condensation [7,8] while the Penston-Larson singular solution is of the first kind (again because it is obtained by balancing all terms in the equations).

General form of self-similar solutions

The solution we are looking after is of the type for the density ρ,

ρ(r, t) = (-t) β R ( r(-t) β/α ) , ( 90
)
and for the radial velocity u,

u(r, t) = (-t) γ U ( r(-t) β/α ) , ( 91
)
where α, β and γ are real exponents to be found. The functions R(.) (different from the function R(t) introduced at the beginning of this paper. We keep this letter to remind that it is the scaled density ρ) and U (.) are numerical functions with values of order one when their argument is of order one as well. They have to satisfy coupled differential equations without small or large parameter (this also concerns the boundary conditions). To represent a solution blowing up at time t = 0 (this time 0 is not the time zero where the saddle-center bifurcation takes place; we have kept the same notation to make the mathematical expressions lighter), one expects that the density at the core diverges. This implies β negative. Moreover this divergence happens in a region of radius tending to zero at t = 0. Therefore α must be positive. Finally, at large distances of the collapsing core the solution must become independent on time. This implies that R(.) and U (.) must behave with

ξ = r(-t) β/α , ( 92
)
as power laws when ξ ≫ 1 such that the final result obtained by combining this power law behavior with the pre-factor (-t) β for R and (-t) γ for U yields functions ρ and u depending on r only, not on time. Therefore one must have

R(ξ) ∼ ξ -α , ( 93
)
and

U (ξ) ∼ ξ -γα/β . ( 94
)
In that case,

{ ρ(r, t) ∝ r -α , u(r, t) ∝ r -γα/β , (95)
for r → +∞ where the proportionality constants are independent on time.

Inserting those scaling assumptions in the dynamical equations, one finds that equation (87) imposes the relation

β α + γ + 1 = 0. ( 96
)
This relation is also the one that yields the same order of magnitude to the two terms u ,t and uu ,r on the left-hand side of equation ( 88). If one assumes, as usually done, that all terms on the right-hand side of equation ( 88) are of the same order of magnitude at t tending to zero, this imposes α = -β = 2 and γ = 0. This scaling corresponds to the Penston-Larson solution. However, let us leave α free (again contrary to what is usually done where α = 2 is selected) and consider the relative importance of the two terms in the right-hand side of equation ( 88), one for the pressure and the other for gravity. The ratio pressure to gravity is of order t 2β/α-β . Therefore the pressure becomes dominant for t tending to zero if α < 2, of the same order as gravity if α = 2 and negligible compared to gravity if α > 2 (in all cases for β negative). For pressure dominating gravity (a case where very likely there is no collapse because the growth of the density in the core yields a large centrifugal force acting against the collapse toward the center), the balance of left and right-hand sides of equation ( 88) gives γ = 0 and β = -α, while in the opposite case, i.e. for α > 2, it gives

β = -2, ( 97
) and γ = 2/α -1. ( 98
)
Therefore the velocity in the collapse region where r ∼ (-t) -β/α diverges only in the case of gravity dominating pressure (α > 2). Our numerical study shows clearly that velocity diverges in the collapse region. We believe that the early numerical work by Larson [5] does not contradict our observation that α is larger than 2: looking at his Figure 1, page 276, in log scale, one sees rather clearly that the slope of the density as a function of r in the outer part of the core is close to -2, but slightly smaller than (-2). The author himself writes that this curve "approaches the form r -2 " without stating that its slope is exactly (-2), and the difference is significant, without being very large. The slope -α = -24/11 derived below fits better the asymptotic behavior in Figure 1 of Larson [5] than the slope (-2) does (the same remarks apply to Figure 1 of Penston [START_REF] Penston | [END_REF]). Therefore we look for a solution with α > 2 for which the gravitational term dominates the pressure in equation ( 88). As shown below, the existence of a solution of the similarity equations requires that α has a well defined value, one of the roots of a second degree polynomial, and the constraint α > 2 allows us to have a velocity field decaying to zero far from the singularity region, as observed in our numerics, although α < 2 yields a velocity field growing to infinity far from the collapse region, something that forbids to match the collapse solution with an outer solution remaining smooth far from the collapse. The case α = 2 imposes a finite velocity at infinity, also something in contradiction with the numerical results.

A new self-similar solution where gravity dominates over pressure

Eigenvalue problem of the second kind

In the following, we assume that gravity dominates over pressure forces, i.e. α > 2. The set of two integro-differential equations (87) and (88) becomes a set of coupled equations for the two numerical functions R(ξ) and U (ξ) such that

ρ(r, t) = (-t) -2 R(r(-t) -2/α ), (99)
and

u(r, t) = (-t) -1+ 2 α U (r(-t) -2/α ), (100)
where ξ = r(-t) -2/α is the scaled radius. As explained previously, we must have

R(ξ) ∼ ξ -α , and U (ξ) ∼ ξ -(α/2-1) , ( 101
)
for ξ → +∞ in order to have a steady profile at large distances. The equations of conservation of mass and momentum become in scaled variables

2R + 2ξ α R ,ξ + 2 ξ RU + (RU ) ,ξ = 0, ( 102
) ( 1 - 2 α ) U + 2 α ξU ,ξ + U U ,ξ = - 4πG ξ 2 ∫ ξ 0 dξ ′ ξ ′2 R(ξ ′ ). ( 103
)
The integro-differential equation ( 103) can be transformed into a differential equation, resulting into the following second order differential equation for U (.), supposing R(.) known,

U ,ξ 2 ( U + 2 α ξ ) + U ,ξ [ 1 + 4 α + 2 ξ U + U ,ξ ] - 2γ ξ U + 4πGR = 0. ( 104
)
From now on, we use the dimensionless variables defined in Sec. 3.3. Concerning the initial conditions (namely the conditions at ξ = 0), they are derived from the possible Taylor expansion of U and R near ξ = 0, like

R = R 0 + R 2 ξ 2 + R 4 ξ 4 + ... ( 105
)
and

U = U 1 ξ + U 3 ξ 3 + U 5 ξ 5 + ... ( 106
)
Putting those expansions in equations ( 102) and ( 103), one finds U 1 = -2/3 and R 0 = 1/(6π). Note that R = R 0 and U = ξU 1 is an exact solution of the equations ( 102) and ( 103), that is not the usual case for such Taylor expansions. This corresponds to the well-known free-fall solution of a homogeneous sphere [START_REF] Penston | [END_REF]. It follows from this peculiarity that, at next order, we obtain a linear homogeneous algebraic relation because the zero value of R 2 and U 3 must be a solution. Inserting the above values of R 0 and U 1 at this order, we obtain the homogeneous relations

3 -α 3α R 2 + 5 24π U 3 = 0, ( 107
) and 4πR 2 + 5 12 -5α 3α U 3 = 0. ( 108
)
This has a non trivial solution (defined up to a global multiplying factor -see below for an explanation) if the determinant of the matrix of the coefficients is zero, namely if α is a root of the second degree polynomial

7 3 α 2 -18α + 24 = 0. ( 109
)
This shows that α cannot be left free and has to have a well defined value. However, it may happen that none of these two values of α is acceptable for the solution R(ξ), U (ξ) we are looking for, so that we should take R 2 = U 3 = 0 and pursue the expansion at next order. This is the case for our problem because one solution of equation ( 109) is α = 12/7 which does not belong to the domain α > 2 we are considering (because we assume that the gravity effects are stronger than the pressure effects) 10 , and the other solution α = 6 is excluded by the argument in section 5.2.2 below. Therefore we have to choose R 2 = U 3 = 0 and consider the next order terms of the expansion, which also provides a homogeneous linear system for the two unknown coefficients R 4 and U 5 . It is could be the relevant one for our problem. In that case, we get β = -2 and γ = -1/12. The density decreases at large distances as r -24/11 and the velocity as r -1/11 (while in the Penston-Larson solution, the density decreases at large distances as r -2 and the velocity tends to a constant value). Of course, we can carry this analysis by beginning the expansion with an arbitrary power k bigger than 2 like

4 3 -α 3α R 4 + 7 12π U 5 = 0, (110)
R = R 0 + R k ξ k + ... and U = U 1 ξ + U k ξ k+1 + ...
with arbitrary k (actually, k must be even for reasons of regularity of the solution). In that case, we find the two exponents

α(k) = 6k 2k + 3 ( 114
)
and α = 3k/(k -1). We note that the first exponent varies between 0 (homogeneous sphere) and 3, while the second exponent is larger than 3 for k > 1 which is unphysical by the argument in section 5.2.2.

In the case considered above, we note that the exponent α(4) = 24/11 is close to 2 so that it is not in contradiction with previous numerical simulations analyzed in terms of the Penston-Larson solution (which has α = 2). Moreover there is obviously a freedom in the solution because, even with α root of the secular equation, R 4 and U 5 are determined up to a multiplicative constant. This is the consequence of a property of symmetry of the equations ( 102) and ( 103

): if (R(ξ), U (ξ)) is a solution, then ( R(ξ/λ), λ -1 U (ξ/λ)
) is also a solution with λ an arbitrary positive number. This freedom translates into the fact that U 5 and R 4 are defined up to a multiplication by the same arbitrary (positive) constant. If U 5 and R 4 are that we are considering) since α = 12/7 < 2. This is why Penston [START_REF] Penston | [END_REF] and Larson [5] considered a self-similar solution of the isothermal Euler-Poisson system (87)-( 89) where both pressure and gravity terms scale the same way. Alternatively, by assuming a more general expansion ρ = ρ0 + ρ k r k + ... with k > 2 close to the origin, we find a new self-similar solution where gravity dominates over pressure.

) = R0 + R4 exp(4yi), V (yi) = U1 + U5 exp(4yi), V,y(yi) = U1 + 4U5 exp(4yi) at yi = -10, with R0 = 1 6π , U1 = -2 3 , R4 = -7(8-3α) 4πα U5 and U5 = 10 2 .
multiplied by λ, the next order coefficients of the Taylor expansion, like U 9 and R 8 (U 7 and R 6 being set to zero) should be multiplied by λ 2 , and more generally the coefficients U 4n+1 and R 4n , n integer, by λ 2n , the coefficients U 2n and R 2n+1 being all zero. The behavior of U (ξ) and R(ξ) at ξ → ∞ was derived in equation (101). As one can see, the power law behavior for R at ξ infinity follows from the assumption that terms linear with respect to R in equation ( 102) become dominant at large ξ. Keeping the terms linear with respect to U in equation ( 103) and canceling them yields U (ξ) ∼ ξ 1-α/2 . This shows that both the perturbation to u and ρ described by the self-similar solution have first a constant amplitude far from the core (defined as the range of radiuses r ∼ (-t) 2/α ) and then an amplitude tending to zero as the distance to the core increases, which justifies that the linear part of the original equation has been kept to derive this asymptotic behavior of the similarity solution. As already said, this large distance behavior of the self-similar solution makes possible the matching of this collapsing solution with an outer solution behaving smoothly with respect to time.

The numerical solution of equations ( 102)-( 103) was actually obtained by using the system (115)-(116) for the coupled variables R, V = U/ξ, then changing the variable ξ into y = ln(ξ). It writes

2R + 2 α R ,y + 3RV + (RV ) ,y = 0, (115)
and where

A ,y (V ) + 3A(V ) + 4πR(y) = 0, (116)
A(V ) = V + 2 α V ,y + V 2 + V V ,y .
The self-similar solutions R(ξ) and -U (ξ) are drawn in log scale in Figs. 9 and 10 respectively together with the corresponding time dependent density and velocity ρ(r, t) and -u(r, t). In Appendix B, by proceeding differently, we obtain the selfsimilar solution of the free-fall analytically, in parametric form. As shown later, the analytical solution is equivalent to the numerical solution of equations ( 115)-( 116), see Fig. 18.

An upper bound for α

We have seen that α must be larger than 2. It is interesting to look at a possible upper bound. Such a bound can be derived as follows. At the end of the collapse, the density and radial velocity follow simple power laws near r = 0, derived from the asymptotics of the self-similar solution. As said below, at the end of the collapse one has precisely ρ(r) ∼ r -α . Therefore, from elementary estimates, the total mass converges if α is less than 3, which gives an upper bound for α. In summary, the exponent α has to be in the range 2 < α < 3 (117) in order for a physically self-similar solution to fulfill the condition that gravity is dominant over pressure.

Homologous solution for general polytropic equations of state

The self-similar solution that we have found is independent on the pressure term in the original equation for momentum. Therefore, it is natural to ask the question of its dependence on the equation of state (namely the pressuredensity relation). Because the density diverges at r = 0 in the similarity solution, it is reasonable to expect that, if the pressure grows too much at large densities, it will become impossible to neglect the pressure term compared to gravity. Let us consider a pressure depending on ρ with a power law of the form p = Kρ Γ with Γ ≡ 1 + 1/n a real exponent and K a positive constant. We know already that, if Γ = 1, the pressure term can be neglected in the collapsing core, and the collapsing solution is characterized by the exponent α = 24/11. The same system of equations ( 102)-( 103) for the self-similar solution will be found whenever the pressure can be neglected. Therefore we expect that the above solution is valid, with the same α, as long as the power Γ in the pressure-density relation leads to negligible pressure effects in the collapsing region.

Putting the power law estimate derived from the similarity solution without pressure, one finds that the marginal exponent Γ is Γ c = 2 -2/α which for α = 24/11 is equal to

Γ c = 13 12 , (n c = 12). ( 118
)
For Γ > Γ c , the pressure becomes formally dominant compared to gravity in the collapse domain (still assuming α = 24/11), although if Γ is less than Γ c the pressure is negligible compared to gravity in the same collapse domain. When the pressure is dominant, either there is no collapse because the outward force it generates cannot physically produce an inward collapse, or other scaling laws with a different α yield a collapsing solution different from the one that we have derived (see below). If Γ is less than Γ c = 13/12 the collapse is driven by dominant gravity forces and the scaling laws derived above apply and are independent on the value of Γ . This occurs because the values of the exponents α = 24/11, β = -2, and γ = -1/12 were deduced from the Euler-Poisson equations after canceling the pressure term in the right-hand side of equation ( 88).

Let us be more general and consider other possible values of α.

If we assume that pressure and gravity forces are of the same order, the exponents are

α = 2 2 -Γ , β = -2, γ = 1 -Γ. ( 119
)
The condition α < 3 (see Section 5.2.2) implies that Γ < 4/3. It is well-known that a polytropic star with index Γ > 4/3 is dynamically stable, so there is no collapse. The critical index Γ = 4/3 corresponds to ultra-relativistic fermion stars such as white dwarfs and neutron stars. In that case, the system collapses and forms a core of mass of the order of the Chandrasekhar mass as studied by Goldreich and Weber [35]. The collapse of polytropic spheres with 6/5 ≤ Γ ≤ 4/3 described by Euler-Poisson equations has been studied by Yahil [31]. For Γ < 4/3, the star collapses in a finite time but since α < 3 the mass at r = 0 at the collapse time t = 0 is zero (in other words, the density profile is integrable at r = 0 and there is no Dirac peak). We can also consider the case where gravity forces overcome pressure forces so that the system experiences a free fall. If we compare the magnitude of the pressure and gravity terms in the Euler-Poisson system when the homologous solutions (90)-( 91) are introduced, we find that the pressure is negligible if α > 2/(2 -Γ ). Therefore, for a given polytropic index Γ , the pressureless homologous solutions are characterized by the exponents

2 2 -Γ < α ≤ 3, (120)
and

β = -2, γ = 2/α -1. (121)
The collapse exponent α is selected by considering the behavior of the solution close to the center. Setting R(ξ

) = R 0 + R k ξ k and U (ξ) = U 1 + U k+1 ξ k+1 , the relation (114)
between α and k leads to the following choice: α will be the smallest value of α(k) satisfying both relations (120) and (114) for k even. If follows that

α = 12 7 for Γ ≤ 5 6 , (122)
which is the exponent derived by Penston [START_REF] Penston | [END_REF] for zero pressure or T = 0 assuming k = 2. Next, we find

α = 24 11 , for 5 6 < Γ ≤ 13 12 , (123)
as obtained above assuming k = 4. Finally, we find that

α = 6k 2k + 3 , for 4k -3 3k < Γ ≤ 4k + 5 3k + 6 , ( 124
)
for any k ≥ 4 even. We note that there is no solution for Γ ≥ 4/3 since the polytropic stars with such indices are stable as recalled above. Finally, when pressure forces dominate gravity forces, the scaling exponents are obtained by introducing the selfsimilar form (90)-(91) into the Euler-Poisson system without gravity forces, yielding

β = - 2 2/α + Γ -1 , γ = - Γ -1 2/α + Γ -1 . ( 125
)
However, this situation is not of physical relevance to our problem since it describes a slow "evaporation" of the system instead of a collapse.

Comparison of the self-similar solution with the numerical results

Invariant profiles and scaling laws

The numerical solutions of the full Euler-Poisson system were obtained using a variant of the centpack program [29] by Balbas and Tadmor. Comparing our theoretical predictions of the self-similar solution just before collapse with the numerical solution of the full Euler-Poisson system, we find that both lead to the same result, namely they give a value of the exponent α slightly larger than two. The numerical solutions of ρ(r, t) and u(r, t) versus the radial To draw the self-similar curves, we may get around the difficult task of the exact determination of the collapse time by proceeding as follows. We define a core radius r 0 (t) such that ρ(0, t)r 0 (t) α = 1 (or any constant value), then we draw ρ(r, t)/ρ(0, t) and u(r, t)/u(r 0 , t) versus r/r 0 (t). The merging of the successive curves should be a signature of the self-similar behavior. The result is shown in Figs. 13 and14 for the density and velocity respectively. The log scale of the density curve illustrates the expected asymptotic behavior (large ξ values) R ∼ ξ -α or ρ(r, t)/ρ(0, t) ∼ (r/r 0 (t)) -α . The asymptotic behavior of the velocity, U ∼ ξ 1-α/2 is less clear on Fig. 14 where the curves display an oscillating behavior below the line with slope 1-α/2. We attribute the progressive decrease of the curves below the expected asymptote to the shock wave clearly visible in the outer part of the velocity curves (in addition, as discussed by Larson [5] p. 294, the velocity profile approaches the self-similar solution much slower than the density). In Figs. [START_REF] Emden | Gaskugeln anwendungen der mechanischen wärmetheorie auf kosmologische und meteorologie probleme[END_REF] In Fig. 13 the merging density curves have all the same ordinate at the origin, since we have plotted ρ(r, t)/ρ(0, t). To complete the comparison between the theory and the simulation for the self-similar stage, we have also drawn the series of self-similar density curves R(ξ) in order to check whether the central behavior of the numerical curves agrees with the expected value R(0) = 1/(6π). To do this we have first to define the collapse time as precisely as possible, then to plot the quantity (t * -t) 2 ρ(r, t) versus r/(t * -t) 2/α . These curves are shown in Fig. 15. They clearly merge except in a close domain around the center. We observe that the numerical value at ξ = 0 is noticeably larger than the expected value R(0) = 1/(6π) ≃ 0.05 (it is also substantially larger than the value 0.133... corresponding to the Penston-Larson solution). This shows that the system has not entered yet deep into the self-similar regime. Therefore, our numerical results should be considered with this limitation in mind. However, a precise study displays a clear decrease of the value of (t * -t) 2 ρ(0, t) during the approach to collapse, as illustrated in Fig. 16, which shows a good trend of the evolution (see below). versus t to show a quasi-linear time dependence of the numerical solution in the Painlevé and in the pre-collapse regimes. The green dots are the numerical results, the blue dotted-dashed curve is the Painlevé solution, the red curve the self-similar one which includes an additional second order term, see text.

In Fig. 16 we compare the numerics with the theory, both in the Painevé regime described in section 4 and in the self-similar regime described here. In these two regimes, the central density is expected to behave as (t * -t) -2 , see equations (11) and (99) for the Painlevé and the homologous regime respectively. Therefore we draw ρ(0, t) -1/2 which should decrease linearly with time (with different slopes). The green dots result from the numerical integration of the full Euler-Poisson equations at constant temperature (actually similar results are obtained with a temperature decreasing with time), with initial condition at temperature T = 0.9T c (out of equilibrium). At the beginning of the integration, in the Painlevé regime, the den-sity ρ c (t) = ρ(0, t) is expected to evolves as ρ c (0)+ρ 1 A(t), where A(t) is the solution of the modified version of equation (75), valid for constant temperature, which writes

Ä = ( 1 - T T c ) γ + KA 2 , ( 126
)
with γ = 120.2, and K = 12.32, as in section 4. The dotted-dashed blue line displays the function ρ p (t) -1/2 with ρ p (t) = 36A(t) + 26.85, where the coefficients are fitted to the numerical Euler-Poisson solution, and the initial conditions for the Painlevé equation are A(0) = Ȧ(0) = 0. Close to the collapse time (t * = 0.55 in the numerics), the numerical solution ρ(0, t) is expected to behave as 1 6π (t * -t) -2 , up to an additional second order term. A term of order (t * -t) -4/3 was chosen because it is the perturbation associated to the eigenvalue λ = -2/3 of the linear analysis around the fixed point C = [R 0 = 1/(6π); U 1 = -2/3] 11 and fits well the numerical results. The red curve displays the function ρ f (t) -1/2 with ρ f (t) = 1 6π (t *t) -2 + 6.5(t *t) -4/3 , which agrees well with the numerical dots, indicating that the Euler-Poisson solution tends to converge towards the self-similar form close to the center, whereas with some delay. In the following subsection we show that the fixed point C is a saddle point, with one stable direction but another unstable. It follows that the numerical solution has a priori no reason to reach C. However we observe that it clearly tends towards this fixed point as the collapse is approached.

Dynamical behavior close to the center

Recall that we have derived the theoretical value of the exponent α = 24/11 by expanding the density as R(ξ) = R 0 + R 4 ξ 4 + ... close to ξ = 0, with R 0 = 1/(6π) and U = U 1 ξ +U 5 ξ 5 +... with U 1 = -2/3 (R 4 and U 5 being defined up to a multiplicative coefficient). In order to explain the discrepancy between the numerics and the theoretical value R(0) = 1/(6π), we look at the stability of the selfsimilar solution close to ξ = 0. Let us assume here that R and U are functions of ξ and time, with ξ = r(-t) -2/α and define the time dependent variable [36]:

s = -ln (-t). ( 127
)
We set ρ(r, t) = (-t) -2 R(ξ, s), (128) and u(r, t) = (-t) -1+ 2 α U (ξ, s).

(129) where the variable s is positive for small t, increasing up to infinity as collapse is approached. Substituting this ansatz in equations (87)-(88) which include the terms due to pressure and gravity, yields the dynamical equations for R and U :

R ,s + R ,ξ ( U + 2 α ξ ) + RU ,ξ + 2 ξ R(U + ξ) = 0, ( 130
)
11 See the next subsection where the change of variable in equation ( 127) is a trick converting a problem with algebraic decay into exponential decay permitting spectral analysis. and

U ,s + U ,ξ ( U + 2 α ξ ) -γU + T R ,ξ R e 2γs + 4π ξ 2 ∫ ξ 0 dξ ′ ξ ′ 2 R(ξ ′ , s) = 0, ( 131
)
where γ is negative, see equation ( 98). These two coupled equations generalize the self-similar study of Larson [5], Penston [START_REF] Penston | [END_REF] and Brenner-Witelski [6] to the case of an exponent α different from 2. Besides the fact that in equations ( 130)-( 131) the α-dependent coefficients are slightly different from theirs, the main difference with previous works is that here the prefactor e 2γs of the pressure term decreases as s increases (as the collapse is approached), while this factor was unity in their case.

The self-similar functions U and R can be expanded as R(ξ, s) = R 0 (s) + R 2 (s)ξ 2 + R 4 (s)ξ 4 + ... and U = U 1 (s)ξ + U 3 (s)ξ 3 + ... close to ξ = 0. Writing R i (s) = R i + r i (s) and U i (s) = U i + u i (s), for i = 0, 1, 2... one gets the asymptotic relations R 0 = 1/(6π) and U 1 = -2/3 at lowest order, which is strictly the steady-state values found above in the equations without pressure, because asymptotically the pressure term vanishes. However these asymptotic values are not stable, as we shall prove now.

Because we are interested in what happens just before the collapse time, we can neglect the pressure term in equation ( 131). It becomes

U ,s + U ,ξ ( U + 2 α ξ ) -γU + 4π ξ 2 ∫ ξ 0 dξ ′ ξ ′ 2 R(ξ ′ , s) = 0.
(132) The autonomous system (130) and (132) has the useful property to reduce itself to a closed set of ODE's for R 0 (s) and U 1 (s). This set reads Writing R = R 0 + δre λs , and U = ξ(U 1 + δue λs ), the linear stability analysis of equations ( 133)-(134) in the vicinity of the fixed point [R 0 , U 1 ] gives the eigenvalues equation

U 1,s + U 1 ( U 1 + 2 α ) + ( 1 - 2 α ) U 1 + 4π 3 R 0 = 0, (133)
λ 2 + (5U 1 + 3)λ + (2U 1 + 1)(3U 1 + 2) -4πR 0 = 0. (135)
It follows that the fixed point C has one unstable and one stable direction in the phase plane, with eigenvalues +1 and -2/3, independently of the α value.

The fixed point R 0 = 0 and U 1 = -1 has two unstable directions with a degenerate eigenvalue +1, although R 0 = U 1 = 0 is stable in all directions, with eigenvalues -1 and -2. The consequences for the whole solution are not completely clear. This could explain why in the numerical work it seems so hard to get the right value of R 0 . This could be either because the initial condition for this set of ODE's does not allow to reach the fixed point U 1 = -2/3 and R 0 = 1/(6π) or because the numerics does not have the accuracy necessary to reach in logarithmic times the fixed point. Moreover, this fixed point, because it is stable in only one direction and unstable in the other, is reached from special initial conditions, on its stable manifold. Otherwise the solution are attracted either to infinity or to R 0 = U 1 = 0, depending on the initial condition.

Near the stable fixed point

Assuming that the solution approaches the stable fixed point R 0 = U 1 = 0, one may write R(s, ξ) = δr(s, ξ) and U (s, ξ) = ξδu(s, ξ), where δr and δu are smaller than unity. Setting x = -ln(ξ), the functions δr(s, x) and δu(s, x) are solutions of a linear autonomous system derived from equations ( 130) and (132). We obtain

δr ,s (s, x) - 2 α δr ,x (s, x) + 2δr(s, x) = 0, (136)
δr(s, x) = e -2s r ( 2 α s + x ) , ( 138
)
where r = δr(s, 0) is the profile of the density at the initial time t 0 of the collapse regime, with s = -ln(t 0 -t * ) by definition. It follows that the solution of the linear equation (136) decreases exponentially to zero as the collapse is approached.

Beyond the singularity: post-collapse

The question of the post-collapse was considered by Yahil [31] in his study of Euler-Poisson equations with a polytropic equation of state p = Kρ Γ with 6/5 ≤ Γ ≤ 4/3. For the critical index Γ = 4/3, corresponding to ultrarelativistic neutron stars, during the homologous collapse all the mass in the core contracts towards the center, such that at the singularity time there is a non-zero mass, of the order of the Chandrasekhar mass, at r = 0 [35]. In that case, the post-collapse regime begins with a non-zero mass at r = 0, represented in the equations by a Dirac peak at r = 0. This is not what happens for polytropic equations of state with Γ < 4/3 when pressure and gravity are of the same order [START_REF] Penston | [END_REF]5,31], or in our description of the self-similar collapse where gravity overcomes pressure forces (free fall), because, at the singularity time t = 0, as we have seen, the density does not write as a Dirac distribution but as a power law ρ(r, 0) ∝ r -α which yields for α < 3 a mass converging at r = 0 (the large distance behavior is to be matched with an outer solution to make the total mass finite). Because we do not expect a Dirac peak of finite mass at r = 0 at the time of the singularity, our post-collapse situation looks (mathematically) like the one of the dynamics of the Bose-Einstein condensation where the mass of the condensate begins to grow from zero after the time of the singularity [7,8] 12 .

Let us derive the equations for the self-similar dynamics after the collapse. As in the case of the post-collapse dynamics of self-gravitating Brownian particles [37] and of the Bose-Einstein condensation [7,8], we have to add to the equations of density and momentum conservation an equation for the mass at the center. Let M c (t) be this mass. It is such that M c (0) = 0. We need an equation for its growth. The mass flux across a sphere of radius r is J = 4πr 2 ρ(r)u(r). Therefore the equation for M c (t) is

M c,t = [ -4πr 2 ρ(r)u(r) ] r→0 . ( 139
)
To have a non zero limit of [ -4πr 2 ρ(r)u(r) ] as r tends to zero constrains the behavior of u(r) and ρ(r) near r = 0. The velocity near r = 0 is a free-fall velocity. At r very small, it is completely dominated by the attraction of the mass at r = 0. Therefore it can be estimated by taking the relation of energy conservation in free-fall, with a zero total energy, because at such short distances the initial velocity is negligible compared to the velocity of free-fall. This yields u ≈ -(2M c /r) 1/2 , which shall define the limit behavior of u(r, t) near r = 0. Because r 2 ρ(r)u(r) must tend to a finite value at r = 0, one must have ρ(r) ∼ r -3/2 . Note that this gives an infinite density at r = 0 for t > 0 while ρ(0) was finite before the collapse time; but close to r = 0 the density ρ(r) decreases (versus r) less rapidly for positive t than it did for negative t.

The equations one has to solve now are the same as before plus the attraction by the mass M c (t) at r = 0 included (the pressure being again considered as negligible, which is to be checked at the end),

ρ ,t + 1 r 2 ( r 2 ρu ) ,r = 0, ( 140
)
u ,t + uu ,r = - GM (r, t) r 2 , ( 141
)
and

M (r, t) = 4π ∫ r 0 dr ′ r ′2 ρ(r ′ , t) + M c (t). ( 142
)
The equation (139) for M c (t) with the initial condition M c (0) = 0 has to be added to the set of equations of motion. The scaling laws of this system are derived as was done for the self-similar dynamics before the singularity. Because the equations after singularity include the whole set of equations leading to the singularity, the scaling laws are the same as before, with a free exponent like the one denoted as α (this assuming, as we shall check it, that the scaling laws have as much freedom as they had before collapse, which is not necessarily true because one has another equation (139) for another unknown function, M c (t)). But the free exponent has to be the same as before collapse, because the asymptotic behavior of the solution remains the same before and after collapse: at very short times after collapse only the solution very close to r = 0 is changed by the occurrence of a finite mass at r = 0, a mass which is very small at short positive time. Therefore we look for a self-similar solution of the equations above with the same scaling laws as before collapse for ρ(r, t) and u(r, t) plus another scaling for M c (t):

ρ(r, t) = t -2 R + (rt -2/α ), ( 143
) u(r, t) = t -1+ 2 α U + (rt -2/α ), (144)
and

M c (t) = K M t b , (145)
where α = 24/11 and b is a positive exponent to be found. Moreover there has been a change of sign from (-t) to t in the scaled functions, which is obviously due to the fact that we are looking for positive times after the singularity, this one taking place at t = 0. To have the two terms on the right-hand side of equation (142) of the same order of magnitude with respect to t imposes

b = 6 α -2, (146)
a positive exponent as it should be (recall the condition that α is less than 3). For α = 24/11, we get b = 3/4. This yields the following set of definitions of the self similar unknowns after collapse,

ρ(r, t) = t -2 R + (ξ + ), ( 147
) u(r, t) = t 2/α-1 U + (ξ + ), (148)
and

M c (t) = K M t 6/α-2 . ( 149
)
The stretched radius is ξ + = rt -2/α . The equations to be satisfied by the scaled functions are

-2R + - 2ξ + α R +,ξ+ + 2 ξ + R + U + + (R + U + ) ,ξ+ = 0, (150) and ( 1 - 2 α ) U + + 2 α ξ + U +,ξ+ -U + U +,ξ+ = G ξ 2 + ( 4π
∫ ξ+ 0 dξ ′ + ξ ′ 2 + R + (ξ ′ + ) + K M ) . ( 151
)
The coefficient K M in equation ( 151) is related to the limit values of R + and U + near ξ + = 0. The solution of the two equations near ξ + = 0 are

R + ≈ K R ξ -3/2 + , (152)
and

U + ≈ K U ξ -1/2 + . ( 153
)
Equation ( 150) does not constrain the coefficients K's. By setting to zero the coefficient of the leading order term, of order ξ -5/2 + near ξ + = 0, in equation ( 151) yields a relationship between the K's,

K U = -(2GK M ) 1/2 . ( 154
)
Another relation comes from equation (139). It yields

K M = - 2π 3/α -1 K U K R . ( 155
)
Therefore there is only one free parameter among the three coefficients K's. This free parameter is fixed by the matching with the large distance behavior of R + and U + , which is defined itself by the matching with the outside of the collapse domain. The system (150)-( 151) was solved numerically by using the coupled variables R, V = U/ξ and y = ln(ξ) (dropping the + indices) as in subsection 5.2, that gives the coupled equations analogous to equations (115)-(116),

{ -2R -2 α R ,y + 3RV + (RV ) ,y = 0 A ,y (V ) + 3A(V ) -4πR(y) = 0 (156) with A(V ) = V + 2 α V ,y -V 2 -V V ,y
, which are free of the inner core mass term. In Appendix B, by proceeding differently, we obtain an analytical solution of the postcollapse dynamics which agrees with the numerical solution of equations (156), see Fig. 17 for comparison.

Conclusion and perspectives

This contribution introduced a theory of the early stage of supernova explosion which assumes that this belongs to the wide class of saddle-center bifurcations with a slow sweeping of the parameters across the bifurcation range. This explains well the suddenness of the explosion occurring after aeons of slow evolution. The hugely different time scales combine into a single intermediate time scale for the slow-to-fast transition which could be of the order of several hours. This transition is described by a "universal" dynamical equation, the Painlevé I equation. Comparing this prediction with a model of star presenting a saddle-center bifurcation shows a quantitative agreement with the predictions based on general arguments of bifurcation theory.

This shows at least one thing, namely that the collapse of the star by the loss of equilibrium between pressure and gravitational forces is a global phenomenon depending on the full structure of the star and cannot be ascribed, for instance, to an instability of the core reaching the Landau-Chandrasekhar limit mass, as often assumed. We also looked at the evolution of the star following the onset of instability, namely when the amplitude of the perturbations grows to finite values and cannot be described by the Painlevé I equation anymore. In our equation of state model, the pressure becomes proportional to the density in the large density limit. The pressure increase is likely less steep than what is expected for the inner core of stars, even though there are big uncertainties on the interior of stars, particularly the ones yielding supernovae: showing no early warning on the incoming explosion they are not scrutinized spectroscopically. Nevertheless, an analysis of this situation teaches us a few interesting lessons. First, we do not consider self-similar (or homologous) collapse in the usual sense (where pressure and gravity scale the same way) because our numerical results and our analysis lead us to claim that the pressure becomes negligible in the core. Secondly, we find a new self-similar free-fall solution toward the center.

Our numerical results together with physical considerations about the velocity field make us argue that besides the mathematically correct Penston-Larson solution, our new self-similar (free-fall) solution is relevant to describe the collapse. In other words, writing self-similar equations is not enough to guaranty their relevance for a given problem because there can be more than one such kind of solution, like in the present case, where Zel'dovich type 2 solution corresponds to the numerical results, although a type 1 solution also exists, but is not relevant.

The numerical results presented here were obtained by starting from the equilibrium state of the star at the saddle-center, then decreasing slowly the temperature. However we notice that the same conclusions are obtained when starting slightly away from the saddle-center point and performing the numerical integration at constant temperature. We point out that the previous numerical studies of gravitational collapse by Penston [START_REF] Penston | [END_REF], Larson [5] and later by others [6] were performed starting from a uniform density initial state (and finite radius), that represents initial conditions which are very far from ours and from any physical situation; nevertheless these authors did find a density behaving asymptotically (at large distance) as r -α , with α larger than 2, as we find here.

The free-fall solution we found is not the free-fall solution studied for many years, because the exponents of our self-similar solution are not the ones usually found. This conclusion is based upon a detailed comparison between the direct numerical solution of the evolution equations and the solution of the simpler equations for the self-similar problem. As far as we are aware, although the self-similar paradigm is often invoked in this field, such a detailed comparison between dynamical solutions of the full Euler-Poisson system and the full self-similar solution has not yet been done (the merging of the curves before the collapse time was not shown). We show that it is a relatively non trivial endeavor to perform such a comparison. Moreover we point out that our self-similar pressure-free solution is more tricky to derive than the standard Penston-Larson homologous solution including the pressure for which standard scaling laws (Zel'dovich first kind) can be derived formally without any difficulty. Finally we have mentioned that the center is a saddle point for our self-similar solution. Numerically this property is manifested in the behavior towards r = 0 of the density profile ρ(r, t) -ρ(0, t) which should pass from r 2 to r 4 in the self-similar regime (for generic initial conditions). The mechanism of this change of exponent, if it really occurs, has not been clearly identified and requires a deeper study.

This work leaves open many questions. One central issue is how the scenario we outlined, namely slow starting in the universality class Painlevé I, and later finite time collapse toward the central core, is dependent on the pressure/density relation. We suspect that, if the pressure increases more rapidly with the density than linearly at large densities, there will be no finite time singularity. Likely, because shock waves will form, irreversible transformations will take place in those shock waves and another equation of state will become relevant for the star.

We greatly acknowledge the "Fondation des Treilles" where this work was initiated, and Paul Clavin for many very stimulating discussions.

A Boundary conditions to derive the normal form

Let us derive the boundary conditions to solve the integral equation (79) by transforming it into the differential equation (82). We have to cancel the terms [ g (c) ζM (2) ,r

] rc 0 , [ (g (c) ζ) ,r M (2)
] rc 0 , and [ bζM (2) ] rc 0 .

(i) At r c we have g (c) (r c ) = 0 and M (2) (r c ) = 0 that ensure the cancelation of the terms g (c) ζM (2) ,r , g (c) ζ ,r M (2) , bζM (2) , and g (c)

,r ζM (2) at r = r c (while g (c)

,r and ζ are both non zero at r = r c , see Fig. 6). This suppresses all the terms taken at r = r c .

(ii) At r = 0 we impose ζ = 0 that cancels the terms g (c) ζM (2) ,r , g (c)

,r ζM (2) , and bζM (2) . The last term g (c) ζ ,r M (2) vanishes under the condition M (2) (0) = 0 (because g (c) and ζ ,r are both non zero at r = 0). This suppresses all the terms taken at r = 0.

B Analytical self-similar solutions for the free-fall

Penston [START_REF] Penston | [END_REF] has given an exact solution of the free-fall problem without thermodynamic pressure (p = 0). It could seem that, because of the absence of thermodynamic pressure, this is irrelevant for the problem of singularity in the evolution of the collapsing core of models of stars. However, this is not quite true because we have shown that during the collapse this thermodynamic pressure becomes negligible, and so the evolution of the system is essentially like a free-fall. By analyzing the equations for this pressureless collapse we have shown that, actually, a discrete set of solutions exists, with different singularity exponents. The free-fall solution found by Penston corresponds to the exponent α = 12/7. Since this exponent is smaller than 2 pressure effects become important at a certain point of the evolution (Penston obtains the estimate δt/t f ∼ 10 -4 ) and this is why he considers in a second step the case where pressure and gravity forces are of the same order leading to another self-similar solution (the Penston-Larson solution) with α = 2. Actually, we propose another possibility which is in agreement with our numerical results (and actually with many others). We show below that other exponents than 12/7 are possible for the free-fall, some of them being larger than 2 and providing therefore a possible solution of the initial problem in which gravity always dominates over pressure forces 13 . Our solutions are based on the choice of initial conditions for the radial dependence of the density ρ(a) = ρ 0 (1 -a k /A k ) where a is the radial variable (same notations as in Penston [START_REF] Penston | [END_REF]). The exponent k is left free, although Penston takes k = 2 with the comment: "we are 'almost always' correct in taking the form ρ(a) = ρ 0 (1 -a 2 /A 2 )".

We consider a sphere of gas initially at rest and call M (a, 0) the mass of gas contained within the sphere of 13 It does not mean that the Penston-Larson solution is incorrect. It represents a mathematically exact (type I) self-similar solution of the isothermal Euler-Poisson equations. However, we argue that other (type II) self-similar solutions exist in which gravity overcomes pressure. They are characterized by α > 2 and by a density behaving as ρ0 + ρ k r k with k > 3 close to the origin (see below), while the Penston-Larson solution has α = 2 and the density behaves as ρ0 + ρ2r 2 close to the origin. Our numerical work (despite its limitations because we follow the collapse only over a few decades in density) together with important physical considerations (e.g. the fact that the velocity profile in our solution decreases to zero instead of tending to a constant value) suggest that these new solutions are relevant to describe the collapse. radius a and ρ(a) = 3M (a, 0)/4πa 3 the average density of that sphere. Using Gauss theorem, the Euler equation [START_REF] Emden | Gaskugeln anwendungen der mechanischen wärmetheorie auf kosmologische und meteorologie probleme[END_REF] with the pressure neglected is equivalent to

d 2 r dt 2 = du dt = - GM (a, 0) r 2 , ( 157
)
where r and u are the position and the velocity at time t of a fluid particle located at r = a at t = 0. This equation can be solved analytically [33] and the solution can be expressed in parametric form as r = a cos 2 θ, (158)

t = ( 3 8πGρ(a) ) 1/2 ( θ + 1 2 sin(2θ) ) , ( 159
)
where θ runs between 0 and π/2. Taking θ = π/2, we find that a particle initially at r = a arrives at r = 0 at a time t(a) = (3π/32Gρ(a)) 1/2 . Setting a = 0 + in the foregoing expression, we find that the first particle reaches the center at the time

t f = ( 3π 32Gρ 0 ) 1/2 , ( 160
)
where ρ 0 = ρ(0). This is called the free-fall time. At t = t f , the central density becomes infinite (ρ c = +∞).

Using the equation of motion ( 158)-( 159) giving r = r(a) and the relation ρ(r, t)r 2 dr = ρ(a, 0)a 2 da, which is equivalent to the equation of continuity [START_REF] Landau | Statistical physics[END_REF], we can determine the evolution of the density profile ρ(r, t) and of the velocity profile u(r, t) in the pre-and post-collapse regimes. For t → t f and r not too large, they have a selfsimilar form. The derivation of this self-similar solution follows rather closely the one by Penston with the only difference that his assumption ρ(a) = ρ 0 (1 -a 2 /A 2 ) is replaced by ρ(a) = ρ 0 (1 -a k /A k ). Therefore, we skip the details of the derivation and directly give the final results.

B.1 The pre-collapse regime

In the pre-collapse regime (t < t f ), the self-similar density and velocity profiles are given in parametric form by ρ(r, t) ρ c (t) = where y = 1 2 ( a A ) k t f δt goes from 0 to +∞ (here δt = t ft). For k = 4, the curves ρ(r, t)/ρ c (t) and -u(r, t)/u 0 (t) drawn in Fig. 18, solid lines, coincide with the self-similar numerical solution (dashed line) of equations ( 102)-(103) derived in section 5.2.

In the above parametric representation the central density is given by the relation

ρ c (t) = ( 4 3π
) 2 ρ 0 ( t f t f -t ) 2 . ( 164
)
Using equation (160) it can be rewritten as

ρ c (t) = 1 6πG 1 (t f -t) 2 , ( 165
)
which agrees with the result of Sec. 5.2. Moreover, we have

r 0 (t) = ( 3π 4
) 2/3 2 1/k A t f -t t f (2k+3)/3k , ( 166
)
u 0 (t) = π 2 (k-1)/k ( 4 3π
) 1/3 A t f t f -t t f (3-k)/3k
. (167)

For r → 0, we get ρ(r, t) ∼ ρ c (t)

α k = 6k 2k + 3 , ν k = 3 -k 2k + 3 . ( 172
)
The expressions (170) and (171) also give the density and velocity profiles for all r at t = t f . For k = 2, we get α 2 = 12/7 and ν 2 = 1/7; for k → +∞, we get α ∞ = 3 and ν ∞ = -1/2; for k = 4, we get α 4 = 24/11 and ν 4 = -1/11. The exponent α achieves the critical value 2 for k = 3. For k < 3, i.e. α < 2, the pressure wins over gravity as we approach the collapse time t f , and the free-fall solution is not valid anymore. For k > 3, i.e. α > 2, the gravity always wins over pressure so the free-fall solution may be valid for all times.

Let us discuss the form of the density and velocity profiles depending on k.

For any k, the density profile ρ(r, t) starts from a finite value (for t < t f ) and decreases with the distance r. The central density ρ c (t) increases with time and diverges at the collapse time t f . At t = t f , the density profile is singular at the origin.

For k < 3, i.e. α < 2, the velocity profile -u(r, t) starts from zero at r = 0 and increases with the distance r. The magnitude of the velocity u 0 (t) decreases with time and tends to zero at the collapse time t f . At t = t f , the velocity is still zero at the origin.

For k = 3, i.e. α = 2, the velocity profile -u(r, t) starts from zero at r = 0 (for t < t f ), increases with the distance r, and reaches an asymptotic value u 0 (the prefactor u 0 (t) is constant). At t = t f , the velocity profile has a constant non-zero value u 0 .

For k > 3, i.e. α > 2, the velocity profile -u(r, t) starts from zero at r = 0, increases with the distance r, reaches a maximum, and decreases towards zero at large distances. The prefactor u 0 (t) increases with time and diverges at the collapse time t f . At t = t f , the velocity profile is singular at the origin.

B.2 The post-collapse regime

In the post-collapse regime (t > t f ), the self-similar density and velocity profiles are given in parametric form by ρ(r, t) ρ c (t) = where y goes from -∞ to -1, and ρ c (t), r 0 (t) and u 0 (t) are defined by equations ( 164)-(167) as in the pre-collapse regime. For r → +∞, the behavior is the same as in the pre-collapse regime, but for t > t f and r → 0, we get ρ(r, t) ∼ ρ c (t) 3 2k

( r 0 (t) r

) 3/2 , ( 176
)
u(r, t) ∼ -u 0 (t) ( r 0 (t) r

) 1/2 . ( 177
)
We note that the density and the velocity are always singular at r = 0. For any k, the density profile ρ(r, t) is decreasing, as illustrated in Fig. 19-(a) . For k < 3, the velocity profile -u(r, t) decreases, reaches a minimum value, and increases. For k = 3 it decreases towards an asymptotic value u 0 and for k > 3 it decreases towards zero, see Finally, the mass contained in the Dirac peak ρ D (r, t) = M D (t)δ(r) at time t > t f is

M D (t) = 8π 3 2 (3-k)/k ρ 0 A 3 ( t -t f t f ) 3/k . ( 178
)
The mass in the core grows algebraically with an exponent b k = 3/k. For k = 2, we get b 2 = 3/2; for k → +∞, we get b ∞ = 0; for k = 3, we get b 3 = 1; for k = 4, we get b 4 = 3/4.

30 VFig. 1 .

 301 Fig. 1. Potential evolution close to a saddle-center, equation (2) with b = c = 1 and two values of a = -ct; t = -2 for the blue curve, t = 2 for the red dashed curve.

10 RFig. 2 .

 102 Fig. 2. Numerical solution of equation (5), or equation (3) with b = c = 1, for two different initial conditions taken at time ti = -20; (i) relation (8) for the blue curve without any oscillation; (ii) R(ti) = √ -ti + 0.5 and R ′ (ti) = -1 2 √ -t i for the red oscillating curve.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Density ρ(r) versus the radial variable at the saddlecenter (T = Tc, or ĥ0 = 2.296). The density vanishes at the edge of the star indicated by the arrow (r = r0).

Fig. 5 . 3 0

 53 Fig. 5. Radius r0 = r0 M -1/3 0

  Fig. 7.Comparison between theory (thick black curves) and the numerics (thin colored curves) for the first order terms: (a) mass M(1) (r, t) (b) density ρ(1) (r, t), in scaled variables. The numerical curves correspond to the times t = 0.2 to 0.6 in Fig.8.

Fig. 8 .

 8 Fig.8. Comparison between normal form (solid blue curve) and numerical solution (red points) for the maximum of M(1) (r, t) versus time. In the numerical simulations of the Euler-Poisson system we start from the critical profile Mc(r) at t = 0 and decrease the temperature as T (t) = 1 -γ ′ t with γ ′ = 0.1.

Fig. 9 .

 9 Fig. 9. Density of the self-similar problem obtained by solving equations (115)-(116) with α = 24/11. (a) R(ξ); (b) ρ(r, t) versus r at times 1, 0.1, 0.05, 0.01, 0.001. The initial conditions are R(yi) = R0 + R4 exp(4yi), V (yi) = U1 + U5 exp(4yi), V,y(yi) = U1 + 4U5 exp(4yi) at yi = -10, with R0 = 1 6π , U1 = -23 , R4 = -7(8-3α) 4πα U5 and U5 = 10 2 .

uFig. 10 .

 10 Fig. 10. Velocity of the self-similar problem, obtained by solving equations (115)-(116) with α = 24/11. (a) -U (ξ), (b)-u(r, t) versus r at same times and with same initial conditions as in Fig.9.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Density ρ(r, t) versus the radial variable r in log 10 scale: numerical solutions of the full Euler-Poisson system, equations (87)-(88) at different times before the collapse. The solid line with slope -24/11 fits better the asymptotic behavior (large r) of the curves than the dotted-dashed line with slope -2.

Fig. 13 .Fig. 14 .

 1314 Fig.13. Self-similar density curves ρ(r, t)/ρ(0, t) versus r/r0(t) in log scale with r0(t) defined in the text and α = 24/11.

4 RFig. 15 . 20 ρ c 1 2 Fig. 16 .

 41520216 Fig. 15. Numerical self-similar density curves ρ(r, t)(-t) 2 versus ξ = r(-t) -2/α for α = 24/11.

and R 0

 0 ,s + 3R 0 U 1 + 2R 0 = 0. (134) This system has three fixed points (namely solutions independent on s): (i) the point C=[R 0 = 1/(6π); U 1 = -2/3] defined in the previous subsection (the values at ξ = 0 of R and U , solution of the similarity equations already derived); (ii) also [R 0 = U 1 = 0]; (iii) and finally [R 0 = 0; U 1 = -1].

Fig. 17 .

 17 Fig.17. Self-similar density R+(ξ+) and velocity U+(ξ+) in the post-collapse regime, solution of equations (156), in ln scale, to be compared with solutions in the pre-collapse regime drawn in Figs. 9-10 also in ln scale.

  y)1/3 , (163)

Fig. 18 .

 18 Fig. 18. Parametric solutions (161)-(163) compared with the self-similar solutions of section 5.2 for α = 24/11. (a) Density ρ(r, t)/ρc(t) versus r/r0(t) for k = 4 in solid line, R(ξ)/R(0) versus ξ = r/r0(t) for R4 = -4 in dashed line. (b) Velocity u(r, t)/u0(t) versus r/r0(t) in solid line; -1.6U (ξ) in dashed line.

) 6 /( 8πρ 0 G 3 ) 1

 631 r → +∞, we get ρ(r) ∼ ρ 0 3 (2k+3) 1 r 6k/(3+2k) , (170) u(r) ∼on time as it should. We have ρ ∼ r -α k and u ∼ r ν k with

3 3 + 2 ( 3 +

 323 k)y + (3 + 2k)y 2 , (173) r r 0 (t) = |y| 1/k |1 + y| 2/3 , (174) u(r, t) u 0 (t) = -|y| 1/k |1 + y| 1/3 , (175)

Fig. 19 -

 19 (b).

Fig. 19 .

 19 Fig. 19. Parametric solutions of equations (173)-(175) and self-similar solutions of equations (156) in the post-collapse regime. (a) density ρ(r, t)/ρc(t) versus r/r0(t) for k = 4, or α = 24/11, in solid line, 15R(ξ) versus ξ = r/r0(t) for KU = -1 in dashed line, (b) velocity -u(r, t)/u0(t) versus r/r0(t) in solid line, -U (ξ) in dashed line.

This equation of state is inspired by the study of selfgravitating boson stars in general relativity[22][23][24]. Such an equation of state could hold in the core of neutron stars because of its superfluid properties[24]. The neutrons (fermions) could form Cooper pairs and behave as bosons. In this context ρc 2 represents the energy density and the parameter T has an interpretation different from the temperature (in the core of neutron stars T is much less than the Fermi temperature or than the Bose-Einstein condensation temperature so it can be taken as T = 0). We use here this equation of state with a different interpretation.

This temperature-radius relation T (R) is the counterpart of the mass-radius relation M (R) of boson stars in general relativity, that also presents a spiralling behavior[24]. The dynamical stability of the configurations may be determined from the theory of Poincaré on the linear series of equilibria as explained in[START_REF] Chavanis | [END_REF]. If we plot the temperature T as a function of the parameter ĥ0, a change of stability can occur only at a turning point of temperature. Since the system is stable at high temperatures (or low ĥ0) because it is equivalent to a

The authors of[START_REF] Chavanis | [END_REF] study the dynamics of Smoluchowski-Poisson equations close to a saddle-node but for a fixed value of the temperature T → T - c .

We note that the exponent α = 12/7 was previously found by Penston[START_REF] Penston | [END_REF] for the free-fall of a pressureless gas (T = 0) by assuming a regular Taylor expansion ρ = ρ0 + ρ2r 2 + ... close to the origin. This solution is valid if T is exactly zero but, when T > 0, as it is in reality, this solution cannot describe a situation where gravity dominates over pressure (the situation

Some analogies between the post-collapse dynamics of selfgravitating Brownian particles[37] and the Bose-Einstein condensation have been discussed in[8].

B.3 The homogeneous sphere

Finally, for completeness, we recall the solution corresponding to the collapse of a homogeneous sphere with mass M , initial density ρ 0 and initial radius R 0 . Since ρ(a) = ρ 0 , we find from equations ( 158)-( 160) that all the particles collapse at r = 0 at the same time t f . Therefore, a Dirac peak ρ D (r) = M δ(r) is formed at t = t f . The evolution of the radius R(t) of the homogeneous sphere is given by

where θ runs between 0 and π/2. For t → t f , we get

The density ρ c (t) = 3M/4πR(t) 3 increases as

The velocity field is u(r, t) = -H(t)r with