Single Stars

6.1 Single Stars

The numerical study of the evolution of stars in general relativity involves two parts: first, one must solve for the evolution of matter (relativistic hydrodynamics, see [77

]), and second, one must compute the new configuration of the gravitational field. Whereas spectral-methods based codes are now able to study quite well the second part (see Section 6.2), the first has not benefited from the efforts of groups using spectral methods in the past decade. Thus, one faces a paradox: spectral methods have been primarily developed for the simulation of hydrodynamic systems (see Section 1.2), but they are not often used for relativistic hydrodynamics. This might be understood as a consequence of the general problem of spectral methods to deal with discontinuous fields and supersonic flows: the Gibbs phenomenon (see Section 2.4.4). Relativistic flows in astrophysics are often supersonic and therefore contain shocks. Although some techniques have been devised to deal with them in one-dimensional studies (see, e.g., [45 ]), there has been no convincing multidimensional convincing work. Another problem of multidimensional relativistic hydrodynamics that can spoil rapid convergence properties is sharp density profiles near neutron star surfaces. These can imply a diverging or discontinuous radial derivative of the density, thus slowing down the convergence of the spectral series.

6.1.1 Supernova core collapse

The physical scenario studied here is the formation of a neutron star from the gravitational collapse of a degenerate stellar core. This core can be thought of as the iron core of a massive star at the end of its evolution (standard mechanism of type II supernova). The degeneracy pressure of the electrons can no longer support the gravity and the collapse occurs. When the central density reaches nuclear values, the strong interaction stiffens the equation of state, stopping the collapse in the central region and generating a strong shock. This mechanism has long been thought to be a powerful source of gravitational radiation, but recent simulations show that the efficiency is much lower than previously estimated [70 , 196 ]. The first numerical study of this problem was the spherically-symmetric approach by May and White [149 ] using artificial viscosity to damp the spurious numerical oscillations caused by the presence of shock waves in the flow solution. Currently, state-of-the-art codes use more sophisticated High-Resolution Shock-Capturing (HRSC) schemes or High-Resolution Central (HRC) schemes (for details on these techniques, see the review by Font [77 ]). The first axisymmetric fully (general) relativistic simulations of the core collapse scenario were done by Shibata [192 ] and Shibata and Sekiguchi [196 ], which used HRSC schemes and a parametric equation of state. More recently, magnetohydrodynamic effects have been taken into account in the axisymmetric core collapse by Shibata et al. [194 ] using HRC schemes. Three-dimensional core collapse simulations, including a more realistic equation of state and deleptonization scheme have been performed within the cactus-carpet-whisky [148 , 17 ] framework by Ott et al. [165 , 164 ]. These simulations have been compared with those of the CoCoNuT code (see [68 , 69 ] and later in this section). A more detailed historical presentation can be found in the Living Review by Fryer and New [84 ].

The appearance of a strong hydrodynamic shock is, in principle, a serious problem to numerical models using spectral methods. Nevertheless, a first preliminary study in spherical symmetry and the Newtonian theory of gravity was undertaken in 1986 by Bonazzola and Marck [43 ], with the use of “natural” viscosity. The authors show a mass conservation to a level better than 10^–4 using one domain with only 33 Chebyshev polynomials. In 1993, the same authors performed the first three-dimensional simulation (still in Newtonian theory) of the pre-bounce phase [46 ], giving a computation of the gravitational wave amplitude, which was shown to be lower than standard estimates. Moreover, they showed that for a given mass, the gravitational wave amplitude depends only on the deformation of the core. These three-dimensional simulations were made possible thanks to the use of spectral methods, particularly for the solution of the Poisson equation for the gravitational potential.

Thus, shock waves pose a problem to spectral codes and have either been smoothed with spectral vanishing viscosity [112 ] or ignored by the code stopping before their appearance. Another idea developed first between the Meudon and Valencia groups was to use more appropriate techniques for the simulation of shock waves: namely the High-Resolution Shock-Capturing techniques, also known as Godunov methods (see the Living Reviews by Martí and Müller [145 ] and by Font [77 ]). On the other hand, one wants to keep the fewest degrees of freedom required by spectral methods for an accurate-enough description of functions, in particular for the solution of elliptic equations or for the representation of more regular fields, such as gravitational ones. Indeed, even in the case where a hydrodynamic shock is present, since it only appears as a source for the metric in Einstein’s equations, the resulting gravitational field is at least $𝒞1$ and the spectral series do converge, although slower than in the smooth case. Moreover, in a multidomain approach, if the shock is contained within only one such domain, it is then necessary to increase resolution in only this particular domain and it is still possible to keep the overall number of coefficients lower than the number of points for the HRSC scheme. The combination of both types of methods (HRSC and spectral) was first achieved in 2000 by Novak and Ibáñez [159 ]. They studied a spherically-symmetric core collapse in tensor-scalar theory of gravity, which is a more general theory than general relativity and allows a priori for monopolar gravitational waves. The system of PDEs to be solved resembles that of general relativity, with the addition of a scalar nonlinear wave equation for the monopolar dynamic degree of freedom. It was solved by spectral methods, whereas the relativistic hydrodynamics equations were solved by Godunov techniques. Two grids were used, associated to each numerical technique, and interpolations between the two were done at every timestep. Although strong shocks were present in this simulation, they were sharply resolved with HRSC techniques and the gravitational field represented through spectral methods did not exhibit any Gibbs-like oscillations. Monopolar gravitational waveforms could, thus, be given. In collaboration with the Garching relativistic hydrodynamics group, this numerical technique was extended in 2005 to three dimensions, but in the conformally-flat approximation of general relativity (see Sections 5.5 and 5.6) by Dimmelmeier et al. [69 ]. This approach using spectral methods for the gravitational field computation is now sometimes referred to as the “Marriage des Maillages” (French for “grid wedding”) and is currently employed by the core-collapse code CoCoNuT of Dimmelmeier et al. [68 , 69 ] to study general relativistic simulations of protoneutron stars, with a microphysical equation of state as well as an approximate description of deleptonization [70 ].

6.1.2 Collapse to a black hole

Stellar collapse to a black hole has been widely studied, starting with spherically-symmetric computations; in the case of dust (matter with no pressure), an analytical solution by Oppenheimer and Snyder [163 ] was found in 1939. Pioneering numerical works by Nakamura and Sato [154 , 155 ] studied the axisymmetric general relativistic collapse to a black hole; Stark and Piran [204 ] gave the gravitational wave emission from such a collapse in the formalism of Bardeen and Piran [19 ]. Fully general relativistic collapse simulations in axisymmetry have also been performed by Shibata [191 ], and the first three-dimensional calculations of gravitational-wave emission from the collapse of rotating stars to black holes was done by Baiotti et al. [17 ]. Recently, Stephens et al. [205 ] developed an evolution code for the coupled Einstein–Maxwell-MHD equations, with application to the collapse to a black hole of a magnetized, differentially-rotating neutron star.

To our knowledge, all studies of the collapse to a black hole, which use spectral methods are currently restricted to spherical symmetry. However, in this case and contrary to the core-collapse scenario, there is a priori no shock wave appearing in the evolution of the system and spectral methods are highly accurate at modeling the hydrodynamics as well. Thus, assuming spherical symmetry, the equations giving the gravitational field are very simple, first because of Birkhoff’s theorem, which gives the gravitational field outside the star, and then from the absence of any dynamic degree of freedom in the gravitational field. For example, when choosing the radial (Schwarzschild) gauge and polar slicing, Einstein’s equations, expressed within 3+1 formalism, turn into two first-order constraints, which are simply solved by integrating with respect to the radial coordinate (see [95 ]).

In the work of Gourgoulhon [95 ] a Chebyshev-tau method is used. The evolution equations for the relativistic fluid variables are integrated with a semi-implicit time scheme and a quasi-Lagrangian grid: the boundary of the grid is comoving with the surface of the star, but the grid points remain the usual Gauss–Lobatto collocation points (Section 2.3.2). Due to the singularity-avoiding gauge choice, the collapsing star ends in the “frozen-star” state, with the collapse of the lapse. This induces strong gradients on the metric potentials, but the code is able to follow the collapse down to very small values of the lapse, at less than 10^–6. The code is very accurate at determining whether a star at equilibrium is unstable, by triggering the physical instability from numerical noise at very low level. This property was later used by Gourgoulhon et al. [102 ] to study the stability of equilibrium configurations of neutron stars near the maximal mass, taking into account the effect of weak interaction processes. The addition of an inward velocity field to the initial equilibrium configurations enabled Gourgoulhon [96 ] to partially answer the question of the minimal mass of black holes: can the effective mass-energy potential barrier associated with stable equilibrium states be penetrated by stars with substantial inward radial kinetic energy? In [96 ], Gourgoulhon was able to form a black hole with a starting neutron star that was 10% less massive than the usual maximal mass.

The spectral numerical code developed by Gourgoulhon [95 ] has been extended to also simulate the propagation of neutrinos, coming from the thermal effect and nonequilibrium weak interaction processes. With this tool, Gourgoulhon and Haensel [101 ] have simulated the neutrino bursts coming from the collapse of neutron stars, with different equations of state. Another modification of this spectral code has been done by Novak [157 ], extending the theory of gravity to tensor-scalar theories. This allowed for the simulation of monopolar gravitational waves coming from the spherically-symmetric collapse of a neutron star to a black hole [157 ]. From a technical point of view, the solution of a nonlinear wave equation on curved spacetime has been added to the numerical model. It uses an implicit second-order Crank–Nicolson scheme for the linear terms and an explicit scheme for the nonlinear part. In addition, as for the hydrodynamics, the wave equation is solved on a grid, partly comoving with the fluid. The evolution of the scalar field shows that the collapsing neutron star has “expelled” all of its scalar charge before the appearance of the black hole.

6.1.3 Relativistic stellar pulsations

Oscillations of relativistic stars are often studied as a time-independent, linear eigenvalue problem [130 ]. Nevertheless, numerical approaches via time evolutions have proved to bring interesting results, as obtained by Font et al. [78 ] for the first quasiradial mode frequencies of rapidly-rotating stars in full general relativity. Nonlinear evolution of the gravitational-radiation–driven instability in the $r$ -modes of neutron stars has been studied by many authors (for a presentation of the physical problem, see Section 13 of [5 ]). In particular, the first study of nonlinear $r$ -modes in rapidly-rotating relativistic stars, via three-dimensional general-relativistic hydrodynamic evolutions has been done by Stergioulas and Font [207 ]. Different approaches to numerical hydrodynamic simulations in Newtonian gravity have been attempted by Lindblom et al. [139 ], with an additional braking term, as by Villain and Bonazzola [224 ] (see the following).

Because of their very high accuracy, spectral methods are able to track dynamic instabilities in the evolution of equilibrium neutron star configurations, as shown in section 6.1.2 by the work of Gourgoulhon et al. [95 , 102 ]. In this work, when the initial data represents a stable neutron star, some oscillations appear, which corresponds to the first fundamental mode of the star. As another illustration of the accuracy, let us mention the work of Novak [156 ], who followed the dynamic evolution of unstable neutron stars in the tensor-scalar theory of gravity. The instability is linked with the possibility for these stars of undergoing some “spontaneous scalarization”, meaning that they could gain a very high scalar charge, whereas the scalar field would be very weak (or even null) outside the star. Thus, for a given number of baryons, there would be three equilibria for a star: two stable ones with high scalar charges (opposite in sign) and an unstable one with a weak scalar charge. Starting from this last one, the evolution code described in [157 ] was able to follow the transition to a stable equilibrium, with several hundreds of damped oscillations for the star. This damping is due to the emission of monopolar gravitational waves, which carry away the star’s kinetic energy. The final state corresponds to the equilibrium configuration, independently computed by a simple code solving the generalized Tolman–Oppenheimer–Volkoff system with a scalar field, up to 1% error, after more than 50,000 timesteps. These studies could be undertaken with spectral methods because in these scenarios the flow remains subsonic and one does not expect any shock to be formed.

It is therefore quite natural to try and simulate stellar pulsations using spectral methods. Unfortunately, there have been only a few such studies, which are detailed in the following. Lockitch et al. [142 ] have studied the inertial modes of slowly-rotating stars in full general relativity. They wrote down perturbation equations in the form of a system of ordinary differential equations, thanks to a decomposition into vector and tensor spherical harmonics. This system is then a nonlinear eigenvalue problem for the dimensionless mode frequency in the rotating frame. Equilibrium and perturbation variables are then expanded in terms of a basis of Chebyshev polynomials, taking into account the coordinate singularity at the origin and parity requirements. The authors were therefore able to determine the values of the mode frequency, making the whole system singular and looked for eigenfunctions, through their spectral decomposition. They found that inertial modes were slightly stabilized by relativistic effects.

A different and maybe more natural approach, namely the time integration of the evolution equations, has been undertaken by Villain et al. [224 , 225 ] with a spectral magnetohydrodynamics code in spherical coordinates. The code solves the linearized Euler or Navier–Stokes equations, with the anelastic approximation. This approximation, which is widely used in other fields of astrophysics and atmospheric physics, consists in neglecting acoustic waves by assuming that time derivatives of the pressure and the density perturbations are negligible. It allows for a characteristic time, which is not set by acoustic propagation time, but is much longer and the timestep can be chosen so as to follow the inertial modes themselves. In their 2002 paper [224 ], Villain et al. study inertial modes (i.e. modes whose restoring force is the Coriolis force, among which the $r$ -modes [5 ]) in slowly rotating polytropes with $γ = 2$ in the linear regime. First, this is done in the framework of Newtonian gravity, where the anelastic approximation implies that the Eulerian perturbations of the gravitational potential do not play any role in the velocity perturbations. Second, they study the relativistic case, but with the Cowling approximation, meaning again that the metric perturbations are discarded. In both regimes and trying different boundary conditions for the velocity field at the surface of the star, they note the appearance of a polar part of the mode and the “concentration of the motion” near the surface, showing up in less than 15 periods of the linear $r$ -mode. A more recent work [225 ] deals with the study of gravity modes, in addition to inertial modes, in neutron stars. The interesting point of this work is the use of a quite realistic equation of state for nuclear matter, which is valid even when beta equilibrium is broken. The authors were, thus, able to show that the coupling between polar and axial modes is increasing with the rotation of the star, and that the coupling of inertial modes with gravity modes in nonbarotropic stars can produce fast energy exchanges between polar and axial parts of the fluid motion. From a numerical point of view, one of the key ingredients is the solution of the vector heat equation, coming from the Euler or Navier–Stokes equations. This is done by a poloidal-toroidal [47 ] decomposition of the velocity field into two scalar potentials, which is very natural in spectral methods. Moreover, to ensure correct analytical behavior at the origin, all scalar quantities are projected at each timestep to a modified Legendre function basis.

More recently, a complete nonlinear study of rotating star pulsations has been set by Dimmelmeier et al. [71 ]. They used the general relativistic code CoCoNuT (see above, Section 6.1.1) in axial symmetry, with an HRSC hydrodynamic solver, and spectral methods for the simplified Einstein equations (conformally-flat three metric). They noted that the conformal flatness condition did not have much effect on the dynamics when comparing with the Cowling approximation. Nevertheless, they found that differential rotation was shifting the modes to lower frequencies and they confirmed the existence of the mass-shedding–induced damping of pulsations.