Improved simulations of relativistic stellar core collapse

1. Improved simulations of relativistic stellar core collapse José A. Font Departamento de Astronomía y Astrofísica Universidad de Valencia (Spain) • Collaborators: • P. Cerdá-Durán, J.M. Ibáñez (UVEG) • H. Dimmelmeier, F. Siebel, E. Müller (MPA) • G. Faye (IAP), G. Schäfer (Jena) • J. Novak (LUTH-Meudon) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

2. Outline of the talk • Numerical simulations of rotational stellar core collapse: gravitational waveforms • Relativistic hydrodynamics equations in conservation form (Godunov-type schemes) • Approximations for the gravitational field equations (elliptic equations – finite-difference schemes, pseudo-spectral methods) • CFC (2D/3D) • CFC+ (2D) • Axisymmetric core collapse in characteristic numerical relativity • Improved means: • Treatment of gravity: from CFC to CFC+, and Bondi-Sachs • Modified CFC equations (high-density NS, BH formation) • Dimensionality: from 2D to 3D • Collapse dynamics: inclusion of magnetic fields Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

3. Astrophysical motivation General relativity and relativistic hydrodynamics play a major role in the description of gravitational collapse leading to the formation of compact objects (neutron stars and black holes):Core-collapse supernovae, black hole formation (and accretion), coalescing compact binaries (NS/NS, BH/NS, BH/BH), gamma-ray bursts. Time-dependent evolutions of fluid flow coupled to the spacetime geometry only possible through accurate, large-scale numerical simulations. Some scenarios can be described in the test-fluid approximation: hydrodynamical computations in curved backgrounds (highly mature nowadays). (see e.g. Font 2003 online article: relativity.livingreviews.org/Articles/lrr-2003-4/index.html). The (GR) hydrodynamic equations constitute a nonlinear hyperbolic system. Solid mathematical foundations and accurate numerical methodology imported from CFD. A “preferred” choice: high-resolution shock-capturing schemes written in conservation form. The study of gravitational stellar collapse has traditionally been one of the primary problems in relativistic astrophysics (for about 40 years now). It is a distinctive example of a research field in astrophysics where essential progress has been accomplished through numerical modelling with gradually increasing levels of complexity in the input physics/mathematics. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

4. Introduction: supernova core collapse in a nutshell The study of gravitational collapse of massive stars largely pursued numerically over the years. Main motivation in May and White’s 1967 first one-dimensional numerical relativity code. • Current standard model for a core collapse (type II/Ib/Ic) supernova: (from simulations! [Wilson et al (late 1980s), MPA, Oak Ridge, University of Arizona (ongoing)]) • Nuclear burning in massive star yields shell structure. Iron core with 1.4 solar masses and 1000 km radius develops in center. EoS: relativistic degenerate fermion gas, =4/3. • Instability due to photo-disintegration and e- capture. Collapse to nuclear matter densities in ~100ms. • Stiffening of EoS, bounce, and formation of prompt shock. • Stalled shock revived by neutrinos depositing energy behind it (Wilson 1985). Delayed shock propagates out and disrupts envelope of star. • Nucleosynthesis, explosion expands into interstellar matter. Proto-neutron starcools and shrinks to neutron star. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

5. Introduction (continued) Romero et al 1996 (radial gauge polar slicing). Purely hydrodynamical (prompt mechanism) explosion. No microphysics or -transport included! May & White’s formulation and 1d code used by many groups to study core collapse. Most investigations used artificial viscosity terms in the (Newtonian) hydro equations to handle shock waves. The use of HRSC schemes started in 1989 with the Newtonian simulations of Fryxell, Müller & Arnett (Eulerian PPM code). Basic dynamics of the collapse at a glance: 1d core collapse simulations Relativistic simulations of core collapse with HRSC schemes are still scarce. • Nonspherical core collapse simulations in GR very important: • To produce and extract gravitational waves consistently. • To explain rotation of newborn NS. • Collapse to NS is intrinsically relativistic (2M/R~0.2-0.4) (let alone to BH!) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

6. Multidimensional core collapse & gravitational waves Numerical simulations of stellar core collapse are nowadays highly motivated by the prospects of direct detection of the gravitational waves (GWs) emitted. GWs, ripples in spacetime generated by aspherical concentrations of accelerating matter, were predicted by Einstein in his theory of general relativity. Their amplitude on Earth is so small (about 1/100th of the size of an atomic nucleus!) that they remain elusive to direct detection (only indirectly “detected” in the theoretical explanation of the orbital dynamics of the binary pulsar PSR 1913+16 by Hulse & Taylor (Nobel laureates in physics in 1992). International network of resonant bar detectorsInternational network of interferometer detectors Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

7. Core collapse & gravitational waves (continued) • GWs are dominated by a burst associated with the bounce. If rotation is present, the GWs large amplitude oscillations associated with pulsations in the collapsed core(Mönchmeyer et al 1991; Yamada & Sato 1991; Zwerger & Müller 1997; Rampp et al 1998 (3D!)). • GWs from convection dominant on longer timescales(Müller et al 2004). • Müller (1982): first numerical evidence of the low gravitational wave efficiency of the core collapse scenario: E<10-6Mc2 radiated as gravitational waves. (2D simulations, Newtonian, finite-difference hydro code). • Bonazzola & Marck (1993): first 3D simulations of the infall phase using pseudo-spectral methods. Still, low amount of energy is radiated in gravitational waves, with little dependence on the initial conditions. • Zwerger & Müller (1997): general relativity counteracts the stabilizing effect of rotation. A bounce caused by rotation will occur at larger densities than in the Newtonian case •  need for relativistic simulations: • Dimmelmeier et al 2001, 2002; Siebel et al 2003; Shibata & Sekiguchi 2004, 2005; Cerdá-Durán et al 2005. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

8. Supernova codes vs core collapse numerical relativity codes • State-of-the art supernova codes are (mostly) based on Newtonian hydrodynamics(e.g. MPA group, Oak Ridge National Laboratory group). • Strong focus on microphysics (elaborate EoS, transport schemes for neutrinos – computationally challenging). • Often use of the most advanced initial models from stellar evolution. • Simple treatment of gravity (Newtonian, possibly relativistic corrections). However… no generic explosions yet obtained! (even with most sophisticated multi-dimensional models) • Core collapse numerical relativity codes (mostly) originate from vacuum Einstein codes (e.g. Whisky (EU), Shibata’s). • No microphysics: matter often restricted to ideal fluid EoS. • Simple initial (core collapse) models(uniformly or differentially rotating polytropes). • Exact or approximate Einstein equations for spacetime metric (inherit the usual complications found in numerical relativity: formulations of the field equations, gauge freedom, long-term numerical stability, etc). Our approach: flux-conservative hyperbolic formulation for the hydrodynamics CFC, CFC+, and Bondi-Sachs for the Einstein equations Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

9. 3+1 General Relativistic Hydrodynamics equations (1) Equations of motion: local conservation laws of density current (continuity equation) and stress-energy (Bianchi identities) Perfect fluid stress-energy tensor Introducing an explicit coordinate chart: • Wilson (1972) wrote the system as a set of advection equation within the 3+1 formalism. Non-conservative. • Conservative formulations well-adapted to numerical methodology are more recent: • Martí, Ibáñez & Miralles (1991): 1+1, general EOS • Eulderink & Mellema (1995): covariant, perfect fluid • Banyuls et al (1997): 3+1, general EOS • Papadopoulos & Font (2000): covariant, general EOS • Different formulations exist depending on: • The choice of time-slicing: the level surfaces of can be spatial (3+1) or null (characteristic) • The choice of physical (primitive) variables (, , ui …) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

10. 3+1 General Relativistic Hydrodynamics equations (2) Einstein’s equations Foliate the spacetime witht=constspatial hypersurfacesSt  n t Letnbe the unit timelike 4-vectororthogonalto St such that Eulerian observers u:fluid’s 4-velocity,p:isotropic pressure,r : rest-mass density e :specific internal energy density,e=r( 1+e ): energy density Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

11. 3+1 General Relativistic Hydrodynamics equations (3) Replace the “primitive variables” in terms of the “conserved variables” : First-order flux-conservative hyperbolic system Banyuls et al, ApJ, 476, 221 (1997) Font et al, PRD, 61, 044011 (2000) where is the vector of conserved variables fluxes sources Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

12. Nonlinear hyperbolic systems of conservation laws (1) For nonlinear hyperbolic systems classical solutions do not exist in general even for smooth initial data. Discontinuities develop after a finite time. For hyperbolic systems of conservation laws, schemes written inconservation formguarantee that the convergence (if it exists) is to one of the weak solutions of the original system of equations (Lax-Wendroff theorem 1960). A scheme written in conservation form reads: where is a consistentnumerical fluxfunction: Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

13. Nonlinear hyperbolic systems of conservation laws (2) The conservation form of the scheme is ensured by starting with the integral version of the PDE in conservation form. By integrating the PDE within a spacetime computational cell the numerical flux function is an approximation to the time-averaged flux across the interface: The flux integral depends on the solution at the numerical interfaces during the time step When a Cauchy problem described by a set of continuous PDEs is solved in adiscretized formthe numerical solution ispiecewise constant(collection of local Riemann problems). Key idea:a possible procedure is to calculate bysolving Riemann problemsat every cell interface (Godunov) Riemann solution for the left and right states along therayx/t=0. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

14. Nonlinear hyperbolic systems of conservation laws (3) Any FD scheme must be able tohandlediscontinuities in asatisfactoryway. • 1st order accurate schemes (Lax-Friedrich): Non-oscillatory but inaccurate across discontinuities (excessive diffusion) • (standard) 2nd order accurate schemes (Lax-Wendroff): Oscillatory across discontinuities • 2nd order accurate schemes with artificial viscosity • Godunov-type schemes (upwind High Resolution Shock Capturing schemes) Lax-Wendroff numerical solution of Burger’s equation at t=0.2 (left) and t=1.0 (right) 2nd order TVD numerical solution of Burger’s equation at t=0.2 (left) and t=1.0 (right) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

15. Nonlinear hyperbolic systems of conservation laws (4) rarefaction wave shock front Solution attime n+1of the two Riemann problems at the cell boundariesxj+1/2and xj-1/2 Spacetime evolution of the two Riemann problems at the cell boundariesxj+1/2 and xj-1/2.Each problem leads to a shock wave and a rarefaction wave moving in opposite directions (Piecewise constant) Initial data attime nfor the two Riemann problems at the cell boundariesxj+1/2and xj-1/2 cell boundaries where fluxes are required Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

16. Approximate Riemann solvers • In Godunov’s method the structure of the Riemann solution is “lost” in thecell averagingprocess (1st order in space). • Theexact solutionof a Riemann problem iscomputationally expensive, particularly in multidimensions and for complicated EoS. • Relativistic multidimensional problems: coupling of all flow velocity components through the Lorentz factor. • Shocks: increase in the number of algebraic jump (RH) conditions. • Rarefactions: solving a system of ODEs. Roe-type SRRS (Martí et al 1991; Font et al 1994) HLLE SRRS (Schneider et al 1993) Exact SRRS (Martí & Müller 1994; Pons et al 2000) Two-shock approximation (Balsara 1994) ENO SRRS (Dolezal & Wong 1995) Roe GRRS (Eulderink & Mellema 1995) Upwind SRRS (Falle & Komissarov 1996) Glimm SRRS (Wen et al 1997) Iterative SRRS (Dai & Woodward 1997) Marquina’s FF (Donat et al 1998) This motivated thedevelopment of approximate (linearized) Riemann solvers. Based on the exact solution of Riemann problems corresponding to a new system of equations obtained by a suitable linearization of the original one. The spectral decomposition of the Jacobian matrices is on the basis of all solvers. Approach followed by an important subset of shock-capturing schemes, the so-calledGodunov-type methods(Harten & Lax 1983; Einfeldt 1988). Martí & Müller, 2003 Living Reviews in Relativity www.livingreviews.org Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

17. A standard implementation of a HRSC scheme 1. Time update:Conservation form algorithm In practice: 2nd or 3rd order time accurate, conservative Runge-Kutta schemes (Shu & Osher 1989) 2.Cell reconstruction:Piecewise constant (Godunov), linear (MUSCL, MC, van Leer), parabolic (PPM, Colella & Woodward 1984)interpolation proceduresof state-vector variables from cell centers to cell interfaces. 3. Numerical fluxes:Approximate Riemann solvers (Roe, HLLE, Marquina). Explicit use of the spectral information of the system Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

18. HRSC schemes: numerical assessment Shock tube test Relativistic shock reflection • Stable and sharp discrete shock profiles • Accurate propagation speed of discontinuities • Accurate resolution of multiple nonlinear structures: discontinuities, raraefaction waves, vortices, etc V=0.99999c (W=224) Wind accretion onto a Kerr black hole (a=0.999M) Font et al, MNRAS, 305, 920 (1999) Simulation of a extragalactic relativistic jet Scheck et al, MNRAS, 331, 615 (2002) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

19. Relativistic Rotational Core Collapse (CFC) Dimmelmeier, Font & Müller, ApJ, 560, L163 (2001); A&A, 388, 917 (2002a); A&A, 393, 523 (2002b) Goals extend to GR previous results on Newtonian rotational core collapse (Zwerger & Müller 1997) determine the importance of relativistic effects on the collapse dynamics (angular momentum) compute the associatedgravitational radiation(waveforms) Model assumptions axisymmetry and equatorial plane symmetry (uniformly or differentially) rotating 4/3 polytropes in equilibrium as initial models (Komatsu, Eriguchi & Hachisu 1989). Central density 1010 g cm-3 and radius 1500 km. Various rotation profiles and rotation rates simplified EoS: P = Ppoly + Pth (neglect complicated microphysics and allows proper treatment of shocks) constrained system of the Einstein equations (IWM conformally flat condition) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

20. CFC metric equations In the CFC approximation (Isenberg 1985; Wilson & Mathews 1996) the ADM 3+1 equations reduce to a system of five coupled, nonlinear elliptic equations for the lapse function, conformal factor, and the shift vector: Solver 1: Newton-Raphson iteration. Discretize equations and define root-finding strategy. Solver 2: Conventional integral Poisson iteration.Exploits Poisson-like structure of metric equations, uk=S(ul). Keep r.h.s. fixed, obtain linear Poisson equations, solve associated integrals, then iterate until nonlinear equations converge. Both solvers feasible in axisymmetry but no extension to 3D possible. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

21. Animation of a representative rotating core collapse simulation For movies of additional models visit: www.mpa-garching.mpg.de/rel_hydro/axi_core_collapse/movies.shtml Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

22. Central Density Gravitational Waveform HRSC scheme: PPM + Marquina flux-formula Solid line:relativistic simulation Dashed line:Newtonian Larger central densities in relativistic models Similar gravitational radiation amplitudes (or smaller in the GR case) GR effects do not improve the chances for detection (at least in axisymmetry) “transition”Type II“multiple bounce”Type I“regular” Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

23. Gravitational Wave Signals www.mpa-garching.mpg.de/Hydro/RGRAV/index.html Influence of relativistic effects on signals: Investigate amplitude-frequency diagram Spread of the 26 modelsdoes not change much Signal of agalactic supernova detectable On average:Amplitude → Frequency ↑ If close to detection threshold: Signal could fall out of the sensitivity window! Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

24. CFC+ metric equations Cerdá-Duran, Faye, Dimmelmeier, Font, Ibáñez, Müller, and Schäfer, A&A, in press (2005) (ADM gauge) CFC+ metric: The second post-Newtonian deviation from isotropy is the solution of: (Schäfer 1990) (complicated) transverse, traceless projection operatorNewtonian potential Modified equations for , i and  (with respect to CFC): Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

25. CFC+ metric equations (2) We can solve the equations by introducing some intermediate potentials: 16 elliptic linear equations Linear solver: LU decomposition using standard LAPACK routines Boundary conditions Multipole development in compact-supported integrals Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

26. CFC+ results: rotating neutron stars Initial models (KEH method) Study the time-evolution of equilibrium models under the effect of a small amplitude perturbation. Computation of radial and quasi-radial mode-frequencies (code validation: comparison between CFC and CFC+ results, and with those of an independent full GR code) Equatorial profiles of the non-vanishing components of hij for the sequence of rigidly rotating models RNS0 to RNS5 Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

27. CFC+:radial modes of spherical NS quasi-radial modes of rotating NS spherical NS rotating NS No noticeable differences between CFC and CFC+ Good agreement in the mode frequencies (better than 2%), also with results from a full GR 3D code(Font et al 2002) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

28. CFC+: core collapse dynamics (1) Type I (regular collapse) Type III (rapid collapse) Relative differences between CFC and CFC+ for the central density and the lapse remain of the order of 10-4 or smaller throughout the collapse and bounce. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

29. CFC+: core collapse dynamics (2) Extreme case (torus-like structure) Type II (multiple bounce) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

30. CFC+: core collapse waveforms Two distinct ways to extract waveforms: From the quadrupole formula: From the metrichij: Offset correction (dashed line) Absolute differences between CFC and CFC+ waveforms. No significant differences found. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

31. “Mariage des Maillages”HRSC schemes for hydrodynamics and spectral methods for metric Reference: Dimmelmeier, Novak, Font, Müller, Ibáñez, PRD 71, 064023 (2005) The extension of our code to 3D has been possible thanks to the use of a metric solver based on integral Poisson iteration (as solver 2) but using spectral methods. MdM idea: Use spectral methods for the metric (smooth functions, no discontinuities) and HRSC schemes for the hydro (discontinuous functions). Valencia/Meudon/Garching collaboration. New metric solver uses publicly available package in C++ from Meudon group (LORENE). Communication between finite-difference grid and spectral grid necessary (high-order interpolators). It works! • Spectral solver uses several (3-6) radial domains (easy with LORENE package): • Nucleus limited by rd (domain radius parameter) roughly at largest density gradient. • Several shells up to rfd. • Compactified radial vacuum domain out to spatial infinity. In contraction phase of core collapse, inner domain boundaries are allowed to move (controlled by mass fraction or sonic point). The relation between the FD grid and the spectral grid changes dynamically due to moving domains. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

32. MdM code: Importance of a moving spectral grid In core collapse relevant radial scale contracts by a factor 100. Spectral grid setup with moving domains allows to put resolution where needed. Example: influence of bad spectral grid setup on collapse dynamics. • Domain radius rdmust follow contraction. • Domain radius rd should stay fixed at roughly rpns after core bounce. • More than 3 domains needed in dynamical core collapse. • Compare with previous solvers in axisymmetry: 33 collocation points per domain sufficient. rd held fixed (10% initial rse) wrong result! final rd too large bounce time only 3 radial domains Gibbs-type oscillations Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

33. MdM code: Oscillations of rotating neutron stars Another stringent test: can code keep rotating neutron stars in equilibrium? Test criterion: preservation of rotation velocity profile (here shown after 10 ms). Compactified grid essential if rfd close to rstar (profile deteriorates only negligibly) 3d low resolution2d high resolution 3d low resolution without artificial perturbation Axisymmetric oscillations in rotating neutron stars can be evolved as in other codes. No important differences between running the code in 2d or 3d modes. Proof of principle: code is ready for simulations of dynamical triaxial instabilities. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

34. MdM code: Generic nonaxisymmetric configurations Explore nonaxisymmetric configurations in 3d. Extension from axisymmetry to 3d trivial with LORENE. Even in axisymmetry spectral solver uses  coordinate with 4 collocation points (shift vector Poisson equation is calculated for Cartesian components). Setup: rotating NS with strong (unphysical) nonaxisymmetric “bar” perturbation. Rotation generates spiral arms Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

35. MdM code: Comparison with full GR core collapse simulations by Shibata and Sekiguchi Full GR simulations of axisymmetric core collapse available recently (Shibata & Sekiguchi 2004). Comparison between CFC and full GR possible! Shibata & Sekiguchi used rotational core collapse models with parameters close (but not equal) to the ones used by Dimmelmeier, Font & Müller (2002). Disagreements in the GW amplitude of about 20% at the peak (core bounce) and up to a factor 2 in the ringdown. Most plausible reason for discrepancy:different functional form of the density used in the wave extraction method (W6) and the formulation (stress formulation vs first moment of momentum density formulation). A3B2G4 (DFM 02)  W6 (20% gain at bounce!) A3B2G4 (A/rse=0.32) Shibata & Sekiguchi (A/rse=0.25) A3B2G4 The qualitative difference found by Shibata & Sekiguchi (2004) is due to the differences in the collapse initial model, notably the small decrease of the differential rotation length scale in their model. Shibata & Sekiguchi Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

36. CFC metric equations: modification to allow for black hole formation Original CFC equations It turns out to be essential to rescale some of the hydro quantities with the appropriate power of the conformal factor for the elliptic solvers to converge to the correct solution: To obtain one needs to first compute the conformal factor, which is obtained from the evolution equation Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

37. Black hole formation (spherical symmetry) with rescaling without rescaling High central density TOV solution Collapse of a (perturbed) unstable neutron star to a black hole in spherical symmetry. Collapse can be followed well beyond formation of an apparent horizon. Central density grows by 6 orders of magnitude, central lapse function drops to 0.0002. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

38. Rotational core collapse to high-density NS: CFC vs CFC+ • Model M7C5 (Shibata and Sekiguchi 2005): • Differential rotation parameter A/R=0.1 • Baryon rest mass M*=2.464 • Angular momentum J/M2=0.664 • Polytropic EOS (=4/3, k=7x1014 (cgs)) • Excellent agreement with the full BSSN simulations of Shibata & Sekiguchi (2005) max=1.4x1015 GW amplitude larger at bounce with CFC+ min=0.42 Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

39. Core collapse simulations using the Einstein equations for the Bondi metric Reference: Siebel, Font, Müller, and Papadopoulos, PRD 67, 124018 (2003) Ricci tensor Hypersurface equations: hierarchical set for Evolution equation for Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

40. The light-cone problem is formulated in the region of spacetime between a timelike worldtube at the origin of the radial coordinate and future null infinity. Initial data for  are prescribed on an initial light cone u=0. Boundary data for , U, V and  are also required on the worldtube. For the general relativistic hydrodynamics equations we use a covariant (form invariant respect to the spacetime foliation) formulation developed by Papadopoulos and Font (PRD, 61, 024015 (2000)) which casts the equations in flux-conservative, first-order form. • Gravitational waves at null infinity: • Bondi news function (from the metric variables expansion at scri) • Approximate gravitational waves (Winicour 1983, 1984, 1987): • Quadrupole news • First moment of momentum formula Null code test: time of bounce Time of bounce: 39.45 ms (null code 1), 38.32 ms (CFC code), 38.92 ms (null code 2). Good agreement between independent codes (less than 1% difference). Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

41. Gravitational waves: consistency & disagreement Good agreementin the computation of the GW strain using the quadrupole moment and the first moment of momentum formula. Equivalence valid in the Minkowskian limit and for small velocities, which explains the small differences. But Siebel et al (2002) found excellent agreement between the quadrupole news and the Bondi news when calculating GWs from pulsating relativistic stars. Quadrupole news rescaled by a factor 50. A possible explanation: different velocities involved in both scenarios, 10-5-10-4c for a pulsating NS and 0.2c in core collapse. Functional form for the quadrupole moment established in the slow motion limit on the light cone may not be valid. Large disagreement between Bondi news and quadrupole news, both in amplitude and frequency of the signal. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

42. 3+1 General Relativistic (Ideal) Magnetohydrodynamics equations (1) GRMHD: Dynamics of relativistic, electrically conducting fluids in the presence of magnetic fields. Ideal GRMHD: Absence of viscosity effects and heat conduction in the limit of infinite conductivity (perfect conductor fluid). The stress-energy tensor includes the contribution from the perfect fluid and from the magnetic field measured by the observer comoving with the fluid. with the definitions: Ideal MHD condition: electric four-corrent must be finite. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

43. 3+1 General Relativistic (Ideal) Magnetohydrodynamics equations (2) Adding all up:first-order, flux-conservative, hyperbolic system of balance laws + constraint(divergence-free condition) • Conservation of mass: • Conservation of energy and momentum: • Maxwell’s equations: • Induction equation: • Divergence-free constraint: Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

44. Solution procedure of the GRMHD equations Constrained transport scheme(Evans & Hawley 1988, Tóth 2000). Field components defined at cell interfaces. Zone-centered vector (needed for primitive recovery and cell reconstruction & Riemann problem) obtained from staggered field components: • Same HRSC schemes as for GRHD equations (HLL, Kurganov-Tadmor, Roe-type) • Wave structure information obtained • Primitive variable recovery more involved Details: Antón, Zanotti, Miralles, Martí, Ibáñez, Font & Pons, in preparation (2005) Update of field components: These equations conserve the discretization of Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

45. GRMHD equations: code tests (1) 1D Relativistic Brio-Wu shock tube test (van Putten 1993, Balsara 2001) Dashed line: wave structure in Minkowski spacetime at time t=0.4 Open circles: nonvanishing lapse function (2), at time t=0.2 Open squares: nonvanishing shift vector (0.4), at time t=0.16 HLL solver 1600 zones CFL 0.5 Agreement with previous authors (Balsara 2001) regarding wave locations, maximum Lorentz factor achieved, and numerical smearing of the solution. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

46. GRMHD equations: code tests (2) Magnetized spherical accretion onto a Schwarzschild BH density Test difficulty: keeping the stationarity of the solution Used in the literature (Gammie et al 2003, De Villiers & Hawley 2003) Initial data: Magnetic field of the type on top of the hydrodynamic (Michel) accretion solution. Radial magnetic field component chosen to satisfy divergence-free condition, and its strength is parametrized by the ratio: internal energy radial velocity HLL solver 100 zones =1 Solid lines: analytic solution Circles: numerical solution (t=350M) Increasing the grid resolution shows that code is second-order convergent irrespective of the value of  radial magnetic field Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

47. GRMHD equations: code tests (3) Magnetized equatorial Kerr accretion (Takahashi et al 1990, Gammie 1999) density Test difficulty: keeping the stationarity of the solution (algebraic complexity augmented, Kerr metric) Used in the literature (Gammie et al 2003, De Villiers & Hawley 2003) Inflow solution determined by specifying 4 conserved quantities: the mass flux FM, the angular momentum flux FL, the energy flux FE, and the component F of the electromagnetic tensor. azimuthal velocity a=0.5 FM=-1.0 FL=-2.815344 FE=-0.908382 F=0.5 HLL solver Solid lines: analytic solution Circles: numerical solution (t=200M) radial magnetic field azimuthal magnetic field second order convergence Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

48. GRMHD: spherical core collapse simulation • As a first step towards relativistic magnetized core collapse simulations we employ the test (passive) field approximation for weak magnetic field. • magnetic field attached to the fluid (does not backreact into the Euler-Einstein equations). • eigenvalues (fluid + magnetic field) reduce to the fluid eigenvalues only. HLL solver + PPM, Flux-CT, 200x10 zones The divergence-free condition is fulfilled to good precision during the simulation. The amplification factor of the initial magnetic field during the collapse is 1370. Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005

49. Summary of the talk • Multidimensional simulations of relativistic core collapse feasible nowadays with current formulations of hydrodynamics and Einstein’s equations. • Results from CFC and CFC+ relativistic simulations of rotational core collapse to NS in axisymmetry. Comparisons with full GR simulations show that CFC is a sufficiently accurate approach. • Modification of the original CFC equations to allow for collapse to high density NS and BH formation. • Ongoing work towards extending the SQF for GW extraction (1PN quadrupole formula). • Axisymmetric core collapse simulations using characteristic numerical relativity show important disagreement in the gravitational waveforms between the Bondi news and the quadrupole news. • 3d extension of the CFC core collapse code through the MdM approach (HRSC schemes for the hydro and spectral methods for the spacetime). • First steps towards GRMHD core collapse simulations (ongoing work) Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005