GR-hydro

Document Sample
GR-hydro Powered By Docstoc
					            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
                                          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
                                         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
          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 star cools
    and shrinks to neutron star.


    Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                                      Introduction (continued)
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:

1. To produce and extract gravitational waves
   consistently.
2. To explain rotation of newborn NS.
3. Collapse to NS is intrinsically relativistic
   (2M/R ~0.2-0.4)         (let alone to BH!)                       Romero et al 1996 (radial gauge polar slicing).
                                                                    Purely hydrodynamical (prompt mechanism) explosion. No
                                                                    microphysics or -transport included!


     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
          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 detectors              International network of interferometer detectors




      Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
               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-6 Mc2 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
            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
          3+1 General Relativistic Hydrodynamics equations (1)

                                                              Equations of motion:
        ( u )  0                            [1]
                                                               local conservation laws of density current
      T   0                                [4]            (continuity equation) and stress-energy (Bianchi
      p  p(  ,  )
                                                               identities)
                                                [1]

     Perfect fluid stress-energy tensor                               
                                                                    (  g u  )  0
           T         hu u  pg                                     x 
                                                                      
     Introducing an explicit coordinate chart:                                                
                                                                          (  g T  )   g T 
                                                                     x 
Different formulations exist depending on:              Wilson (1972) wrote the system as a set of advection equation
                                                        within the 3+1 formalism. Non-conservative.
1.    The choice of time-slicing: the level
      surfaces of     can be spatial (3+1) or           Conservative formulations well-adapted to numerical
      null (characteristic) x
                              0                         methodology are more recent:

2.    The choice of physical (primitive)                • Martí, Ibáñez & Miralles (1991): 1+1, general EOS
      variables (, , ui …)                            • Eulderink & Mellema (1995): covariant, perfect fluid
                                                        • Banyuls et al (1997): 3+1, general EOS
                                                        • Papadopoulos & Font (2000): covariant, general EOS

       Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
      3+1 General Relativistic Hydrodynamics equations (2)

   ( u  )  0                                                   1
                                          [1]                 R   g  R  8T [10]
  T   0                              [4]                       2
  p  p(  ,  )                          [1]                  Einstein’s equations


                                                                      Foliate the spacetime with t=const
                                                                      spatial hypersurfaces          St
                    
                                                                                                                         j
              n                                                    ds 2  ( 2    i )dt 2  2  dxi dt   dxi dx
                                                                                    i               i          ij
                     t
                                                                       Let n be the unit timelike 4-vector
                                                                       orthogonal to St such that

                                                                                                 1
                                                                                            n       ( t   i  i )
                                     n  i            1 u     
                                                                i
                                                                                                 
Eulerian observers            v              vi        t i
                                     nu                u
                                                         
                                                                
                                                                
u: fluid’s 4-velocity, p: isotropic pressure,  : rest-mass density
 : specific internal energy density, e=( 1+ ): energy density

   Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
          3+1 General Relativistic Hydrodynamics equations (3)
Replace the “primitive variables” in terms of the “conserved variables” :

                                             D  W                                       W 2  1 /(1  v j v j )
                                            
              
                       
              w   , vi ,               S j  hW 2 v j                                                  p
                                                                                              h  1  
                                             E  hW 2  p                                                   
          First-order flux-conservative hyperbolic system
                                                            i 
              1    u ( w)                             g f ( w)                            Banyuls et al, ApJ, 476,
                                                                      s ( w)                    221 (1997)
                 
              g     t                                   x i      
                                                                     
                                                                                                   Font et al, PRD, 61,
                                                                                                   044011 (2000)


              u ( w)  D, S j , E  D 
               
  where                                                      is the vector of conserved variables

  i        i i   i i                         i i         
            D v 
  f ( w)           , S j  v      p j , E  D v 
                                            i
                                                               pv 
                                                                   i                                      fluxes
                       
                                  
                                                                
                                                                 
                    g                        ln      0  
           0, T   j   gj ,   T  0
  
 s ( w)            x 
                            
                                                      T                                            sources
                                                               
                                                    
                                             x              
     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
       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 in conservation form
guarantee 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:



   n 1  n t   n  n
                   ˆ                      n           n
                                                      ˆ              n         n
  u j  u j  ( f (u j r , u j r 1 ,, u j  q )  f (u j r 1 , u j r ,, u j  q 1 ))
              x
        
        ˆ                                                   
                                                           ˆ                        
where f is a consistent numerical flux function: f (u , u ,, u )  f (u )
    Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
         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
                                    [ x j 1/ 2 , x j 1/ 2 ] to n , t n1 ]
the numerical flux function is an approximation [tthe time-averaged flux across the interface:



ˆ          1 t n1                                  The flux integral depends on the solution at the
                                                                                    
f j 1/ 2   n f (u ( x j 1/ 2 , t )) dt            numerical interfaces          u ( x j1/ 2 , t ) time step
                                                                                        during the
           t t
                                                           When a Cauchy problem described by a set of
                
Key idea: a possible procedure is to
                                                           continuous PDEs is solved in a discretized form the
calculate                / 2 t)
                u ( x j1by,solving                        numerical solution is piecewise constant (collection of
Riemann problems at every cell                             local Riemann problems).
interface (Godunov)

                        n n
  u ( x j 1/ 2 , t )  u (0; u j , u j 1 )

Riemann solution for the left and right
states along the ray x/t=0.


     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
        Nonlinear hyperbolic systems of conservation laws (3)
Any FD scheme must be able to handle discontinuities in a satisfactory way.

1.     1st order accurate schemes (Lax-Friedrich): Non-oscillatory but inaccurate across
       discontinuities (excessive diffusion)
2.      (standard) 2nd order accurate schemes (Lax-Wendroff): Oscillatory across discontinuities
3.     2nd order accurate schemes with artificial viscosity
4.     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
   Nonlinear hyperbolic systems of conservation laws (4)
                     rarefaction wave          shock front

                                                             Solution at time n+1 of the two Riemann
                                                             problems at the cell boundaries xj+1/2 and xj-1/2



                                                             Spacetime evolution of the two Riemann
                                                             problems at the cell boundaries xj+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 at time n for the
                                                             two Riemann problems at the cell boundaries
                                                             xj+1/2 and xj-1/2

                                                                   n 1  n t   n
                                                                                  ˆ           n 
                                                                                              ˆ
                                                                  uj  uj       f j 1/ 2  f j 1/ 2 
                                                                             x                        
   cell boundaries where fluxes are required


Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                                 Approximate Riemann solvers
In Godunov’s method the structure of the Riemann solution is “lost” in the cell averaging process (1st order in
space).
The exact solution of a Riemann problem is computationally 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.


This motivated the development of approximate                     Roe-type SRRS (Martí et al 1991; Font et al 1994)
(linearized) Riemann solvers.                                     HLLE SRRS (Schneider et al 1993)

Based on the exact solution of Riemann problems                   Exact SRRS (Martí & Müller 1994; Pons et al 2000)
corresponding to a new system of equations                        Two-shock approximation (Balsara 1994)
obtained by a suitable linearization of the original
                                                                  ENO SRRS (Dolezal & Wong 1995)
one. The spectral decomposition of the Jacobian
matrices is on the basis of all solvers.                          Roe GRRS (Eulderink & Mellema 1995)
                                                                  Upwind SRRS (Falle & Komissarov 1996)
Approach followed by an important subset of shock-
capturing schemes, the so-called Godunov-type                     Glimm SRRS (Wen et al 1997)
methods (Harten & Lax 1983; Einfeldt 1988).                       Iterative SRRS (Dai & Woodward 1997)
                                                                  Marquina’s FF (Donat et al 1998)
                                                                                      Martí & Müller, 2003
                                                                   Living Reviews in Relativity    www.livingreviews.org
     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                  A standard implementation of a HRSC scheme

1. Time update: Conservation form algorithm

      n1  n t   n
                     ˆ         n 
                               ˆ
     u j  u j   f j 1/ 2  f j 1/ 2 
                x                      
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)
3. Numerical fluxes: Approximate Riemann                      interpolation procedures of state-vector variables
solvers (Roe, HLLE, Marquina). Explicit use                   from cell centers to cell interfaces.
of the spectral information of the system


  1
 ˆ                                 5
                                        ~ ~ ~
 f i   f i ( wR )  f i ( wL )   n  n Rn 
      2                           n 1         
                                  5
        U( wR )  U( wL )   n Rn
                               ~ ~
                                 n 1


        Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                         HRSC schemes: numerical assessment
                                                  Relativistic shock reflection                       Shock tube test
• 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)




 Simulation of a extragalactic relativistic jet                Wind accretion onto a Kerr black hole (a=0.999M)
 Scheck et al, MNRAS, 331, 615 (2002)                          Font et al, MNRAS, 305, 920 (1999)

      Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                 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 associated gravitational 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
                               CFC metric equations                                      ij   4 ij
                                                                                           CFC


In the CFC approximation (Isenberg 1985; Wilson & Mathews 1996) the ADM 3+1 equations


      t ij  2K ij   i  j   j  i
                                                            
      t K ij   i  j   Rij  KK ij  2 K im K m   m m K ij  K im  j  m  K jm  i  m  8T ij
                                                     j

     R  K 2  K ij K ij  16 2T 00  0
                      
      i K ij   ij K  8S j  0
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
            
           5
                              K ij K ij 
  2 W  P   2                                            equations and define root-finding strategy.
                               16 
                                                                Solver 2: Conventional integral Poisson
                                                                  iteration. Exploits Poisson-like structure of metric
                                         7 K ij K ij 
             5
                           
   2  h 3W  2  5 P 
               
                          2
                                           16 
                                                                 equations, uk=S(ul). Keep r.h.s. fixed, obtain linear
                                                                  Poisson equations, solve associated integrals, then
                                                                iterate until nonlinear equations converge.

                                   1
 i  16  4 S i  2 K ij  j  6    i  k  k
                                   3
                                                                  Both solvers feasible in axisymmetry but no
                                                                  extension to 3D possible.
                                  
    Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
            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
                                                                                   Central Density Gravitational Waveform
HRSC scheme:
PPM + Marquina flux-formula




                                                       Type I “regular”
Solid line: relativistic simulation
Dashed line: Newtonian




                                                       Type II “multiple bounce”
Larger central densities in
    relativistic models
Similar gravitational radiation
    amplitudes (or smaller in the
    GR case)                                           “transition”




GR effects do not improve the
   chances for detection (at least
   in axisymmetry)


     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                                  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 models does not change much
Signal of a galactic 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
                                        CFC+ metric equations
         Cerdá-Duran, Faye, Dimmelmeier, Font, Ibáñez, Müller, and Schäfer, A&A, in press (2005)



CFC+ metric:         ij    ij  hij ,
                       CFC      CFC  TT
                                                                   tr TT  0
                                                                    h                  (ADM gauge)


The second post-Newtonian deviation from isotropy is the solution of:
                                   1 TTkl                                1
                     hij 
                       TT
                                       ij (16 vk vl  4 kU lU )   6                        (Schäfer 1990)
                                   c4                                    c 
                (complicated) transverse, traceless projection operator   Newtonian potential

Modified equations for , i and  (with respect to CFC):

                                                    K ij K ij 
                      2  hW  P 
                                       5  2                    
                                                     16 
                                                              
                                                                7 K ij K ij  1 TT
                                   
                                                   
                       2  5  h 3W 2  2  5 P           16  c
                                                                               2 hij  ijU
                                                                            
                                                       1
                     i  16  4 S i  2 K ij  j  6    i  k  k
                                                       3
                                                      
   Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                                     CFC+ metric equations (2)

                                                                                   S i  (4 vi v j   iU jU ) x j
                        TT
 We can solve the      hij
                        equations by introducing some
 intermediate potentials:
                                                                                   S ij  4 vi v j   iU jU
        1      7                                             1
hij 
 TT
          Sij   ij S kk  3x k  (i S j ) k  3 (i S j )  x j  i S kk         S        4 vi v j x i x j
        2      4                                             4
         1       1                 1                       1        1
         iT j  x k  ij S k  x k x l  ij S kl   ij S   i R j
                                                                                   T i    (4 v j v j   jU jU ) xi
         4       2                 4                       4        4              R i     i  kU lU x k x l

 S    v vk x  x U  x  kU d x 
    1                                    M2 i    1                               16 elliptic linear equations
  i        i    k   i      k        3
                                            n  O 2 
    r                                    2r      r                               Linear solver: LU decomposition
    1  i j  ij                      1                                         using standard LAPACK routines
 S    v v 
  ij
                  x j  iU d 3 x  O 2 
                             
    r          2                     r 
         1                          1
         r
                                                                   Boundary conditions
 S         vk vl x k x l d 3 x  O 2 
                                    r 
                                                                   Multipole development in compact-
            vk v k x i  x iU d 3 x 
                                                 1
                                            2
        1                                M i                       supported integrals
        r
 Ti                                        n  O 2 
                                         2r      r 
     M2 i   1
 R 
   i
       n  O 2 
     r      r 
       Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                           CFC+ results: rotating neutron stars
Initial models (KEH method)
                                  Mass
  Model   Axis ratio   /K
                                  (sun)

  RNS0      1.00        0.00       1.40

  RNS1       0.95       0.42       1.44
  RNS2       0.85       0.70       1.51
  RNS3       0.75       0.87       1.59
  RNS4       0.70       0.93       1.63
  RNS5       0.65       0.98       1.67


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
                                               Equatorial profiles of the non-vanishing components of hij for the
code)
                                               sequence of rigidly rotating models RNS0 to RNS5

     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
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
                         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
                         CFC+: core collapse dynamics (2)
      Type II (multiple bounce)                              Extreme case (torus-like structure)




Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                         CFC+: core collapse waveforms
                                                                                   Two distinct ways to extract
                                                                                   waveforms:

                                                                                   From the quadrupole formula:
                                                                                      quad                     AE 2 t  r / c 
                                                                                     h  x , t  
                                                                                                    1 15
                                                                                                         sin  2 20
                                                                                                    8                 r

                                                                                    From the metric hij:
                                                                                                                         AE2 t  r / c 
                                                                                    h PN x , t  r / c  
                                                                                                             1 15
                                                                                     2
                                                                                                                    sin 2 20
                                                                                                             64                 r
                                                                                                             1 quad 
                                                                                                            h x , t 
                                                                                                             8

                                                                                    Offset correction (dashed line)


                                                                                          h PN-corrected  h PN  a ij hij
                                                                                           2                 2              TT




                                                                                      Absolute differences
                                                                                      between CFC and CFC+
                                                                                      waveforms.

                                                                                      No significant differences
                                                                                      found.
Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                                  “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
           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.



          rd held fixed (10% initial rse)                                      1. Domain radius rd must follow
                                                                                  contraction.

                           wrong result!                                       2. Domain radius rd should stay fixed
                                                                                  at roughly rpns after core bounce.
                        final rd too large
                                                      bounce time              3. More than 3 domains needed in
                                                                                  dynamical core collapse.

                      only 3 radial domains                                    4. Compare with previous solvers in
                                                                                  axisymmetry: 33 collocation points
                                                                                  per domain sufficient.
                                                     Gibbs-type oscillations




   Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
             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 resolution            2d high resolution
Axisymmetric oscillations in rotating neutron stars
                                                                     3d low resolution without artificial perturbation
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
         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
             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                                                   A3B2G4 (DFM 02)
models with parameters close (but not equal) to the                              
ones used by Dimmelmeier, Font & Müller (2002).                                               W6 (20% gain at bounce!)
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                                                     A3B2G4 (A/rse=0.32)
extraction method (W6) and the formulation (stress                                                Shibata & Sekiguchi
formulation vs first moment of momentum density                                                     (A/rse=0.25)
formulation).


                           A3B2G4                         The qualitative difference found by Shibata &
                                                          Sekiguchi (2004) is due to the differences in the
                                                          collapse initial model, notably the small decrease of
                       Shibata & Sekiguchi                the differential rotation length scale in their model.


      Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                             CFC metric equations:
                  modification to allow for black hole formation
         Original CFC equations
                                                                 It turns out to be essential to rescale some
                             K K       ij
                                                                 of the hydro quantities with the appropriate
  2    W 2  P  ij
              5                       
                               16 
                                                                 power of the conformal factor for the elliptic
                                                               solvers to converge to the correct solution:
                                        7 K ij K ij 
               
                        
   2  5  h 3W 2  2  5 P       16 
                                                     
                                                                                  *   6
                                                    
                                                                                  S i*   6 S i
                                   1
 i  16  4 S i  2 K ij  j  6    i  k  k
                                   3                                           P*   6 P
                                  

                                                                 1              K ij K ij 
To obtain          * , S i* , P *
                               one
                                                                
                                                                        2
                                                                           
                                                         2    hW  P  
                                                                                5

                                                                                   16 
                                                                                             
                                                                                             
needs to first compute the                                                                  
conformal factor, which is
                                                                     1                                                
obtained from the evolution
equation                                                            
                                                                                
                                                          2    h 3W  5  5 P  
                                                                                        2
                                                                                                   
                                                                                                   
                                                                                                        5
                                                                                                           7 K ij K ij
                                                                                                              16
                                                                                                                         
                                                                                                                         
                                                                                                                        
                                                                                          1
              k  k                                   i  16   2 S i  2 K ij  j  6    i  k  k
           t 6                                                                              3
                                                                                            
    Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                   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
              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
                                      min=0.42                                             bounce with CFC+




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

            V                     
                                                                                                  
    ds 2   e 2   U 2 r 2 e 2 du 2  2e 2  dudr  2Ur 2 e 2 dud  r 2 e 2 d 2  e  2 sin 2 d 2                       
            r                     
The metric functionsV ,U ,  and  only depend on the coordinates u, r and 
            1                                      1              
Gab  Rab  g ab R            Rab    h uaub  g ab   pgab  Ricci tensor
                                                                    
            2                                      2              

                   Rrr   ,r   ,r 
                 r             r       2

                 4             2
                2r 2 Rr  r 4 e 2(   )U ,r ,r  2r 2   ,r   ,r  2 ,r ,   ,  2 ,r cot 
                                                                                       2                  
                                                                                       r                  
                                      1
     r 2 e 2  g AB RAB  2V,r  r 4 e 2(   ) (U ,r ) 2  r 2U ,r  4rU ,  r 2U ,r cot
                                      2
                                                         
                            4rU cot  2e 2 (   )  1  (3 ,   , ) cot   ,   ,  (  , ) 2  2 , ( ,   , )    
      r 2 e 2  g   R   2r (r ) ,ur  (1  r ,r )V,r  (r ,rr   ,r )V  r (1  r ,r )U ,
                             r 2 (cot   , )U ,r  e 2 (   )  1  (3 ,  2  , ) cot   ,  2 , ( ,   , )
                             rU 2r ,r  2 ,  r ,r cot  3 cot 
Hypersurface equations: hierarchical set for                              ,r ,U ,r and V,r
Evolution equation for                (r ),ur
 Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
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:                                           Null code test: time of bounce
• Bondi news function (from the metric variables expansion
at scri)
• Approximate gravitational waves (Winicour 1983, 1984,
1987):
      • Quadrupole news                                                                          150 x
                                                                                  Grid 1 : r 
                                                                                                 1 x4
      • First moment of momentum formula                                                                  
                                                                                  Grid 2 : r  100 tan    x
                                                                                                      2    

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
            Gravitational waves: consistency & disagreement
                                                  Good agreement in 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             A possible explanation: different velocities involved in both scenarios, 10-5-
                       rescaled by a factor        10-4c for a pulsating NS and 0.2c in core collapse.
                       50.
                                                   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
                              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 b  observer comoving with the fluid.
                     the

T   TPF  TEM
             
                                                                    T   hu  u  p g   b  b
  
TPF  hu  u  pg 
                                                                                         with the definitions:
                  1                         1     
TEM      F  F  g  F  F   u  u  g  b 2  b  b                             b 2  b b
                    4                         2     
                                                                                                      b2
                                    F     u b                                       p  p
                                                                                                       2
Ideal MHD condition:                F  u  0
                                                                                                     b2
electric four-corrent must
                                                                                               h  h
be finite.             
                         J  qu   F  u                                                           

       Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                             3+1 General Relativistic (Ideal)
                           Magnetohydrodynamics equations (2)
1.     Conservation of mass:             ( u  )  0                                             Adding all up: first-
                                                                                                    order, flux-
2.     Conservation of energy and momentum:                T   0                               conservative,
                                                                                                    hyperbolic system of
                                                  0 , F   u  B  u B  
                                                                1
3.     Maxwell’s equations:            F                                                       balance laws
                                                                W  
       •     Induction equation:
                                            1 
                                             t
                                                      
                                                                     
                                                     B    v    B
                                                                                                    + constraint
                                                                                                    (divergence-free
                                                                                                  condition)
       •     Divergence-free constraint:            B  0
                        i
     1  u
         
                   g f  
                            s
                                                                             Bi
                                                                                 0
      g  t       x                                                     x
                       i                                                       i
                           
                                                                            0              
        D                         ~i                                                
                                                                      gj
                                  Dv
                                                           T     gj  
       S j  i          2 ~i
                        h W v j v  p  j  b b j
                                         i    i
                                                                     x              
     u   f   2 i                                     s                         
         
         k       h W v  p   i /   Dv i  b 0bi 
                         ~                   ~
                                                         
                                                                   0  ln 
                                                                T                  0 
                                                                                  T   
                                                                                            
        B                ~          ~                               x 
                                                                                          
                         v i B k  v k Bi                                 k            
                                                                           0               
      Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                          Solution procedure of the GRMHD equations
                           i
         1  u     g f                                                                       Bi
                              s                                                                       0
             t
         g           x i                                                                       x i
                               
                                                                               Constrained transport scheme (Evans & Hawley 1988,
 • Same HRSC schemes as for GRHD equations                                     Tóth 2000).
  (HLL, Kurganov-Tadmor, Roe-type)
 • Wave structure information obtained                                         Field components defined at cell interfaces. Zone-
 • Primitive variable recovery more involved                                   centered vector (needed for primitive recovery and cell
                                                                               reconstruction & Riemann problem) obtained from
Details: Antón, Zanotti, Miralles, Martí, Ibáñez, Font & Pons,
          in preparation (2005)                                                staggered field components:

Update of field components:                                                                  Bi , j 
                                                                                                x       1
                                                                                                        2
                                                                                                          
                                                                                                          Bix j  Bix 1, j
                                                                                                            ,               
                                t
B x n 1
   i, j        
              Bix j
                 ,        
                          n
                              
                                y
                                   i , j 1  i , j 
                                t
B y n 1
   i, j        
              Biy j
                 ,        
                          n
                              
                                x
                                   i 1, j  i , j 
i , j 
        4
              
        1 ˆy x           ˆ
                                ,     ˆ          ˆ
           f Bi 1, j  f y Bix j  f x Biy j 1  f x Biy j
                                             ,           ,                      
                                                      
 These equations conserve the discretization of       B
                              Bix 1, j  Bix j          Biy j 1  Biy j
                  
                              
           B                                     
                                           ,              ,          ,

                                      xi , j               yi , j
                   i, j


       Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                       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                              HLL solver
Open circles: nonvanishing lapse function (2), at time t=0.2                                  1600 zones
Open squares: nonvanishing shift vector (0.4), at time t=0.16                                 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
                       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
                                                                                         
                                                                                              
                                                                                b  b , b ,0,of the
                                                                                      on top 0
                                                                                                  t   r
                                                                                                          
                                    hydrodynamic (Michel) accretion solution.           Radial magnetic
        internal energy
                                    field component chosen to satisfy divergence-free condition, and its
                                    strength is parametrized by the ratio:
                                                                   2p
                                                              
                                                                   b2
 radial velocity                    HLL solver
                                    100 zones
                                    =1

                                    Solid lines: analytic solution
                                    Circles: numerical solution
    radial magnetic field           (t=350M)

                                    Increasing the grid resolution
                                    shows that code is second-order
                                    convergent irrespective of the
                                    value of 

Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                       GRMHD equations: code tests (3)
                                     Magnetized equatorial Kerr accretion (Takahashi
   density                           et al 1990, Gammie 1999)
                                     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
                  azimuthal          mass flux FM, the angular momentum flux FL, the energy flux FE, and the
                  velocity           component F of the electromagnetic tensor.


                                     a=0.5
                                     FM=-1.0
                                     FL=-2.815344
  radial magnetic
                                     FE=-0.908382
  field
                                     F=0.5

                                     HLL solver

             azimuthal
             magnetic field          Solid lines: analytic
                                     solution
                                     Circles: numerical
                                     solution (t=200M)
                                                                             second order convergence

Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                   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                  The amplification factor of the initial magnetic field
precision during the simulation.                                    during the collapse is 1370.

     Institute for Pure and Applied Mathematics, University of California, Los Angeles, May 2005
                                       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

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:7
posted:7/31/2011
language:English
pages:49