A&A 418, 283-294 (2004)
DOI: 10.1051/0004-6361:20035619

Evolutionary sequences of rotating protoneutron stars

L. Villain ^1,2,3 - J. A. Pons ^1,4 - P. Cerdá-Durán ¹ - E. Gourgoulhon ³

1 - Departament d'Astronomia i Astrofísica, Universitat de València, 46100 Burjassot, Spain,
2 - Copernicus Astronomical Center (CAMK), Polish Academy of Sciences, Bartycka 18, 00-716 Warszawa, Poland
3 - Laboratoire de l'Univers et de ses Théories, UMR 8102 du CNRS, Observatoire de Paris - Section de Meudon, 92195 Meudon Cedex, France
4 - Departament de Física Aplicada, Universitat d'Alacant, Apartat de correus 99, 03080 Alacant, Spain

Received 3 November 2003 / Accepted 19 January 2004

Abstract
We investigate the evolution of rigidly and differentially rotating protoneutron stars during the first twenty seconds of their life. We solve the equations describing stationary axisymmetric configurations in general relativity coupled to a finite temperature, relativistic equation of state, to obtain a sequence of quasi-equilibrium configurations describing the evolution of newly born neutron stars. The initial rotation profiles have been taken to mimic the situation found immediately following the gravitational collapse of rotating stellar cores. By analyzing the output of several models, we estimate that the scale of variation of the angular velocity in a newly born neutron star is of the order of 7-10 km. We obtain the maximum rotation frequency that can be reached as the protoneutron stars deleptonizes and cools down, as well as other relevant parameters such as total angular momentum or the instability parameter |T/W|. Our study shows that imposing physical constraints (conservation of baryonic mass and angular momentum) and choosing reasonable thermodynamical profiles as the star evolves gives results consistent with the energetics of more complex simulations of non-rotating protoneutron stars. It appears to be unlikely that newly born protoneutron stars formed in nearly axisymmetric core collapses reach the critical angular velocity to undergo the bar mode instability. They could, however, undergo secular or low |T/W| rotational instabilities a few seconds after birth, resulting in a strong emission of gravitational waves retarded with respect to the neutrino luminosity peak. We also found that the geometry of strongly differentially rotating protoneutron stars can become toroidal-like for large values of the angular velocity, before reaching the mass shedding limit.

Key words: stars: neutron - stars: rotation - stars: evolution

1 Introduction

The first minute of life of non-rotating protoneutron stars (PNS) has been studied in detail since the mid 80's by several authors (Burrows & Lattimer 1986; Keil & Janka 1986; Pons et al. 1999,2000,2001). The initial hot, lepton rich remnant, left behind following a successful core collapse supernova is known to evolve in a timescale of tens of seconds to form a cold (T<10¹⁰ K) catalyzed neutron star. This process is usually classified in three main stages: i) the mantle contraction during which fast cooling of the outer regions takes place in about 0.5 s, with probably significant accretion; ii) deleptonization and consequently heating of the internal core as energetic neutrinos diffuse out leaving most of their chemical energy (analogously to Joule effect) on the way out; and iii) cooling by means of diffusion of (mostly) thermal neutrinos, resulting in a decrease of temperature from about 40-50 MeV to below 2-4 MeV, point at which the star becomes transparent to neutrinos.

Less is known about the early evolution of PNSs that are rapidly rotating, because a fully consistent study is a formidable task. It requires solving the neutrino transport equations in at least two-dimensions, coupled to general relativistic hydrodynamics that describes the fluid motion in rotating relativistic stars. In addition, the microphysical inputs also need to be treated carefully. As has been shown previously (Pons et al. 1999), internal consistency between the equation of state of dense matter and the neutrino opacities describing the interaction between neutrinos and matter is needed to achieve reliable results.

A few attempts to treat the problem in a simplified way exist (Romero et al. 1992; Goussard et al. 1997,1998; Sumiyoshi et al. 1999; Strobel et al. 1999; Yuan & Heyl 2003) based on the analysis of temporal sequences of quasi-equilibrium models of PNS, imposing some ad-hoc thermodynamical profiles (i.e., constant temperature, constant entropy and/or constant neutrino fraction). The assumption of quasi-equilibrium is well justified, since the hydrodynamical timescale (10^-3 s) is much smaller than the timescale in which substantial thermodynamical changes occur (diffusion timescale $\approx$1 s). However, it is unclear whether or not an isothermal or isentropic profile can effectively mimic reality. In addition, there are some physical constraints such as conservation of energy (the total gravitational mass must decrease consistently with the neutrino luminosity), or conservation of angular momentum, that need to be satisfied when one studies a temporal sequence of quasi-equilibrium models of PNS. These constraints help to reduce the parameter space as one calculates the thermodynamical evolution, although there is much uncertainty about the transport of angular momentum. In general angular momentum losses by neutrinos are supposed to be small (Kazana 1977). Neutrinos could redistribute or take away a fraction of the initial angular momentum, but quantifying this amount needs of multidimensional transport simulations. It is equally controversial how important is turbulent transport, or the role of gravitational radiation. Because of these many unresolved issues it is usually assumed that the angular momentum is approximately conserved during the Kelvin-Helmholtz stage. In addition, the presence of large magnetic field could be crucial. It has been recently shown that turbulent mean-field dynamo action can be effective for PNSs with periods shorter than 1 s, which would generate very strong magnetic fields in the interior (Bonanno et al. 2003).

The whole problem is far to be solved, and will require some serious computational effort but, in the meantime, the intention of this paper is to begin to understand qualitatively what kind of outcome will be revealed by such numerical simulations, prior to engaging large scale simulations. In this line, we improve in several respects the few existing previous works. First, instead of isentropic or isothermal models, we use realistic profiles coming from 1D simulations, conveniently rescaling the temperature and chemical profiles as functions of density; second, we do not restrict ourselves to the rigid rotation case, as most of previous studies, and include differential rotation in our analysis; last but not least, we check that fixing the baryonic mass and the angular momentum results in some reasonable time dependence to the integrated relevant physical quantities such as the gravitational mass. In this way, we can calculate a complete evolutionary sequence of a PNS.

The plan of the paper is as follows. In Sect. 2 we briefly review the relativistic description of rotating stars, and the numerical method used to obtain equilibrium configurations. A similar brief description of the thermodynamics and the equation of state is presented in Sect. 3, since for both issues the details have been published elsewhere. Section 4 discusses the choice of parameters for the initial models. Our results are presented in Sect. 5. In the first part we describe rigidly rotating PNSs, while the description of differentially rotating PNSs is given in Sect. 5.2. In both cases, we give all the relevant parameters for maximally rotating configurations as well as the evolution of more realistic sequences. Finally, our main conclusions and findings are summarized and discussed in Sect. 6.

   
2 Numerical code for stationary axisymmetric spacetime

Since the time scale for the matter and the gravitational field to settle in stationary equilibrium is much shorter than all other involved times, a simplifying but plausible approximation consists of mimicking the time evolution of the spacetime by obtaining a collection of stationary, axisymmetric and asymptotically flat snapshots, all of them being solutions of the Einstein equations, but with a source term that depends on the snapshot.

This approximation has already been adopted by several authors in the past (Goussard et al. 1997,1998; Strobel et al. 1999; Yuan & Heyl 2003), our code being moreover very similar to the one of Goussard et al. (1997,1998), since we also employ a fully relativistic spectral code based on the work of Bonazzola et al. (1993). In the following, this article will be referred to as BGSM, and we address the interested reader to it for more details about the algorithms. Let us just remind that the BGSM approach uses the 3+1 formulation of general relativity, maximal slicing-quasi-isotropic coordinates, with the assumption that matter is free of convective motions. The latter hypothesis is probably too strong, since PNSs are known to be convectively unstable (Keil et al. 1996; Miralles et al. 2002), but it implies the circularity of the stationary axisymmetric spacetime (Carter 1969) and greatly simplifies the calculations.

Under these approximations, the metric is written

 

$${{\rm d}} s^2 ~ \widehat{=} ~ g_{\mu\nu}~ {{\rm d}} x^{\mu} {{\rm d}} x^{\nu}~$$

$$\widehat{=} ~ - \left (N^2~-~N_{\varphi}~{N^{\varphi}}\right) {{\rm d}} t^2 - 2~N_{\varphi}~ {{\rm d}} t ~{{\rm d}\varphi}~$$

$$+~\frac{A^4}{B^2}\left ({{\rm d}} r^2~+~r^2~{{\rm d}\th}^2~+~B^4~r^2~ \sin^2 \th~{{\rm d}\varphi}^2\right),$$ (1)

where we have introduced the usual notation of the 3+1 formalism: N is the lapse, $N^{\varphi}$ is the third component of the shift 3-vector, with $N_{\varphi}~\widehat{=}~g_{\varphi i}~N^i~=~A^4~B^2~ N^{\varphi}r^2~\sin^2 \th$ being the covariant $\varphi$-component of the latter. Thanks to the hypothesis of axisymmetry and stationarity, the functions  $N, N^{\varphi}, A$ and B depend only on r and $\th$.

Despite the general formalism and structure of the code is very similar to that of Goussard et al. (1997,1998), there are a few differences that should be mentioned:

-

we use a more recent version of the spectral BGSM code (see Gourgoulhon et al. 1999), in which at least two shells map the interior of the star, in order to better describe the non-spherical surface. This code is based on the C++ library LORENE, a software package for numerical relativity freely available under GNU license;

-

the equation of state (EOS) is written in tabular form and interpolated using cubic Hermite polynomials, as defined in Nozawa et al. (1998), according to a general technique presented in Swesty (1996) that preserves thermodynamical consistency (see also further remarks below);

-

during the evolution of PNSs we do not restrict ourselves to 3 or 4 typical evolutionary snapshots that are just isentropic or isothermal profiles. Instead we adopt more realistic profiles, rescaled from the results of 1D simulations with neutrino transport (Pons et al. 1999, more details are given in the next section).

For clarity, let us be more specific about the last remark. From the hydrodynamical point of view, one must deal with the conditions for the existence of stationary solutions of the Einstein equations in presence of a rotating perfect fluid (see also BGSM and Goussard et al.). The main issue is that the relativistic Euler equations, derived from the conservation of the energy-momentum tensor projected into a space-like 3-surface, are reduced in the case of stationary rotating motion to

$$\frac{\partial_i~P}{e~+~P}~+~\partial_i~\ln \left[\frac{N}{\Gamma}\right]~ =~ -~F~\partial_i~\Omega~,$$ (2)

where:

-

P is the pressure and e the total energy density, both measured in the comoving frame;

-

N is the already defined lapse function, $\Gamma$ the Lorentz factor linking the Eulerian observer and the comoving observer, $\Omega~\widehat{=}~{u^\varphi}/{u^t}$ the angular velocity, and $F~\widehat{=}~u_\varphi~u^t$with u being the 4-velocity of the fluid;

-

the index i is either 1 or 2, since all the previous quantities depend only on r and $\th$.

As explained in BGSM Eq. (2) does not admit a first integral for arbitrary thermodynamics or rotation profiles  $\Omega[r,\th]$, but it does in some simple cases.

Concerning rotation, the right hand side of Eq. (2) is a first integral if either $\Omega$ is constant in space (rigid rotation), or the function F can be locally written as $F~=~F[\Omega]$. The latter leads to the relation

$$F[\Omega]~-~\frac{A^4~B^2~r^2~\sin^{2}\th~(\Omega~-~N^\varphi)} {N^2~-~A^4~B^2~r^2~\sin^{2}\th~(\Omega~-~N^\varphi)^2}~=~0~,$$ (3)

between $\Omega$ and the coefficients of the metric, which allows the determination of the profile of angular velocity.

For the differential rotation, we adopt the law proposed by Komatsu et al. (1989) and used by many authors (Komatsu et al. 1989; Goussard et al. 1997,1998; Baumgarte et al. 2000)

$$F[\Omega]~=~R_{0}^2~(\Omega_{\rm c}~-~\Omega)~,$$ (4)

where $\Omega _{\rm c}$ is the limit of the function $\Omega[r,\th]$ on the rotation axis and R₀² is a parameter with the dimension of the square of a length, noted A in Komatsu et al. (1989). We shall discuss later our choices for the values of these parameters. It must be remarked that despite its simplicity, this law gives reasonable rotation profiles qualitatively similar to those obtained from core collapse simulations (see Sect. 4).

Regarding the EOS, it can be shown that isentropic or isothermal profiles are sufficient conditions for Eq. (2) to be integrable. Yet another possibility which is less restrictive exists: assuming that the total energy density ecan be written as an effective function of the pressure e = e[P]. Indeed, in this case, even if the actual EOS is not barotropic, assuming that the EOS is a one-parameter function ensures the existence of the first integral

$$\widehat{H}~+~\ln\left[\frac{N}{\Gamma}\right]~ +~\int{F[\Omega]~{{\rm d}}\Omega}~=~\textrm{const.}~,$$ (5)

where the function $\widehat{H}$ is defined by

$$\widehat{H}~\widehat{=}~\int{\frac{{{\rm d}} P}{e[P]~+~P}},$$ (6)

with suitable boundary conditions. In this work, we have rescaled the temperature and composition profiles assuming, at each evolutionary time, that they depend only on the baryon density $n_{\rm B}$. In this way, the EOS behaves effectively as a barotrope, and the pseudo-enthalpy function (6) can be written as

$$\widehat{H}~\widehat{=}~\int_0^{n_{\rm B}} {\frac{{{\rm d}} ... ...}}~ \frac{1}{e[n_{\rm B}]~+~P[n_{\rm B}]}{{\rm d}} n_{\rm B}},$$ (7)

with the effective functions $P[n_{\rm B}]$ and $e[n_{\rm B}]$ depending only on the baryonic density $n_{\rm B}$.

At zero temperature, $\widehat{H}$ is equal to the logarithm of the specific enthalpy

$$H~=~\ln\left[\frac{e~+~P}{n_{\rm B}~m_{\rm B}}\right],$$ (8)

where $m_{\rm B}$ is the baryon mass. As a basic test of the numerical method of integration, we have checked, in the case of tabulated EOS at T = 0, the agreement between  $\widehat{H}$ calculated numerically and the analytical result (8). Due to the finite number of points in the tabular EOS the agreement is not complete, but the discrepancy was always found to be smaller than $1\%$ in all tested cases (even with strong differential rotation).

   
3 Equation of state and thermodynamics

Most of previous studies of structure of rotating neutron stars have only been performed with polytropes, or at most zero temperature EOS. This is not enough for the purpose of studying PNS evolution, since thermal and chemical effects are relevant and need a careful treatment. Only in a few works (Goussard et al. 1997,1998; Strobel et al. 1999) a realistic, finite temperature EOS has been used. The works by Goussard et al. (1997,1998) used a Skyrme-like potential model (Lattimer & Swesty 1991), while the more recent work by Strobel et al. (1999) used a Thomas-Fermi model with momentum and density dependent effective nucleon-nucleon interaction developed by Myers & Swiatecky (1990).

In this paper, we employ an EOS in the framework of relativistic field theory, in which the interaction between nucleons is mediated by the exchange of three meson fields. The coupling constants are chosen to reproduce some bulk properties of nuclear matter such as the nuclear saturation density, symmetry energy or nuclear compressibility. The formalism is widely used, and to avoid unnecessary repetition we refer to the complete review by Prakash et al. (1997) and references therein. The particular choice of the coupling constants is detailed in Pons et al. (1999).

	Figure 1: Temporal variations of the temperature ( top), entropy ( middle) and lepton fraction ( bottom) as functions of the baryon density, for evolutionary times corresponding to  t=0.1, 0.5, 1 and 8 s after formation of the PNS. Results are for a $M_{\rm B}=1.6~M_\odot$ star in spherical symmetry.
Open with DEXTER

Figure 1 shows profiles from 1D simulations of a $M_{\rm B}=1.6~M_\odot$ spherical star from Pons et al. (1999). In previous works (Goussard et al. 1997,1998; Strobel et al. 1999) the evolution of the star was approximated by choosing either simplified profiles of constant entropy or temperature, or the combination of a constant low entropy core (s=1-2) surrounded by a high entropy mantle (s=4-5) to mimic the situation in the early Kelvin-Helmholtz stage. It was shown that the choice of the critical density that separates the hot envelope from the inner core affected substantially the results. To reduce the uncertainties associated to the choice of this critical density and the fact that PNSs are not isothermal or isentropic configurations, we have chosen to use the thermodynamical structure from 1D simulations as a better approximation to reality, and assume that at each time the temperature and lepton fraction have the same dependence on the baryon density as the non-rotating models, $T=T[n_{\rm B}]$ and $Y_L=Y_L[n_{\rm B}]$. More explicitly: given the results from 1D simulations, the profiles of density, temperature, lepton fraction, entropy and pressure are tabulated at each time step, and these tables are taken as effective one-parameter EOSs where the independent variable is the baryon density ($n_{\rm B}$) and the rest of variables are functions of only $n_{\rm B}$.

Notice that, as we explained in the previous section, this (pseudo-barotropicity) assumption is a sufficient condition to ensure that the equation of motion has a first integral, making affordable the computation of the models. Of course, to know the exact profiles consistently, one should solve the neutrino transport equations in 2D along the whole evolution, but this requires a large technical and computational effort that is out of the scope of this paper.

   
4 Choice of the initial models

One of the first issues that must be addressed to study the early evolution of neutron stars is the initial model. In Pons et al. (1999) it was shown that for spherical models the overall evolution is not very sensitive to modifications of the initial thermodynamical profiles, while the total mass of the star is the parameter that mostly affects the subsequent evolution. It seems therefore natural to chose a canonical star with a baryon mass 1.6 $M_\odot$, that corresponds to a gravitational mass of the old, cold configuration of about 1.44 $M_\odot$. For completeness, we have also studied a model with a baryonic mass of 1.2 $M_\odot$. Notice that one cannot use the effective barotropes obtained for a PNSs with a given mass in the calculations of sequences with a different mass, because the thermodynamic properties are different.

Figure 2: Equatorial profiles of angular velocity ( top), and density ( bottom) of PNSs from axisymmetric simulations of stellar core collapse (DFM). The four models correspond to a different amount of angular momentum of the iron core, namely, $\vert T/W\vert= 0.9\%$ (solid), 0.5% (dashes), 0.25% (dash-dot), and 0.05% (dash-3 dots). The dotted lines on the upper panel are fits to the simple law (9) with R02=50 km2.

Open with DEXTER

	Figure 3: Surfaces of constant angular velocity (projected in a transversal plane) of a differentially rotating star with R0=100 km, for t=1 s, and rotating with $\Omega = 354$ rad/s. ( $\Omega /\Omega _{\rm K}=0.995$). The surface is marked by the thick solid line. Because of relativistic effects, the stationary solution does not show exact cylindrical symmetry.
Open with DEXTER

More difficult is to guess the particular rotation properties of a newly born neutron star. Newly born neutron stars are expected to rotate non rigidly, as opposite to old neutron stars in which the various viscous mechanisms had enough time to act. Even in the case that the iron core of the Supernova progenitor was rigidly rotating (a quite reasonable assumption, since it is an old massive star), the process of collapse would generate a significant amount of differential rotation. Indeed, recent calculations of stellar evolution indicate that the iron core of massive stars is rotating almost rigidly (Heger et al. 2000,2003) and this is the approximation used in the most recent simulations of rotating stellar core collapse (Müller et al. 2003). Based on this works, we have decided to study a few models from Dimmelmeier et al. (2002, DFM hereafter) and to run some more simulations of similar models of stellar core collapse varying the total initial angular momentum. Our models correspond to those labeled by A1 in DFM, i.e., the iron core of the progenitor is almost rigidly rotating. In Fig. 2, we show profiles of the angular velocity (upper panel) and density (lower panel) of the inner 150 km, as a function of the equatorial radius, immediately after the PNS is formed. The different lines correspond to models with a different amount of initial angular momentum, with values of the ratio between rotational to gravitational potential energy (|T/W|) of 0.9% (models B3 in DFM), 0.5% (B2), 0.25% (B1) and 0.05%, respectively. We included the last lower value because it is more realistic (Heger et al. 2003; Müller et al. 2003) than those studied in DFM. In the figure we also show with dots the simple fit

$$\Omega[r,\th]~=\Omega[r~\sin \th]~=~\frac{\Omega_{\rm c} R_{0}^2}{ R_{0}^2 + r^2 \sin^2 \theta}$$ (9)

which corresponds to the Newtonian limit of $\Omega$ consistent with Eq. (4). In this limiting case, we now see that R₀ is the value of the radius at which, in the equatorial plane ( $\th~\widehat{=}~\pi/2$), the angular velocity is half its limiting central value, $\Omega _{\rm c}$. Even if this remark is valid only in the Newtonian case, it is nicely verified for relativistic stars as well. For all four models we have taken  R₀² = 50 km², while  $\Omega _{\rm c}$ takes the values of 4500, 3500, 2600, and 1300 rad/s, from top to bottom. The agreement between the simple law (9) and the results from simulations is quite acceptable. In all the models the spatial scale of variation of the angular velocity (R₀) is the same because the density profiles were very similar (see Fig. 2, bottom). Varying the EOS might lead to a different density distribution, in such a way that one should expect smaller R₀ for softer EOSs (more compact PNS), and larger R₀ for stiffer EOSs. We must remark again that our assumption of having stationary axisymmetric solutions implied that the function F on the rhs of the equation of stationary motion (2) is only a function of the angular velocity. In the Newtonian case, this is equivalent to say that the angular velocity depends only on the distance to the axis ($r \sin\th$). In the relativistic case, there are small corrections but the overall distribution of angular velocity is nearly cylindrical, as can be seen in Fig. 3.

Table 1: Properties of non-rotating PNSs and rigidly rotating PNSs at the limiting frequency, for a fixed baryonic mass $M_{\rm B} = 1.6~M_{\hbox{$\odot$ }}$. The models are labeled by the evolutionary time of non-rotating PNSs at which the thermodynamical profile was calculated. The entries in the table are: central baryon number density, $n_{\rm c}$; gravitational mass, $M_{\rm G}$; circumferential equatorial radius, $R_{\rm eq}$; Kepler frequency, $\Omega _{\rm K}$; angular momentum, J; axis ratio r (polar to equatorial) and the rotation parameter |T/W| (rotational energy on absolute value of gravitational energy). For 4 time steps at early times we give results using the LS EOS at low densities for comparison (first four rows, between parenthesis).

In reality, immediately after core collapse and bounce, there is no reason to expect that the distribution of angular velocity corresponds to that of stationary equilibrium, and it may take several rotation periods for the star to relax to some stationary solution. A more detailed analysis of all these issues is out of the scope of this paper, and requires a parametric study of rotating core collapse simulations. However, for our purpose, which is just to understand in which range of parameters one should move when using Eq. (4), the conclusion is that a reasonable realistic model with differential rotation should be consistent with R₀ of the order of 10 km. Considering that there might be a number of processes that could lead to rigid rotation on a dynamical timescale, such as shear or magneto-rotational instabilities, we will study both cases, with substantial differential rotation and with rigid rotation. Notice that previous works (Goussard et al. 1998) on differentially rotating PNSs have focused on the parameter space region corresponding to  $R_0 \approx 1$ km, which seems too low according to our results.

Table 2: Properties of non-rotating PNSs and rigidly rotating PNSs at the limiting frequency, for a fixed baryonic mass $M_{\rm B} = 1.2~M_{\hbox{$\odot$ }}$. The entries in the table are the same as in Table 1.

A final remark concerns the geometry of the PNSs. Although for progenitors with moderate angular momentum, conservation of angular momentum leads to a characteristic profile of angular velocity that decreases with increasing distance to the rotation axis, in the case of very fast rotation, or unrealistic large differential rotation, the usual oblate shape of the PNS can become toroidal, with a maximum in density off axis. In this case the profile of the angular velocity can exhibit also a maximum at some distance from the axis. We will discuss below for which models a toroidal shape is found, but unfortunately our numerical algorithm does not allow us to compute equilibrium configurations with zero density at the origin.

   
5 Results

We have considered a number of thermodynamical profiles corresponding to snapshots of the early evolution of PNSs, and we have calculated the rotational properties in two different cases: rigidly rotating and differentially rotating PNSs. The left part of Table 1 shows some properties (central baryonic density, gravitational mass, and radius) of the non rotating models with a baryonic mass of 1.6 $M_\odot$. The models are labeled according to the evolutionary time (in s after birth) in the non-rotating case. Notice that for rotating PNSs the evolution time labeling the models is only an indication, since in a self-consistent simulation one would expect faster evolution. The reason of this faster evolution is that the lower densities, and consistently lower temperatures because the entropy is conserved, result in larger neutrino mean free paths. Therefore the neutrino transport would be more effective.

In order to investigate the influence of the low density (<10¹¹ g/cm³) EOS, we have also substituted the zero temperature BPS equation of state (Baym et al. 1971) by the LS (Lattimer & Swesty 1991) EOS in some models at early times. The results for the models with the LS EOS are shown between parenthesis in the first four rows of Table 1. At later times, after the mantle loses most of its thermal and lepton content and contracts, the influence of the low density EOS becomes less important and both EOSs give essentially the same result, because the amount of mass and the thickness of the low density layer is rather small. The left hand side of Table 2 shows the equivalent information but for a non rotating NS with a total baryonic mass of 1.2 $M_\odot$.

   
5.1 Rigidly rotating PNSs

In the right hand side of Tables 1 and 2 we summarize the properties of maximally rotating PNSs following the same thermodynamical sequence as in the non-rotating case. Focusing at early times (t<0.5 s) we see that the central densities of non-rotating and maximally rotating PNSs are not terribly different. The reason is that these objects are only slightly bound and they cannot rotate very fast, reaching the mass shedding limit at angular velocities of hundreds of Hz. Therefore, the rotational energy is not large and the models are not very different from the non-rotating case. This can also be seen in the values of the quantity |T/W|, which is about 0.03-0.05 for the BPS EOS and 0.06-0.07 when using the LS EOS at low density. The differences between these two EOSs come from the fact that the BPS is at zero temperature while the LS is at finite temperature. Due to the more extended mantle, the total angular momentum is systematically larger, of about 30$\%$, for the models with finite temperature at $\rho<10^{11}$ g/cm³.

After about 0.5-1.0 s, the mantle has contracted and lost its thermal and lepton content, so that the PNS becomes more compact and bound. This can be deduced by looking at the central density, which exhibits a sudden increase at the end of this first epoch. For the non-rotating models it happens between 0.5 and 1.0 s, while for the maximally rotating models it seems to be delayed until some time between 1.0 and 2.0 s, due to the role of the centrifugal force. For these more compact configurations, the limiting angular velocities are much larger (3000-6000 rad/s), and their central densities are considerably lower than those of the non-rotating PNSs. The ratio between polar and equatorial radii of maximally rotating PNSs turns out to be rather insensitive to the thermodynamics, being about 0.65 for early models and decreasing to 0.55 for late configurations. With respect to the |T/W| ratio, its maximum value is about 0.1 at the end of the evolution, still below the threshold to undergo dynamical instabilities. As a rule of thumb, we found that the maximum angular velocity can be estimated with good accuracy with the simple fit $\Omega_{\rm K} \approx 0.58~\sqrt{GM/R^3}$, valid for both masses (1.6 and 1.2 $M_\odot$) and all thermodynamical profiles studied. Despite our models are not close to the maximum mass for a given EOS, this result is in agreement with the empirical formula derived by Haensel & Zdunik (1989) relating the maximum rotation frequency to the mass and radius of the maximum mass static configuration. It seems that this empirical relation is valid in general.

The evolution of the angular momentum deserves a separate discussion. The first result is that including thermal effects at low density (LS) allows for significantly larger total angular momentum. However, this discrepancy disappears rapidly, and the maximum angular momentum decreases quite fast and monotonously with time, which has evolutionary implications. Whether or not neutrino transport can redistribute and take away a substantial fraction of the initial angular momentum is unclear and needs of multidimensional transport simulations to be clarified. In general, it is assumed that the angular momentum is approximately constant during the Kelvin-Helmholtz epoch, and it only varies on timescales larger than the diffusion timescale. If this was the case, PNSs born with maximal angular momentum cannot radiate it via neutrino diffusion with the same efficiency as the binding energy, and this may result in significant mass loss after reaching the mass shedding limit. Other possibility is that PNSs are not born maximally rotating, the angular momentum is (almost) conserved in a timescale of seconds, and the star speeds up as it contracts and evolves thermodynamically. We have simulated this scenario by taking a sequence of configurations with constant angular momentum ( $J=1.5~G~M_\odot^2/c$), always with fixed baryonic mass. Results are shown in Fig. 4. It can be seen how a model that initially had a |T/W| ratio of only 0.02 speeds up in a timescale of a second and then slowly approaches the mass shedding limit in a few seconds. We expect that fast (but not maximally) rotating PNSs might reach the threshold of some types of instabilities (CFS, bar-mode, ...) on a diffusion timescale, say within 10 s after birth. It will depend on the efficiency of neutrinos to carry away angular momentum, which is currently under investigation.

	Figure 4: Temporal evolution of the angular velocity $\Omega /\Omega _{\rm K}$ and the rotation parameter |T/W| for a sequence of rigidly rotating PNSs with constant angular momentum ( $1.5~GM_{\hbox{$\odot$ }}^2 /c$) and a fixed baryonic mass of  $M_{\rm B} = 1.6~M_{\hbox{$\odot$ }}$.
Open with DEXTER

The above discussion was focused on results for the canonical NS with a baryonic mass of 1.6 $M_\odot$. For completeness we also give results for a model with a lower baryonic mass of 1.2 $M_\odot$ in Table 2. The results are qualitatively similar to those of Table 1, with only some quantitative differences. Accordingly to the lower mass, the maximum angular velocities and |T/W| ratio are also lower.

	Figure 5: Iso-enthalpy lines of a differentially rotating star with R0=100 km, for t=1 s, and rotating with $\Omega = 354$ rad/s. ( $\Omega /\Omega _{\rm K}=0.995$). For low differential rotation the structure of the star still looks normal.
Open with DEXTER

   
5.2 Differentially rotating PNSs

As discussed in the previous section we will restrict to realistic values of R₀ in the range from 10 km (significant differential rotation) to 50 km (almost rigid rotation). A different scenario would be a hot NS resulting from the merger of two neutron stars. In this case, the distribution of angular momentum is more complex and depends not only on the initial total angular momentum, but also on a number of different parameters as the mass ratio or the relative orientation of the spins of the two stars. Lower values of R₀might be reached in this latter scenario, but above all, the rotation profile would probably be very different from Eq. (4), which is in agreement with a stationary configuration, the case of a merging being far from a stationary situation.

As before, the baryonic mass of the star is conserved and fixed to 1.6 $M_\odot$for all models. Within the range of parameters mentioned above, we have calculated a large number of models of differentially rotating PNSs with either a fixed R₀ of 10, 20, 50, and 100 km or varying R₀ according to the circumferential equatorial radius of the star ( $R_0~=~0.5~R_{\rm eq}$). In Figs. 5 and 6 we show iso-enthalpy surfaces for two representative models at an evolutionary time of t=1 s. Despite both configurations have similar equatorial radii ($\approx$30 km), they are qualitatively very different. In Fig. 5 the parameter R₀=100 km, i.e., the model has a very low degree of differential rotation. Even very close to the mass shedding limit (0.995 of the maximum angular velocity) the star shows a regular oblate profile. On the contrary, the model in Fig. 6 corresponds to the case of significant differential rotation (R₀=10 km). With strong differential rotation we have found that above some critical angular velocity (950 rad/s, for this model) the shape of the star becomes toroidal, with a maximum density off axis.

The most relevant results are presented in a compact form in Figs. 7 and 8, where we show the total angular momentum and the rotation parameter |T/W|as a function of the central angular velocity  $\Omega _{\rm c}$, for different values of R₀. For comparison, we have also included the results for rigid rotation (solid curve), together with those for R₀=50 km (dash-dotted line), R₀=20 km (thick dashes), $R_0=R_{\rm eq}/2$ (thin dashes) and R₀=10 km (dots). For $R_0\ge100$ km, the results are so close to the rigid rotation case that we do not show them in the figures.

	Figure 6: Iso-enthalpy lines of a differentially rotating star with R0=10 km, for t=1 s, and rotating with $\Omega = 1570$ rad/s. The toroidal shape is clearly visible for all models differentially rotating with R0 smaller than a certain threshold.
Open with DEXTER

	Figure 7: Total angular momentum as a function of the central angular velocity for PNSs with fixed baryonic mass of 1.6 $M_\odot$ at evolutionary times of 0.5 ( top), 1.0 ( middle) and 10 s ( bottom), and for different values of R0: $R_0=\infty$ (solid), 50 km (dash-dot), 20 km (thick dashed), $R_{\rm eq}/2$ (thin dashed), 10 km (dots).
Open with DEXTER

	Figure 8: Rotation parameter, |T/W|, as a function of the central angular velocity for PNS with fixed baryonic mass of 1.6 $M_\odot$ at evolutionary times of 0.5 ( top), 1.0 ( middle) and 10 s ( bottom), and for different values of R0: $R_0=\infty$ (solid), 50 km (dash-dot), 20 km (thick dashed), $R_{\rm eq}/2$ (thin dashed), 10 km (dots).
Open with DEXTER

The thermodynamical profiles for each panel are those corresponding to evolutionary times of the non-rotating models of t=0.5 (top), 1.0 (middle) and 10 s (bottom panel) after formation. From the results of the figure we can see that for R₀ of the order or larger than 50 km the effects of differential rotation are unimportant. Only for R₀ < 20 km and $R_0=0.5~R_{\rm eq}$ the differences are significant. As discussed in 4, results from axisymmetric core collapse seem to indicate that the scale in which one should expect variations of the angular velocity is of about 10 km, therefore the last two cases may be looked at as the most realistic, with the caveat of being constrained to the particular rotation law we used.

In the model corresponding to t=1 s it is visible a sudden increase of the angular momentum at a critical value of the angular velocity, different for each value of R₀. This is an effect of the phase transition from bulk matter to nuclei, and it is associated to the point when the angular velocity becomes so large that the central density is below nuclear saturation density ($\rho_0$). This transition is not observed at earlier times because matter at low density is still very hot, and the thermal effects suppress or smooth out the kink. At later times, the star is more compact and the central density is always above $\rho_0$, even close to the mass shedding limit. For the models with strong differential rotation, we had to stop our calculations before reaching the mass shedding limit because of numerical limitations associated to the special toroidal geometry that arises (see Fig. 6). These cases are marked by an asterisk in Figs. 7 and 8. In principle, there could be configurations with slightly larger angular momentum.

An expected result is that the central angular velocity can be up to a factor 5 to 10 larger in the case of differential rotation, but the maximum angular momentum is only about a fifty percent larger than for rigid rotation. More important is the variation in the maximum of |T/W|. While for rigidly rotating stars it is about 0.1, for PNSs with a large degree of differential rotation it can be as high as 0.2. The fact that differential rotation allows for the existence of equilibrium models at large |T/W| have also been found by Baumgarte et al. (2000) and Yoshida et al. (2002) for polytropes in the relativistic Cowling approximation using the same rotation law. A maximum value of |T/W|=0.2 is not enough to reach the dynamical instability threshold ($\approx$0.27), but it is sufficient to allow either the secular, gravitational wave driven, instability (Chandrasekhar 1970; Friedman & Schutz 1978) that happens at  0.1-0.15 or the recently proposed low |T/W| dynamical instabilities (Shibata et al. 2002,2003; Watts et al. 2003). Notice that despite our calculations cannot be carried up to the mass shedding limit for some models, we are not far from it and the numbers are close to the maximum value. This can be seen in Fig. 9, when we show the polar to equatorial axis ratio as a function of  $\Omega _{\rm c}$. For rigid rotation, or moderate differential rotation the calculations end at the Keplerian limit, while for strong differential rotation the calculation ends because of numerical problems (marked with asterisks). The slope of the curves is becoming more and more vertical, and the axis ratio is as low as 0.4, smaller than the minimum in the rigid case (0.6). This is an indication that the shape of the star is changing fast and only slightly larger angular velocities can be reached before the star starts loosing mass.

	Figure 9: Axis ratio (polar to equatorial) as a function of the central angular velocity for a PNS with a baryonic mass of 1.6 $M_\odot$ at t=10 s, and for different values of R0: $R_0=\infty$ (solid), 50 km (dash-dot), 20 km (thick dashed), $R_{\rm eq}/2$(thin dashed), 10 km (dots).
Open with DEXTER

The next question one can pose is whether or not dynamical instabilities arise in a real case. As for rigidly rotating stars, we now impose that the angular momentum is fixed along the evolution, as well as the baryonic mass, and study how the different variables evolve under such constraints. Figure 10 shows the temporal evolution of the central angular velocity (top) and of the |T/W| ratio (bottom) of a PNS differentially rotating with R₀=10 km and a fixed angular momentum of $J=1.5~G~M_\odot^2/c$. As in the case of rigid rotation, |T/W| increases as the star contracts and looses its binding energy in a neutrino diffusion timescale of about 10 s. Thus, if there are no significant losses of angular momentum, the various types of dynamical instabilities associated to critical values of the rotation parameter may arise not at the very beginning, but several seconds after formation.

	Figure 10: Temporal evolution of the central angular velocity $\Omega _{\rm c}$and of the rotation parameter |T/W| for two sequences of differentially rotating PNSs with R0=10 km (solid line) and R0=20 km (dashed line). In both cases we fix the total angular momentum to  $1.5~GM_{\hbox{$\odot$ }}^2 /c$ and the baryonic mass to $M_{\rm B} = 1.6~M_{\hbox{$\odot$ }}$.
Open with DEXTER

Figure 11: Estimated neutrino luminosity ( $-{\rm d}M_{\rm G}/{\rm d}t$) as a function of time, consistent with the evolution at fixed baryonic mass of $M_{\rm B} = 1.6~M_{\hbox{$\odot$ }}$ and angular momentum of  $1.5~GM_{\hbox{$\odot$ }}^2 /c$. The two models are differentially rotating with R0=10 km (solid line) and R0=20 km (dashed line). The order of magnitude and the exponential decay are similar to the luminosities obtained in simulations with neutrino transport for non-rotating stars.

Open with DEXTER

As a consistency check of our assumptions and many simplifications, one can make a simple estimate of the neutrino luminosity considering that the variation in the gravitational mass of the PNS is considered to be equal to the luminosity, i.e., $L_{\nu} = - \frac{{\rm d} M_G}{{\rm d}t}$. The result is shown in Fig. 11. Both the order of magnitude and the exponential decay of the luminosity in a timescale of 10 s are consistent with the results of detailed simulations with neutrino transport in spherical symmetry (Pons et al. 1999). It is interesting that the difference between the two models studied with R₀=10 km (solid) and R₀=20 km (dashes) is only visible after 10 s of evolution, when the stars become more compact. Therefore, imposing conservation of angular momentum and baryonic mass, and building a sequence keeping fixed R₀ is not violating any basic physical law, such as conservation of energy, and can be taken as a qualitative approach to the real case.

   
6 Conclusions

We have approached the problem of the evolution of rotating protoneutron stars by constructing evolutionary sequences of axisymmetric stationary configurations in General Relativity. The thermodynamical structure and evolution have been taken and extrapolated from simulations in spherical symmetry that included neutrino transport. Although this is a crude simplification, it already gives an interesting insight about how the different relevant quantities can evolve as the rotating PNS loses its lepton content and its excess binding energy, and contracts. Moreover, we have found that the luminosity estimates are not terribly different from what one expects.

A special effort has been made to understand in which space of parameters we should be in a realistic case. The biggest uncertainty concerns the rotation law that PNS have at birth. By analyzing results from simplified core collapse simulations, it seems that a typical scale for variations of the angular velocity is about 10 km, and that conservation of angular momentum during the collapse of a stellar core (initially rigidly rotating) does not seem to allow for angular velocities varying in a length-scale shorter than a few km. Less is known about the angular distribution, except that the most recent simulations show the presence of important meridional currents and some turbulent motion. For simplicity, we restricted ourselves to the stationary case. Stationarity implies a quasi-cylindrical distribution (with deviations due to relativistic corrections) of the angular velocity. This stage only can be reached after several dynamical and rotation periods, after the PNS had time to relax. Therefore one must be aware that the first 0.5 s are probably far from stationarity, but after that evolution proceeds in a quasi-stationary way, except for low velocity convective motions. From our study of quasi-stationary sequences we can draw a few interesting qualitative results.

i)

For rigidly rotating stars, the mass shedding limit is well approximated by the simple law  $\Omega_{\rm K} \approx 0.58 \sqrt{GM/R^3}$, despite complications in the EOSs. This empirical formula is valid for stars of different masses within a 5% accuracy for the EOSs we considered. We believe that this result is general and would be valid for other EOSs or stars with different masses.

ii)

For each instant in the evolution, stars with strong differential rotation can have 5 to 10 times larger central angular velocities, and accommodate about a fifty percent more angular momentum. The maximum specific angular momentum J/M varies between $(1{-}2)~G~M_\odot/c \approx (0.5{-}1) \times 10^{16}$ cm²/s, depending on the degree of differential rotation.

iii)

When $R_0 \le 20$ km there is a critical value of the angular velocity above which the shape of the star becomes toroidal, with a maximum density off axis. We found that this configurations are stationary solutions of the hydrodynamics equations, but it is unclear under which circumstances this geometry is stable against perturbations. This issue deserves a separate study.

iv)

The maximum value of |T/W| obtained in the case of differential rotation is about 0.2, while for rigid rotation this is $\approx$0.11. Thus differentially rotating PNSs might undergo the CFS instability, which arises at $\approx$0.14 and, in any case, the recently discussed low |T/W| instability (Shibata et al. 2002,2003; Watts et al. 2003) is plausible to happen.

v)

More interestingly, we found several situations in which, even if the initial model is not close enough to the critical value of |T/W|, as the star contracts in a neutrino diffusion timescale of 5-10 s, it speeds up entering the window of instability. An observational evidence of this effect could be a temporal shift of a few seconds between the neutrino luminosity peak and a gravitational wave burst in the event of a galactic Supernova. Ultimately, this depends on the initial amount of angular momentum, which is approximately equal to the angular momentum of the iron core of the progenitor. Recent stellar evolution calculations suggest that the specific angular momentum of the inner 1.7 $M_\odot$ of a 15 $M_\odot$ star can be as high as  $3\times10^{16}$ cm²/s if magnetic braking is neglected, or 10¹⁵ cm²/s if magnetic torques are included in the evolution (Heger et al. 2003). This corresponds to J/M in the range $J/M=(0.2-6)~G~M_\odot/c$. If the angular momentum happens to be in the upper region (J/M>2), centrifugal forces would stop the collapse before the PNS is formed. Intermediate values (J/M=1) may result in the formation of a rapidly rotating PNS that enters the instability region several seconds after birth. If magnetic braking is very effective during the evolution of a massive star, J/M < 0.5 and the PNS will be formed after collapse without reaching extreme values of the angular velocities and |T/W|.

The next natural step to improve this work is to include the possible mechanisms to transport angular momentum between the different layers of the star, that may involve neutrino transport, turbulent transport, magnetic fields, neutrino viscosity, convective motion and/or angular momentum losses by gravitational wave emission. Unless the star is born with almost maximal angular velocities, since the different instabilities can arise in a timescale of seconds some of all of these dissipative mechanisms can modify our current vision of PNS evolution.

Acknowledgements

We are grateful to H. Dimmelmeier for providing us with rotation profiles that we have used to constrain the parameters of the rotation law. We thank L. Rezzolla, J. Novak, S. Yoshida and D. Gondek-Rosinska for useful discussions and comments. This work has been supported by the EU Programme "Improving the Human Research Potential and the Socio-Economic Knowledge Base'' (Research Training Network Contract HPRN-CT-2000-00137) and the Spanish Ministerio de Ciencia y Tecnologia grant AYA 2001-3490-C02. JAP is supported by a Ramón y Cajal contract from the Spanish MCyT and LV benefited from the Jumelage PAN-CNRS Astronomy France-Poland.

References

• Baumgarte, T. W., Shapiro, S. L., & Shibata, M. 2000, ApJ, 528, L29 [NASA ADS] [CrossRef] (In the text)
• Baym, G., Pethick, C. J., & Sutherland, J. 1971, ApJ, 170, 299 [NASA ADS] [CrossRef] (In the text)
• Bonazzola, S., Gourgoulhon, E., Salgado, M., & Marck, J.-A. 1993, A&A, 278, 421 [NASA ADS] [MathSciNet] (In the text)
• Bonanno, A., Rezzolla, L., & Urpin, V. 2003, A&A, 410, L33 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Burrows, A., & Lattimer, J. M. 1986, ApJ, 307, 178 [NASA ADS] [CrossRef] (In the text)
• Carter, B. 1969, J. Math. Phys., 10, 70 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Chandrasekhar, S. 1970, Phys. Rev. Lett., 24, 611 [NASA ADS] [CrossRef] (In the text)
• Dimmelmeier, H., Font, J. A., & Müller, E. 2002, A&A, 393, 523 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Friedman, J. L., & Schutz, B. F. 1978, ApJ, 221, 937 [NASA ADS] [CrossRef] (In the text)
• Gourgoulhon, E., Haensel, P., Livine, R., et al. 1999, A&A, 349, 851 [NASA ADS] (In the text)
• Goussard, J. O., Haensel, P., & Zdunik, J. L. 1997, A&A, 321, 822 [NASA ADS] (In the text)
• Goussard, J. O., Haensel, P., & Zdunik, J. L. 1998, A&A, 330, 1005 [NASA ADS]
• Haensel, P., & Zdunik, J. L. 1989, Nature, 340, 617 [NASA ADS] [CrossRef] (In the text)
• Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] (In the text)
• Heger, A., Woosley, S. E., Langer, N., & Spruit, H. 2003, Stellar Collapse, Ap&SS, ed. C. L. Fryer, in press [astro-ph/0301374]
• Kazana, D. 1977, Nature, 267, 501 [NASA ADS] (In the text)
• Keil, W., & Janka, H.-Th. 1995, ApJ, 295, 146 [NASA ADS] (In the text)
• Keil, W., Janka, H.-Th., & Müller, E. 1996, ApJ, 473, L111 [NASA ADS] [CrossRef] (In the text)
• Komatsu, H., Eriguchi, Y., & Hachisu, I. 1989, MNRAS, 239, 153 [NASA ADS] (In the text)
• Lattimer, J. M., & Swesty, F. D. 1991, Nucl. Phys., A535, 331 [NASA ADS] (In the text)
• Miralles, J. A., Pons, J. A., & Urpin, V. 2002, ApJ, 574, 356 [NASA ADS] [CrossRef] (In the text)
• Müller, E., Rampp, M., Buras, R., Janka, H.-Th., & Shoemaker, D. H. 2003, ApJ, submitted [astro-ph/0309833] (In the text)
• Myers, W. D., & Swiatecky, W. J. 1990, Ann. Phys. (NY), 204, 401 [NASA ADS] (In the text)
• Nozawa, T., Stergioulas, N., Gourgoulhon, E., & Eriguchi, Y. 1998, A&ASS, 132, 431 [NASA ADS] (In the text)
• Pons, J. A., Reddy, S., Prakash, M., Lattimer, J. M., & Miralles, J. A. 1999, ApJ, 513, 780 [NASA ADS] [CrossRef] (In the text)
• Pons, J. A., Miralles, J. A., Prakash, M., & Lattimer, J. M. 2000, ApJ, 553, 382 [NASA ADS] (In the text)
• Pons, J. A., Steiner, A. W., Prakash, M., & Lattimer, J. M. 2001, Phys. Rev. Lett., 86, 5223 [NASA ADS] [CrossRef]
• Prakash, M., Bombaci, I., Prakash, M., et al. 1997, Phys. Rep., 280, 1 [NASA ADS] (In the text)
• Romero, J. V., Alonso, J. D., Ibáñez, J. M., Miralles, J. A., & Pérez, A. 1992, ApJ, 395, 612 [NASA ADS] [CrossRef] (In the text)
• Shibata, M., Karino, S., & Eriguchi, Y. 2002, MNRAS, 334, L27 [NASA ADS] [CrossRef] (In the text)
• Shibata, M., Karino, S., & Eriguchi, Y. 2003, MNRAS, 343, 619 [NASA ADS] [CrossRef]
• Strobel, K., Schaab, Ch., & Weigel, M. K. 1999, A&A, 350, 497 [NASA ADS] (In the text)
• Sumiyoshi, K., Ibáñez, J. M., & Romero, J. V. 1999, A&ASS, 134, 39 [NASA ADS] (In the text)
• Swesty, F. D. 1996, J. Comp. Phys., 127, 118 [NASA ADS] (In the text)
• Watts, A. L., Andersson, N., & Jones, D. I. 2003, MNRAS, submitted [astro-ph/0309554] (In the text)
• Yoshida, S., Rezzolla, L., Karino, S., & Eriguchi, Y. 2002, ApJ, 568, L41 [NASA ADS] [CrossRef] (In the text)
• Yuan, Y., & Heyl, J. S. 2003 [astro-ph/0305083] (In the text)