GR-hydro

1
José A. Font
Departamento de Astronomía y Astrofísica
Universidad de Valencia (Spain)
Work done in collaboration with Pablo Cerdá-Durán
(UVEG)
2
Outline of the talk

• Introduction and motivation
• Framework
• General relativistic magneto-hydrodynamics
• Test magnetic field approximation
• Gravitational field (Einstein equations in the
  CFC approximation)
• Numerical approach
• Relativistic magneto-rotational core collapse
  simulations
• Initial data
• Code tests
• Collapse dynamics and B-amplification mechanisms
• Gravitational waveforms
• Magneto-rotational instability
• Long-term evolution of the B-field
• Summary and outlook

3
Introduction and motivation
The study of gravitational stellar collapse has
traditionally been one of the primary problems in
relativistic astrophysics. It is a distinctive
example of a research field in astrophysics where
essential progress has been accomplished through
numerical modeling with increasing levels of
complexity in the input physics/mathematics.
Basic set of equations to solve Hydrodynamics
Gravitational field (
B-fields transport ). The GR-Hydro
equations (as well as the GRMHD counterpart)
constitute a nonlinear hyperbolic system. There
exist solid mathematical foundations and accurate
numerical methodology imported from CFD. In
recent years, the so-called (upwind)
high-resolution shock-capturing schemes (or
Godunov-type methods), based upon approximate
Riemann solvers, have been successfully extended
from classical to relativistic fluid dynamics
(see e.g. Martí Müller 2003 and Font 2003
articles at Living Reviews in Relativity). While
such advances also hold true in the case of the
MHD equations, the development still awaits here
for a thorough numerical exploration. But sure
enough, numerical GRMHD is a hot topic these days

4
Motivation intense work in recent years on
formulating/solving the MHD equations in general
relativistic spacetimes (either background or
dynamical). Pioneers Wilson (1975), Sloan
Smarr (1985), Evans Hawley (1988), Yokosawa
(1993) More recently Koide et al (1998 ), De
Villiers Hawley (2003 ), Baumgarte Shapiro
(2003), Gammie et al (2003), Komissarov (2005),
Duez et al (2005 ), Shibata Sekiguchi (2005),
Antón et al (2006), Neilsen et al (2006). Both,
artificial viscosity and HRSC schemes developed.
Most of the applications are in the field of
black hole accretion and jet formation
Development of the MRI in a magnetised torus
around a Kerr black hole (Gammie, McKinney
Tóth 2003)
Jet formation the twisting of magnetic field
lines around a Kerr black hole. The yellow
surface is the ergosphere (Koide et al 2002)
many others under way (you name it!)
The number of groups working in special
relativistic MHD is even larger Komissarov
Balsara Koldoba et al Del Zanna et al Leisman
et al Exact solution of the SRMHD Riemann
problem found recently Romero et al (2005)
particular case Giacomazzo Rezzolla (2005)
general case.
5
State-of-the-art of core collapse simulations
(table courtesy of Harry Dimmelmeier)
Talks by Dimmelmeier and Ott will cover part of
the ongoing work.
Neutron stars have intense magnetic fields
(1012-1013 G) or even larger at birth (1014-1015
G), inferred from studies of anomalous X-ray
pulsars and soft gamma-ray repeaters (Kouveliotou
et al 1998). In magnetars (extreme case) the
B-field can be so strong to affect the internal
structure of the star (Bocquet et al
1995). Performing GR-MHD core collapse
simulations is the next natural step
6
Status of magneto-rotational core collapse (in a
nutshell)

• Proposed as a possible supernova explosion
  mechanism by Müller Hillebrandt (1979)
  magnetohydrodynamical shock.
• Initial conditions for magneto-rotational core
  collapse (Heger et al 2005) initial rotation
  rates are more than one order of magnitude
  smaller than
• those used in past parameter studies of
  magneto-rotational core collapse (e.g.
  Dimmelmeier, Font Müller 2002)
• those predicted by previous evolutionary
  calculations.
• Strength and distribution of the initial
  magnetic field in the core unknown. If initially
  weak, there exist several amplification
  mechanisms
• differential rotation (?-dynamo transforms
  rotational energy into magnetic energy, winding
  up any seed poloidal field into toroidal field).
• MRI, exponential growth of the field strength.
  It occurs if the radial gradient of angular
  velocity is negative, which holds in core
  collapse simulations.
• Magnetic fields (torques) can affect the
  collapse dynamics in a major way by
  redistributing angular momentum (Meier et al
  1976 Spruit Phinney 1998 Spruit 2002 Heger
  et al 2004), which can slow down the forming PNS.
  Slowly-rotating neutron stars at birth (10-15
  ms).
• Magneto-rotational core collapse may lead to
  jet-like explosions (distinctive GW signature
  Obergaulinger et al 2006).
• No relativistic magneto-rotational core collapse
  simulations yet available.

7

• Ideal MHD simulations, 2D, Newtonian physics
  microphysics very active field in recent years
  Wheeler et al 2003 Akiyama et al 2003, Kotake et
  al 2004a,b, 2005 Takiwaki et al 2004 Wheeler
  Akiyawa 2004 Yamada Sawai 2004 Ardeljan et al
  2005 Sawai et al 2005.
• Magneto-rotational effects on the GW signature
  from core collapse were first studied in detail
  by Kotake et al (2004) and Yamada Sawai (2004).
  Found differences from the signature of purely
  hydrodynamical models for very strong initial
  fields (B1012 G)

• Obergaulinger, Aloy Müller (2006) most
  comprehensive parameter study up-to-date of
  magneto-rotational core collapse.
• rotating polytropes, 2D, Newtonian (modified
  Newtonian)
• SQF to compute gravitational waves
• ad-hoc initial poloidal B fields (no
  consistent solution known)
• Weak initial fields (B1011 G)
• change neither collapse dynamics nor resulting
  GW signal (most relevant case, astrophysics-wise).
  
• Strong initial fields (B1011 G)
• slow down core efficiently (even retrograde
  rotation occurs)
• cause qualitatively different dynamics and GW
  signals
• highly bipolar, jet-like outflows occur

8
Core collapse and gravitational waves
Numerical simulations of stellar core collapse
are nowadays highly motivated by the prospects of
direct detection of the gravitational waves
emitted.
International network of resonant bar detectors
International network of interferometers
Do GRMHD effects modify the existing GW signals
from hydrodynamical rotational core collapse?
9
Early core collapse simulations GR hydrodynamics
Equations of motion local conservation laws of
density current and stress-energy

• Conservative formulations well-adapted to
  numerical methodology
• Banyuls et al (1997) Font et al (2000) 31,
  general EOS

First-order flux-conservative
hyperbolic system
Solved using HRSC schemes (either upwind or
central)
10
Numerical schemes HRSC
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 procedures of 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
High-order symmetric (central) schemes also
available (Kurganov Tadmor 2000)
11
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, rarefaction waves,
  vortices, etc

V0.99999c (W224)
Wind accretion onto a Kerr black hole
(a0.999M) Font et al, MNRAS, 305, 920 (1999)
Simulation of a extragalactic relativistic
jet Scheck et al, MNRAS, 331, 615 (2002)
12
CFC metric equations
In the CFC approximation
(Isenberg 1985 Wilson Mathews 1996) the ADM
31 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, ?ukS(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.
13
Relativistic non-magnetized core collapse

• HRSC scheme PPM Marquina flux-formula
• Solid line (CFC) relativistic simulation
• Dashed line Newtonian
• Larger central density in relativistic models
  (more compact PNS)
• Similar gravitational radiation amplitudes
• Multiple bounce collapse suppressed in GR

Axisymmetric simulations Dimmelmeier, Font
Müller (2002) CFC Cerdá-Durán et al (2005)
CFC Shibata Sekiguchi (2005) BSSN
Excellent agreement Extensions to realistic EoS
and 3D available (see talks by Dimmelmeier and
Ott)
14
CFC versus full General Relativity
CFC is indeed an excellent approximation to GR
for studying core collapse
Dimmelmeier et al (astro-ph/0603760)
N Newtonian potential R TOV potential A
Modified TOV potential
15
General relativistic magneto-hydrodynamics (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-current must
be finite.
16
General relativistic magneto-hydrodynamics (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
• Maxwells equations
• Induction equation
• Divergence-free constraint

Antón et al. (2006)
17
Notation and Riemann problem
GRMHD equations are rewritten in terms of the
conserved quantities inside a coordinate volume
flat Nabla operator
The associated flux-vector Jacobians in every
direction are 7x7 matrices. The solution of the
eigenvalue problem (Antón 2006) leads to seven
types of waves Entropic wave Alfvén
waves Magnetosonic waves (no analytic
expression)
Characteristics of the GRMHD equations
18
Test magnetic field approximation
Test B-field approximation to study the collapse
of weakly magnetized stellar cores, it is a good
approximation to consider that the magnetic field
entering in the energy-momentum tensor is
negligible when compared with the fluid part.
Therefore, the magnetic field evolution
(induction equation) does not affect the dynamics
of the fluid, which is governed solely by the
hydrodynamics equations. However, the evolution
of B does depend on the fluid evolution, due to
the presence of the velocity field in the
induction equation.
Practical numerical advantage the seven
eigenvalues of the GRMHD Riemann problem reduce
to three (purely hydro)
(see Antón et al (2006) for complete expressions)
19
Solution procedure for the GRMHD equations
Constrained transport scheme (Evans Hawley
1988, Tóth 2000) Flux-CT scheme

• Same HRSC schemes as for GRHD equations
  (HLL, Kurganov-Tadmor, Roe-type)
• Wave structure information available
• Primitive variable recovery more involved
• See Noble et al 2006 for methods comparison
• Antón et al 2006 find the roots of an 8-th order
  polynomial using a 2d-Newton-Raphson method.

In the test B-field approximation the primitive
recovery is as in the purely GRHD case (as there
is no B-field components in the momentum and
energy equations)
20
Magnetic field evolution flux-CT (1)
The B-field evolution is given by the induction
equation and the divergence constraint. HRSC
schemes used for the hydrodynamics do not
preserve the divergence-free condition. An
ad-hoc scheme has to be used to update the
magnetic field components. Main physical
implication of divergence constraint is that the
magnetic flux through a closed surface is zero
essential to the constrained transport (CT)
scheme (Evans Hawley 1988, Tóth 2000).
For any given surface, the time variation of the
magnetic flux across the surface is
Induction equation
Stokes theorem
The magnetic flux through a surface can be
calculated as the line integral of the electric
field along its boundary.
21
Magnetic field evolution flux-CT (2)
Numerical implementation in axisymmetry
Assumption B-field components constant at each
cell surface E-field
components constant along each cell edge
Evolution equations for the B-field (CT scheme)
(axisymmetry condition)
The polodial (r and ?) B-field components are
defined at cell interfaces (staggered grid)
The total magnetic flux through the cell
interfaces is given by
If the initial data satisfy the divergence
constraint, it will be preserved during the
evolution
22
Magnetic field evolution flux-CT (3)
Discretisation
Equations used by the code to update the
B-field. The only remaining aspect is an
explicit expression for the E-field.
The E-field components can be calculated from the
numerical fluxes of the conservation equations
for the B-field. Done solving Riemann problems at
cell interfaces (characteristic information of
the flux-vector Jacobians incorporated in the
B-field evolution). This procedure is only
valid for r and ? E-field components.
Balsara Spicer (1999) proposed a practical way
to compute the ? component of the E-field from
the numerical fluxes in adjacent interfaces.
Resulting scheme flux-CT
23
Code tests in strong gravity (black holes)
Magnetised spherical accretion onto a
Schwarzschild BH
density
radial velocity
Test difficulty keep stationarity of the
solution. Used in the literature (Gammie et al
2003, De Villiers Hawley 2003)
radial magnetic field
internal energy
2nd order convergence
density
radial magnetic field
Magnetised equatorial Kerr accretion (Takahashi
et al 1990, Gammie 1999)
Test difficulty keep stationarity of the
solution (algebraic complexity augmented, Kerr
metric) Used in the literature (Gammie et al
2003, De Villiers Hawley 2003)
azimuthal velocity
azimuthal magnetic field
24
Initial data
We work with a vector potential such that the
associated B-field is divergence-free.
Given a vector potential, the B-field components
at cell interfaces can be discretized such that
the numerical values of the magnetic flux over
cells are zero up to round-off error. This value
is preserved during the evolution using the
flux-CT scheme.
1) Homogeneous starred B-field
3) Toroidal starred B-field
magnetic field at r0
2) B-field generated by a circular current loop
25
Toroidal test
Let us consider a rotating axisymmetric
configuration with no meridional flows
Induction equation
?-dynamo amplification mechanism (Meier 1976)
Equilibrium solution can be found in three cases
Purely toroidal B-field
Rigid rotator
Special case
Useful for code validation!
Toroidal test A circular current loop non
evolving rigidly rotating fluid of constant
density
Global error in the toroidal B-field
minmod
Grid resolution f1 80x10 f2 160x20 f4
320x40 Convergence order minmod 2.35 MC
2.16
TTB
minmod º MC
MC
TTA
minmod x MC
PPM
Local convergence order
26
Poloidal test
Let us now consider a fluid with only radial
velocities and a purely poloidal magnetic field
From the induction equation and the continuity
equation, it can be shown that the following
equivalence holds in the equatorial plane
We define a Lagrangian coordinate system and
check whether does not
change with time (as a function of the mass
enclosed in a given radius)
Spherical collapse of a 4/3 polytrope
time evolution of the error during the infall
phase for different reconstructions and grid
resolutions
t35 ms Dotted curve initial profile
travelling shock
lt 1
minmod KT solver
PNS boundary
27
Summary of initial models

• Hydrodynamical models follow the notation of
  Dimmelmeier, Font Müller (2002).
• Magnetised models follow notation of
  Obergaulinger, Aloy Müller (2006).
• Magnetic field configuration
• Poloidal field generated by a current loop (CL)
• Toroidal field (T)
• Initial magnetic field strength 1010 G (since we
  use the test B-field approximation all results
  can be rescaled to other initial field values)
  and rmag400 km.
• Last column indicates how stable is initially the
  B-field configuration.
• We focus our discussion on models A1B3G5
  (representative of the sample Cerdá-Durán
  Font, in preparation).

28
A1B3G5-D3M0 dynamics energy parameters
Their evolution allows to quantify the signatures
of rotation and B-field on the collapse dynamics.

• B-field amplification mechanisms present in our
  axisymmetric simulations
• radial compression radial velocity compress
  B-field lines perpendicular to the velocity only
  during infall phase acts on poloidal and
  toroidal components.
• ?-dynamo winds up poloidal field lines into
  toroidal field lines. The latter are generated
  and amplified even from purely poloidal initial
  configurations. The toroidal B-field dominates
  after core bounce.
• ?-?-dynamo (Spruit 1999), toroidal ? poloidal
  NOT possible (3D effect).
• MRI (Balbus Hawley 1991) NOT possible (test
  B-field approximation)

total
toroidal
poloidal
RC ? ?
The final ?mag 1 (and ?rot) as the B-field
strength is small enough not to affect the
dynamics
29
A1B3G5-D3M0 dynamics morphology
density velocity
poloidal B-field lines
toroidal B-field orientation
time of bounce (t30 ms)
At bounce the initial poloidal configuration is
highly distorted by the radial compression of the
field lines. Maximum values of the poloidal
B-field reached at the center. A toroidal
component has formed surrounding the inner
core. The travelling shock wave is visible in
both the hydrodynamical variables and the
magnetic field.
30
A1B3G5-D3M0 dynamics morphology
density velocity
poloidal B-field lines
toroidal B-field orientation
final time (t60 ms)
At the final time a compact toroidal B-field is
confined in the outer layers of the PNS the
inner region is dominated by the poloidal
B-field. Convective motions (meridional flows)
responsible of the twisting of the poloidal
field, which changes direction.
changes sign
The ?-dynamo mechanism changes the sign of the
toroidal B-field in a shell (B?gt0) surrounded by
toroidal field in the opposite direction (B?lt0).
Such configuration is maintained until the end of
simulation, resulting in a linear amplification
of the toroidal B-field.
31
A1B3G5-D3M0 dynamics field lines
Equations of the magnetic field lines
bounce
final time
The box represents the region above the
equatorial plane. Axes are in km. Colours are
proportional to the magnetic flux in each line
(bluelow, green-redhigh).

• Outside the PNS the B-field is roughly poloidal.
• A shell of twisted toroidal field lines formed
  around the inner core of the PNS.
• Poloidal structure of innermost PNS region
  hidden by the shell.
• Longer time simulations required to clarify the
  stability of such configuration.

32
A1B3G5-T3M0 dynamics B-field morphology

• Purely toroidal B-fiel initially. Since no
  poloidal field present initially, toroidal
  B-field can not be amplified via ?-dynamo.
• ?mag decreases during evolution as the only
  amplification mechanism acting is RC (much slower
  than ?-dynamo). Final PNS less magnetised than
  progenitor.
• Toroidal B-field does not generate poloidal
  B-field due to axisymmetry (no ?-?-dynamo).
  Drawback of the simulation.
• Mixed model (DT3M0) long-term evolution same as
  D3M0 model. The (initial) toroidal B-field does
  not amplify by RC nor transforms into poloidal,
  while the (initial) polodial B-field gets
  amplified by RC and transformed into toroidal.
• If ?-?-dynamo were not important, the B-field
  distribution in the PNS (toroidal) would depend
  almost exclusively on the poloidal B-field
  distribution of the collapse progenitor (and
  barely on the toroidal one).

Black T3M0 Red D3M0 Blue DT3M0
Toroidal B-field PNS
33
A4B5G5-D3M0 extreme case
Faster initial rotation. Differential rotation.
?-dynamo stronger.
Maximum density reached off-center (strong
centrifugal rotation). B-field weak in the inner
hole, since no matter flows transport field
lines into this region. Strongest B-field
reached outside shock location at bounce, the
shock dragging the B-field lines along its path.
Complicated windings and twists of the field
lines by the action of ?-dynamo and meridional
flows. Positive toroidal B-field regions
formed. B-field shows large variations in small
scales finite conductivity effects such as field
line reconnection could be important (magnetic
diffusion proportional to ?B).
34
PNS B-field evolution under ?-dynamo
Fits of ?mag to obtain B-field growth and
saturation timescales.
?rot
A1B3G3-D3M0
?mag
We can estimate the B-field saturation timescale
in the newly-formed PNS under the (unlikely)
assumption that ?-dynamo were the only
amplification mechanism at work (transforming
rotational energy into magnetic energy). During
the quasi-stationary evolution of the PNS, and
assuming fixed angular profiles and poloidal
B-field distribution, the toroidal B-field grows
linearly with time.
As the PNS rotates slower the B-field saturates
because ?-dynamo cannot extract more energy from
rotation complete halting of the PNS, exchange
sense of rotation, or torsion pendulum. However,
MRI is likely to amplify the B-field to
saturation on much shorter timescales ( in a few
slides)
35
Comparison with existing simulations

• Comparisons made with models by Obergaulinger et
  al (2006a,b)
• D3M10 models chosen from their sample, which
  have low magnetisation and hence our passive
  field approximation applies. Obergaulinger et al
  (2006) checked that the dynamics does nor change
  (MHD versus hydro) for such models.
• B-field dynamics and morphology qualitatively
  similar in CFC models than in Newtonian gravity
  and TOV models (as long as the hydrodynamics is
  similar in both cases, which does not apply to
  Newtonian, multiple bounce core collapse models).
• Magnetisation (?mag) of final PNS smaller in CFC
  than in both Newton and TOV. Irrespective of
  numerical scheme (solver, minmod, PHM) and
  resolution. General trend?

36
Gravitational waves
A1B3G5
A4B5G5
A1B3G3
Black hydro GW Red magnetic GW,
D3M0 Blue magnetic GW, T3M0 Green as blue
but counter-rotating
Hydro stress formula
MHD
The magnetic component of the GW is several
orders of magnitude smaller than the
hydrodynamical component. However As the
B-field amplification reaches saturation, the
effects of the field on the waveforms are
expected to be significant. In particular MRI may
noticeably change the waveforms, as it is the
most efficient amplification mechanism.
37
Magneto-rotational instability (1)
MRI is a shear instability that generates
turbulence and an amplification of the magnetic
field in rotating magnetized plasma, transporting
angular momentum in the star (Balbus Hawley
1991). Magnetized collapse models are indeed
susceptible of developing MRI (Akiyama et al
2003 Obergaulinger et al 2005). The (Newtonian)
condition for the MRI to occur (neglecting
buoyancy effects and if the B-field strength is
very low) is
If this condition is fulfilled and the magnetic
field has a poloidal component, MRI grows
exponentially in time.
The timescale of the fastest growing unstable
mode is
independent of B-field configuration and strength
A1B3G3-D3M0
Main caveat of our simulations no back-reaction
of B-field onto the dynamics. Estimate only the
potential effects of MRI.
38
Magneto-rotational instability (2)
A1B3G3-D3M0
A1B3G5-D3M0
A4B5G5-D3M0
time of bounce
end of simulation
Depicted in white are the regions where MRI is
not fulfilled or timescale exceeds 1s. MRI does
not affect the collapse dynamics (infall phase)
but the evolution after core bounce is dominated
by this instability. Passive B-field
approximation no longer valid. Models T3M0 are
not affected by MRI (no poloidal B-field).
However, relaxing axisymmetry condition leads to
poloidal B-fields via ?-?-dynamo, and hence to
MRI.
39
Summary of the talk

• Multidimensional simulations of relativistic
  magneto-rotational core collapse feasible
  nowadays with current formulations of GRMHD and
  Einsteins equations.
• We have designed a method to build weakly
  magnetised stars in GR with either toroidal or
  poloidal (or both) B-field components.
• Tests performed to check accuracy and
  convergence of our GRMHD code.
• Results from CFC relativistic simulations of
  magneto-rotational core collapse to NS in
  axisymmetry presented. Passive B-field
  approximation used.
• Core collapse dynamics, B-field distribution,
  and gravitational waveforms analysed. No major
  differences found with previous Newtonian
  studies. Strength of final B-field smaller in CFC
  than in Newtonian case.
• B-field amplification mechanisms investigated
  (radial compression and ?-dynamo).
• ?-?-dynamo and MRI not included in the
  modelling. Growth time of MRI shows it needs to
  be included as it dominates the post-bounce
  dynamics within a few ms.
• Current GRMHD code already extended to relax
  passive B-field approximation. Tests in progress.