Evolution of the Magnetized, Neutrino-Cooled Accretion Disk in the Aftermath of a Black Hole Neutron Star Binary Merger

Evolution of the Magnetized, Neutrino-Cooled Accretion Disk in the Aftermath of a Black Hole Neutron Star Binary Merger

Fatemeh Hossein Nouri Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164, USA Inter-University Centre for Astronomy and Astrophysics (IUCAA), Post Bag 4, Ganeshkhind, Pune 411 007, India    Matthew D. Duez Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164, USA    Francois Foucart Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA 94720, USA Department of Physics, University of New Hampshire, Durham, New Hampshire 03824, USA    M. Brett Deaton Joint Institute for Nuclear Astrophysics, Michigan State University, East Lansing, MI 48824, USA Department of Physics, North Carolina State University, Raleigh, NC 27695, USA Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164, USA    Roland Haas NCSA, University of Illinois at Urbana-Champaign, Urbana, Illinois, 61801, USA Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam-Golm, 14476, Germany TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA    Milad Haddadi Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164, USA    Lawrence E. Kidder Center for Radiophysics and Space Research, Cornell University, Ithaca, New York, 14853, USA    Christian D. Ott TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA Center for Gravitational Physics and International Research Unit of Advanced Future Studies, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, Kyoto Prefecture 606-8317, Japan    Harald P. Pfeiffer Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, Ontario M5S 3H8, Canada Canadian Institute for Advanced Research, 180 Dundas St. West, Toronto, ON M5G 1Z8, Canada Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam-Golm, 14476, Germany    Mark A. Scheel TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA    Bela Szilagyi TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA

Black hole-torus systems from compact binary mergers are possible engines for gamma-ray bursts (GRBs). During the early evolution of the post-merger remnant, the state of the torus is determined by a combination of neutrino cooling and magnetically-driven heating processes, so realistic models must include both effects. In this paper, we study the post-merger evolution of a magnetized black hole-neutron star binary system using the Spectral Einstein Code (SpEC) from an initial post-merger state provided by previous numerical relativity simulations. We use a finite-temperature nuclear equation of state and incorporate neutrino effects in a leakage approximation. To achieve the needed accuracy, we introduce improvements to SpEC’s implementation of general-relativistic magnetohydrodynamics (MHD), including the use of cubed-sphere multipatch grids and an improved method for dealing with supersonic accretion flows where primitive variable recovery is difficult. We find that a seed magnetic field triggers a sustained source of heating, but its thermal effects are largely cancelled by the accretion and spreading of the torus from MHD-related angular momentum transport. The neutrino luminosity peaks at the start of the simulation, and then drops significantly over the first 20 ms but in roughly the same way for magnetized and nonmagnetized disks. The heating rate and disk’s luminosity decrease much more slowly thereafter. These features of the evolution are insensitive to grid structure and resolution, formulation of the MHD equations, and seed field strength, although turbulent effects are not fully converged.


I Introduction

The cause of short-hard gamma ray bursts (GRBs) remains unknown, but some of the most promising central engine models involve rapid ( s) accretion onto a stellar mass black hole (BH). Such systems are naturally produced by some black hole-neutron star (BHNS) and neutron star-neutron star (NSNS) binary mergers. (For reviews of short GRBs, see Nakar (2007); Berger (2014).)

Given the requisite dense, hot accretion flow, there are several ways energy could be channeled into a baryon-poor ultra-relativistic outflow of the sort needed to explain GRB properties. The accretion gas cools primarily by neutrino emission, and so such systems are classified as neutrino-dominated accretion flows (NDAFs) Popham et al. (1999); Di Matteo et al. (2002); Janiuk et al. (2004). Some emitted neutrino energy can be transferred to a pair fireball through neutrino-antineutrino annihilations outside the disk Paczynski (1991); Jaroszynski (1996); Richers et al. (2015); Just et al. (2016). Magnetic fields can also extract energy from the disk or black hole spin Blandford and Znajek (1977); Meszaros and Rees (1997), and the energy outflow can be Poynting flux dominated.

The lifetime of a short GRB (, presumably related to the disk lifetime ) is much greater than the dynamical timescale ( ms) and perhaps also the thermal timescale ( in the standard alpha viscosity, thin disk model Shakura and Sunyaev (1973)). Therefore, the GRB mechanism is a process that takes place in the accretion system’s dynamical and probably also thermal equilibrium.

The post-merger accretion disks formed in BHNS/NSNS mergers have densities of and temperatures of . Hence, photons are trapped and in equilibrium, and radiation is by neutrinos. For high enough accretion rate , the disk is opaque to neutrinos, which must diffuse out and provide an additional source of pressure. Neutrino luminosities can reach erg s, and this emission will strongly affect the disk (on a secular timescale ) by cooling it and altering the composition, quantified by the electron fraction , the fraction of nucleons that are protons. Unstable entropy or gradients can drive convection in the disk Lee et al. (2005). In addition to these emission effects, there are also neutrino transport effects. Neutrino absorption near the neutrinosphere can drive thermal winds Dessart et al. (2009); Perego et al. (2014); neutrino momentum transport can create a viscosity that slows the growth of the magnetorotational instability Masada and Sano (2008); Guilet et al. (2015) (although probably not for BHNS mergers Foucart et al. (2015); Kiuchi et al. (2015)).

In previous papers Deaton et al. (2013); Foucart et al. (2014, 2015), we simulated BHNS mergers at realistic mass ratios using a finite-temperature nuclear equation of state and incorporating neutrino effects. The latter were modeled in some cases with a leakage scheme (which includes emission but not transport) Ruffert et al. (1996); Rosswog and Liebendörfer (2003); O’Connor and Ott (2010); Deaton et al. (2013); Foucart et al. (2014) and in some cases with an energy-integrated two-moment M1 transport scheme Shibata et al. (2011); Foucart et al. (2015). Comparing to the earlier times of evolution we found that the post-merger accretion disks become cold, and more neutron-rich with dimmer neutrino emission after a few tens of milliseconds. Comparing leakage to M1, we find that the former gives a reasonable estimate for the neutrino emission and global thermal evolution, although it overestimates temperature gradients, and cannot accurately track the evolution in low-density regions. No significant neutrino-driven winds were seen. The cooling and dimming of the disks is unsurprising, given that these simulations included the major cooling mechanisms–neutrino emission and advection of the hot inner gas into the black hole–but contained only one significant heating mechanism (in addition to numerical dissipation heating): shock heating from the circularization and pulsation of the disk gas.

Long-term accretion requires an angular momentum transport process that will naturally release orbital energy and heat the gas. This is probably provided by magnetic fields, which were not included in the above simulations. Weakly magnetized accretion flows are subject to the magnetorotational instability (MRI) Balbus and Hawley (1998), inducing turbulence which dissipates energy at small scales and whose mean (mostly Maxwell) stresses transport angular momentum outward, driving accretion Hawley et al. (1995). Magnetic fields also transport angular momentum through magnetic winding (the effect). Reconnection at current sheets provides a way to convert magnetic energy into plasma thermal and kinetic energy. Simulations of radiatively inefficient magnetized accretion tori find strong winds along disk surfaces and magnetically dominated poles De Villiers and Hawley (2003); McKinney and Gammie (2004). Large-scale fields threading the BH ergosphere enables extraction of the black hole spin energy into a Poynting flux-dominated jet Blandford and Znajek (1977); McKinney and Gammie (2004).

There have been successful GRMHD simulations, neglecting neutrino effects, of BHNS Chawla et al. (2010); Etienne et al. (2012, 2013); Paschalidis et al. (2015); Kiuchi et al. (2015) and NSNS Liu et al. (2008); Anderson et al. (2008); Giacomazzo et al. (2011); Rezzolla et al. (2011); Kiuchi et al. (2014, 2015) mergers. The highest resolution BHNS simulations with an initial seed field confined in the neutron star Kiuchi et al. (2015) find strong winds and Poynting-dominated jets only at very high resolutions (and even here, it is unclear that convergence has been achieved). There are also indications that unconfined seed fields produce jets more readily Paschalidis et al. (2015), consistent with disk studies that find jets but not disk interiors to be very sensitive to the seed field Beckwith et al. (2008). The helicity of the magnetic field may also have subtle long-term effects Wan (2017). These merger simulations used simplified thermal components of the equation of state and neglected neutrino effects; they had the main heating effects but not a major cooling effect.

Clearly, accurate evolution on thermal timescales requires both neutrino cooling and magnetoturbulent heating. The two will influence each other. The neutrino luminosity, and hence the viability of “neutrino” mechanisms for driving a GRB, depends on magnetic heating, while the saturation strength of the magnetic field in an MRI turbulent disk will depend on the temperature of the gas Sano et al. (2004); Shi et al. (2016) set partly by neutrino cooling. NSNS merger simulations with both effects have been performed Neilsen et al. (2014); Palenzuela et al. (2015), but our understanding of long-term post-merger evolution of BHNS (and high-mass NSNS) systems relies on accretion disk models. In most cases, turbulent transport and dissipation is modeled by a phenomenological “alpha” viscosity Shakura and Sunyaev (1973). These include the original one-dimensional (axisymmetric, vertically summed), equilibrium NDAF studies Popham et al. (1999); Di Matteo et al. (2002); Janiuk et al. (2004). One-dimensional NDAFs were evolved by Janiuk et al. Janiuk et al. (2007), who found disks can become visco-thermally unstable in some regions, but only for very high accretion rates (). Evolutions have also been carried out in higher dimensions, again in the alpha viscosity framework, yielding valuable information on neutrino-antineutrino energy release and late-time outflows Setiawan et al. (2004); Lee et al. (2005); Fernández and Metzger (2013); Fernández et al. (2015). Efficient release of energy by radiation requires low (so ), proving Lee et al. (2005) the importance of first-principles, magnetohydrodynamic (MHD) simulations to assess the adequacy of viscosity models and to reveal the actual efficiency of angular momentum transport.

MHD disk simulations with neutrino cooling have been carried out in 2D beginning from analytic, constant angular momentum equilibrium tori by several groups Shibata et al. (2007); Barkov and Baushev (2011); Shibata and Sekiguchi (2012); Janiuk et al. (2013) , and recently in 3D by Siegel and Metzger Siegel and Metzger (2017). They identify the MRI, with associated heating, neutrino emission, and powerful outflows. These studies probably provide the most realistic picture available of the evolution of the post-merger disk, but their artificial disk profiles neglect the strong angular momentum gradients, high compactness, and nonaxisymmetric features seen in merger simulations. These neglected features will most likely have strong effects in the early, and most neutrino luminous, post-merger phase. In addition, the 2D (axisymmetric) simulations Shibata et al. (2007); Barkov and Baushev (2011); Shibata and Sekiguchi (2012); Janiuk et al. (2013) are affected by the known differences between the saturation of the MRI in 2D vs. 3D Goodman and Xu (1994); Obergaulinger et al. (2009), including the impossibility of an axisymmetric dynamo Hawley and Balbus (1992); Brandenburg (1996).

In this paper, we study the effects of magnetic fields on the post-merger evolution of a BHNS binary system. We evolve in 3D using as initial data the BH accretion flow produced by a BHNS merger simulation Foucart et al. (2014). In addition to MHD, we employ a realistic finite-temperature nuclear equation of state and neutrino cooling via a leakage approximation, giving us all the basic ingredients needed for a realistic thermal evolution. For this first study, we restrict ourselves to a simple seed field geometry with high field strength, for which the MRI is resolved with modest grid sizes. Studying a strongly magnetized disk most likely gives us a sense of the maximum effect that magnetic fields can have. Our simulations use the Spectral Einstein Code (SpEC) and required the development of new numerical techniques for SpEC: MHD on a cubed-sphere multipatch grid, coordinate maps to optimize grid use, and an improved technique to control entropy evolution in regions where kinetic energy dominates over internal energy.

Comparing disk evolutions with and without magnetic fields, we find some expected effects. The magnetic field drives strong and sustained accretion, while the late-time accretion rate of a nonmagnetized disk is, by comparison, negligible. Magnetic effects also do increase the disk’s specific entropy, as a result of magnetoturbulent heating and numerical reconnection, leading to a roughly steady entropy in comparison to the secularly decreasing entropy of a nonmagnetized disk. However, at early times the nonmagnetized disk’s cooling rate is significantly slower than neutrino emission would predict, indicating the continued importance of shock heating 30 ms after merger as a heating source of comparable strength to MHD-related heating. The effects of disk depletion and heating on the neutrino luminosity roughly cancel, and the magnetized disk dims at roughly the same rate as the nonmagnetized disk. Thus, for the case we consider, MHD turbulence does little to assist neutrino-related mechanisms for powering a GRB during the most neutrino luminous phase of the accretion, even in the case of an extremely strong seed field.

This paper is organized as follows. In Sec.II, the initial configuration and set up is discussed. Section III briefly describes the numerical methods used. In Sec.IV, numerical results are presented, focusing on the effects of magnetic field on the accretion rate, thermal evolution and general properties of the disk. Finally, Sec.V is devoted to the summary and conclusion. A detailed discussion of new numerical techniques is reserved for the Appendix.

Ii Initial state

ii.1 Input physics

As in our recent BHNS merger studies Deaton et al. (2013); Foucart et al. (2014) we employ the Lattimer-Swesty equation of state Lattimer and Swesty (1991) with nuclear incompressibility (LS220), using the table available at http://www.stellarcollapse.org and described in  O’Connor and Ott (2010).

Neutrino emission effects are captured using a simple leakage scheme, described in Deaton et al. (2013); Foucart et al. (2014). Leakage schemes remove energy and alter lepton number at rates based on the local free-emission and diffusion rates. They account for these emission effects within factors of accuracy (as determined by comparisons with genuine neutrino transport schemes Foucart et al. (2015)) but do not include the effects of neutrino transport and absorption. Our leakage scheme integrates out spectral information, assuming Fermi-Dirac distributions at the local temperature (with chemical potentials estimated as in Foucart et al. (2014)), although we can estimate an average energy of emitted neutrinos from the total luminosity and number emission rate. (See Perego et al. (2014, 2016) for approximate ways, not pursued in this study, to incorporate absorption and spectral information in a leakage framework.) Our leakage scheme includes -capture processes, pair annihilation, plasmon decay and nucleon-nucleon Bremsstrahlung interactions. In optically thick regions, the neutrinos contribute to the pressure.

ii.2 Initial Configuration

For our initial state, we use the BHNS configuration M12-7-S9 presented in Foucart et al. (2014). (See Table 2 of that paper.) The initial masses of the BH and NS are and respectively. The BH is rapidly spinning with . The remnant torus mass is about , with maximum density of g cm, and average temperature of MeV. We restart our simulation using data of this case at after merger. At this time, the spacetime has settled to a nearly stationary BH metric in the coordinate system produced by the numerical relativity simulation, but the disk remains significantly nonaxisymmetric and nonstationary. We therefore evolve only the fluid, keeping the metric at its initial state.

We set up an initially poloidal magnetic field via a toroidal vector potential


where is the axisymmetrized density field (to initiate the field with large poloidal loops), is the cylindrical radius in grid coordinates, sets the overall strength of the resulting -field, and the cutoff density , set to 6% of the maximum density, confines the initial field to regions of high-density matter. We follow the same prescription as that in Noble et al. Noble et al. (2009) to set the initial magnetic field strength, so that the ratio of the volume-weighted integrated gas pressure to the volume-weighted integrated magnetic pressure is about 13 for our strongly magnetized disk. This magnetic field at the maximum value is about . This is likely much stronger than realistic BHNS post-merger magnetic fields. We focus on this extreme case first for two reasons. First, it allows us to resolve the rapidly-growing modes of the MRI very well with modest resolution. Second, an extreme field might be expected to reveal the maximum effect that magnetic fields might have.

Strong seed fields may induce qualitatively different behavior from weaker seeds if it is strong enough to suppress the MRI before the disk can become turbulent. This will certainly be the case where is near or below unity, so the fastest-growing MRI mode wavelength () exceeds the disk height. This is not a danger in most of our disk. However, MRI growth might also be affected if is comparable to the length scale on which itself varies (due to variation in Alfven speed) Piontek and Ostriker (2004); Pino and Mahajan (2008); Obergaulinger et al. (2009) or comparable to the radius of curvature of the field lines, which occurs even in some strong-field, high-density regions. In order to estimate the effect of seed field strength, we carry out another simulation with a weaker seed field, set by . This corresponds to a maximum field strength of  G. This simulation does show weaker heating and less outflow, confirming our expectation (also suppported by 2D strong-seed disk simulations Janiuk et al. (2013)) that a strong field maximizes MHD-related effects.

Name (m) (m) BH singularity Energy evolution
B0-P--L0 213 2880 580 Puncture
B0-P--L1 266 2285 467 Puncture
13-P--L0 213 2880 580 Puncture 13
13-P--L1 266 2285 467 Puncture 13
13-P--L2 332 1806 376 Puncture 13
B0-M--L1 178 2360 573 Multipatch
13-M--L1 178 2360 573 Multipatch 13
B0-M-Ent-L0 138 3128 763 Multipatch Entropy
B0-M-Ent-L1 178 2360 573 Multipatch Entropy
B0-M-Ent-L1r 178 1675 573 Multipatch Entropy
13-M-Ent-L0 138 3128 763 Multipatch Entropy 13
13-M-Ent-L1 178 2360 573 Multipatch Entropy 13
13-M-Ent-L1r 178 1675 573 Multipatch Entropy 13
13-M-Ent-L2 231 1790 426 Multipatch Entropy 13
36-M-Ent-L1 178 2360 573 Multipatch Entropy 36
36-M-Ent-L2 231 1790 426 Multipatch Entropy 36
36-P--L1 266 2285 467 Puncture 36
Table 1: A list of simulations reported. Simulations vary by grids, numerical methods, and strength of seed field.
The cube of is the total number of grid points.
and are the radial and vertical grid spacing, respectively, on the equator at the radius of the initial density maximum.
See Appendix for details.

Iii Numerical Methods

Previous SpEC hydrodynamics simulations evolved fluids on Cartesian grids with points inside a radius inside the BH horizon excised, resulting in an irregular-shaped cubic-sphere or “legosphere” excision region. Points within a stencil of were evolved with one-sided differencing. This proved numerically unstable for MHD evolutions–an unsurprising result given the presence of incoming characteristic speeds on legosphere boundary faces.

Figure 1: Comparison of the total density-averaged entropy for different numerical methods for nonmagnetized and stronger magnetized cases. The method for energy evolution shows extra heating happening at the late time evolution for the nonmagnetized case using both puncture and multipath methods.
Figure 2: Convergence test on the specific entropy in the first of evolution for the magnetized disk with puncture-tau and multipatch-entropy methods. A smoothing of scalar primitive variables after the interpolation onto multipatch grids causes slightly higher initial average entropy in these runs, but the difference quickly decreases and the subsequent evolutions for both methods are in good agreement at different resolutions.
Figure 3: Comparison of the magnetic energy to the thermal energy ratio for different magnetized cases with different methods. The 36-P--L1 case reaches the same saturation level as the stronger field cases, but the magnetic energy starts to dissipate in the 36-M-Ent-L2 case, leading the ratio to decrease over time and finally saturate at a lower level (see table 1 for simulation labels).

We implemented two fixes to enable stable magnetized inflow into the BH. The first is to map to a new coordinate system in which the sphere is mapped to a point, so that the interior of this sphere is not on the grid (”excision by coordinates”). This method is implicitly used in non-vacuum numerical relativity moving puncture evolutions Baiotti and Rezzolla (2006); Faber et al. (2007); Liu et al. (2009) and has been explicitly used for MHD by Etienne et al Etienne et al. (2010a). We then evolve the MHD equations as in Muhlberger et al. (2014) with constrained transport and no explicit excision. We call this a “puncture” method. The second fix is to replace Cartesian grids with cubes deformed so as to fit together and fill the space between inner and outer spherical shells, the so called “cubed-sphere” configuration which has already been successfully applied by other codes to numerical relativity Thornburg (2004); Schnetter et al. (2006); Reisswig et al. (2007); Pazos et al. (2009); Pollney et al. (2011), hydrodynamics Zink et al. (2008); Korobkin et al. (2011, 2012); Reisswig et al. (2013), and MHD Koldoba et al. (2002); Fragile et al. (2009); Romanova et al. (2012). Each deformed sphere is evolved on its local coordinate system. We call this method “multipatch”. For the induction equation, we implement a centered hyperbolic divergence cleaning method. Details of these methods and code tests are provided in the Appendix.

An additional numerical challenge is posed by the nonmagnetized disk which, as it cools, becomes more supersonic. In our conservative MHD formulation, only the total energy and momentum density are evolved, so it becomes difficult to accurately extract temperature information when internal energy is much less than kinetic energy. SpEC has a procedure Muhlberger et al. (2014) for “fixing” energy and momentum evolution variables when they fail to map to any physical temperature and velocity. In previous papers, this fixing was invoked only in unimportant low density regions, but here it leads to glitches in temperature inside the high-density region of the torus and unphysical heating. We cure this problem by introducing an auxiliary entropy variable used to exclude unphysical jumps in temperature, similar to a technique used in the HARM3D code Noble et al. (2009). Details are given in the Appendix.

A list of the combinations of methods and resolutions reported in this paper is provided in Table 1. A comparison of results for the average entropy evolution is given in Fig. 1. Entropy is a particularly useful diagnostic of thermal evolution because it responds only to physical heating and cooling effects. Unlike temperature, entropy is unaffected by adiabatic expansion/compression and by nuclear reactions (if, as here, the gas remains in nuclear statistical equilibrium). We see that the methods give overall agreement, except that only simulations with the new entropy variable can maintain cooling of the nonmagnetized disk.

In Fig. 2, we test convergence of magnetized disk runs by evolving with both grid types at three resolutions. Fortunately, puncture and multipatch runs seem to converge to each other. Puncture grids have more gridpoints for a given resolution of the disk interior, but they also allow larger timesteps (because they don’t have the multipatch code’s concentration of angular grid points near the horizon).

For the nonmagnetized disk evolution, we have investigated the effect of numerical viscosity on the late-time cooling rate. We evolve in multipatch mode at three resolutions (the same as in the magnetized disk convergence test). We also perform a fourth simulation with a radial map that concentrates resolution near the maximum-density ring, increasing resolution there by a factor of 2.5. (See Appendix for details.) In all cases, the entropy curves, and especially the late-time cooling slopes, are nearly identical. We conclude that numerical viscosity cannot be an important part of the energy budget for this disk’s evolution.

In the magnetized disk simulations it is essential to resolve the MHD instabilities to capture all MHD effects. Resolving the MRI requires high resolution ( grid points, to capture the growth of fastest-growing mode, along  Shibata et al. (2006)). We achieve this resolution despite a modest number of grid zones by using strong seed fields and by using coordinate maps to increase the resolution in high density regions near the disk midplane as described in the Appendix. Measuring at the initial time shows that MRI fastest-growing mode is resolvable in over of the magnetized fluid (medium resolution). We find that the thermal evolution is much more sensitive to vertical than to radial resolution, presumably because it is the mode of the axisymmetric MRI with vertical wavenumber that is most significant in the high-density region, so we use grids with .

Although it is simple to check , it is not possible to disentangle MRI-driven field amplification from other effects. Local field amplification on the orbital time is seen–in fact, it is seen even in some regions where the MRI fastest-growing mode is certainly not resolved, as would be expected from nonmodal shearing wave amplification Squire and Bhattacharjee (2014). In our case, there is the additional complication that our initial state is not a hydrodynamic equilibrium but an extremely dynamical configuration.

Fig. 3 shows another comparison of results for the total magnetic energy to the thermal energy ratio for the magnetized cases. There is a good agreement between the puncture and multipatch methods for the stronger field case; The energy ratio grows by more than one order of magnitude and saturates at the same level for both methods. For the weaker field case, the multipatch run does not resolve MRI growth as well, so there is a larger difference. For the puncture run, the magnetic field saturates at the same level as the stronger field, indicating that the saturation state is independent of the initial seeded magnetic field (at least for our range of seed fields). The weakly magnetized-multipatch simulation (case 36-M-Ent-L2 in table 1) on the other hand, tracks the similar puncture simulation 36-P--L1 for about , and then it decreases for about and finally saturates at a level that is lower by a factor of two. This shows that our puncture method can resolve the magnetic field growth better for weakly magnetized case. Based on the methods comparison and convergence studies, we present puncture simulation for the weakly magnetized case (), and multipatch simulations for the nonmagnetized and strongly magnetized () cases in the next section.

Iv Results

We concentrate only on the results of simulations using multipatch grid and auxiliary entropy evolution methods with moderate resolution for the nonmagnetized case (B0-M-Ent-L1), and high resolution for strongly magnetized case (13-M-Ent-L2), and the puncture evolution methods for the weakly magnetized case (36-P--L1) in Table 1. All grids and evolution methods give similar results for the first 25 ms, but these particular runs give more reasonable results in the subsequent evolution (see the detailed discussion in appendix A). At the initial time, the thermal timescale is estimated as  ms. We evolve for about 50 ms, long enough to see the disk altered by thermal effects.

iv.1 Dynamical evolution

In Fig. 4, we plot several global quantities of the disk. As expected, adding a magnetic field enables angular momentum transport, leading to an accretion rate roughly one order of magnitude higher than that of the disk evolved without a magnetic field. The accretion rate does not appear very sensitive to the strength of the seed field, at least for the very limited range studied here. The settled accretion rate of is low enough that a thermal instability is not expected Janiuk et al. (2007). MHD effects can also cause the disk to expand radially and vertically, as is expected from angular momentum transport (see the 2D images of the density profile Figs. 6 and  5 showing the nonmagnetized and magnetized disks at  ms respectively). This transport especially drives matter into the inner radii (), leading to higher densities there. The nonmagnetized disk, on the other hand, contracts vertically and radially, becoming more ring-like as it loses thermal pressure support. Evolution without a magnetic field leads to a significantly denser disk, which explains why the magnetized disks have lower average temperature even though they have higher average entropy (last panel in Fig. 4).

The lower three panels of Fig. 4 show the effect of magnetic fields on the average entropy per baryon , electron fraction , and temperature . Even with no magnetic field, cooling (as measured by ) is delayed 10 ms by shock heating; once the disk has settled, it commences cooling. If a seed field is introduced, increases with time. The slope for the first 10 ms is higher and quite seed field-strength dependent and should perhaps be considered a transient as the field saturates, while subsequent heating is slower and less sensitive to seed field strength.

With no magnetic field, decreases monotonically, continuing the behavior seen in our earlier simulations Foucart et al. (2014), while magnetized runs show a leveling off and slight increase. Siegel and Metzger’s 3D magnetized disk simulations also find that the inner disk remains neutron rich Siegel and Metzger (2017). Radial profiles of , displayed in Fig. 7, show that the magnetized disk has higher mostly in a region around radius  km. This can be understood from the equilibrium electron fraction . In this region, the magnetized disk has lower and higher . As shown in Fig.18 of Foucart et al. (2014), increases with and decreases with , so the higher is consistent with . The outer regions of the disk, on the other hand, are too cool for to equilibrate on the simulated timescale.

Figure 8 shows gas pressure, total and poloidal magnetic pressures versus distance from the black hole at two times in the two configurations (top panel) and (bottem panlel). The toroidal field quicky grows to be the dominant component, contributing about of the total magnetic energy at the late time evolution for our stronger field case (as seen in either puncture or multipatch runs). This figure also shows that the total field pressure exceeds the gas pressure in the inner regions. Because of our rather large seed fields, the field can only grow one to two orders of magnitude before reaching overall equipartition with the internal energy. Strong toroidal fields can suppress the MRI, especially at low wavenumbers Pessah and Psaltis (2005), and this suppression may take place in some regions of our disk.

Figure 4: The evolution of the accretion rate (first panel from the top), the electron fraction (second panel), the specific entropy (third panel), and temperature (last panel) for magnetized 13-P--L1 (dashed dot line), 36-P--L1 (dashed line), and nonmagnetized B0-M-Ent-L1 cases (solid line). The accretion rate is higher by about one order of magnitude for the magnetized cases due to magnetorotational instablity. The entropy grows higher as a result of effective viscous heating, while the temperature decreases over time because of adiabatic cooling for the magnetized cases.
Figure 5: Snapshot of the rest-mass density in the meridional x-z plane at for 13 case. The solid line shows magnetic field magnitude contours correspond to [, ,,,] G.
Figure 6: Snapshot of the rest-mass density in the meridional x-z plane at for the nonmagnetized case.
Figure 7: Vertically and azimuthally averaged electron fraction profiles at the initial time, , and . decreases at densest regions and increases at low-density outer radii regions in nomnagnetized case. In the magnetized case, starts decreasing in the high-density regions only at late times.
Figure 8: Vertically and azimuthally averaged radial profiles of gas pressure, total and poloidal component magnetic pressures at , and for (top panel) and (bottom panel) magnetized cases. The toroidal component becomes dominant after the magnetic field saturates in all the magnetized cases.

iv.2 Neutrino emission and optical depth

Figure 9: Neutrino luminosity evolution for electron-flavor neutrinos and antineutrinos in the and nonmagnetized simulations. The magnetized run has systematically higher electron neutrino and antineutrino luminosities.
Figure 10: The energy-averaged optical depth of electron neutrinos radial profiles at the initial time, , and . The nonmagnetized case becomes more opaque in the high-density region. This makes the neutrino cooling less efficient than in the magnetized case, which becomes more transparent during the evolution.

Fig 9 shows the neutrino luminosity for electron-flavor species. The electron antineutrino luminosity is the strongest in both the magnetized and nonmagnetized cases. The neutrino luminosity is higher in the magnetized case for all the species, but the changes are not as large as might have been expected. In all cases, the total neutrino luminosity drops from about to a few times over about after merger. Radial emission profiles show that the luminosity drops by a comparable factor throughout the high-density region; The drop in emission does not reflect some local effect, but rather the global evolution of the disk: the contributions to the luminosity are distributed smoothly throughout the high-density region.

One possible influence on would be a change in the neutrino optical depth. Figure 10 plots the energy-averaged optical depth of electron neutrinos (the only neutrino flavor with optical depth sometimes greater than unity). The nonmagnetized disk maintains an optical depth of a few, while spreading of the magnetized disk makes it optically thin. Our disk has too low density to show the optically thin to optically thick transition from the inner radii to the outer radii seen in some alpha disk studies Lee et al. (2005); Janiuk et al. (2007). On the other hand, the total neutrino luminosity, and the fact that electron anti-neutrino emission is brightest, are consistent with the literature for disks Lee et al. (2005); Fernández and Metzger (2013).

Based on the accretion rate observations for the strongly magnetized case we estimate the equivalent parameter for our disk is . This is measured from the accretion rate, which in fact, includes the angular momentum transport due to the magnetically-driven winds and the MRI turbulence. The accretion efficiency for the stronger magnetic field case is . This efficiency is a few percent higher than the optically thin NDAF disk models () with high spin black holes as reported by Shibata et al. (2007) Shibata et al. (2007)

iv.3 Thermal evolution

The transport mechanism in an accretion disk affects the luminosity in two ways. By heating the disk, it tends to increase the luminosity. By spreading the disk to larger radii and lower densities and by facilitating higher accretion rates onto the black hole, it tends to decrease the luminosity. For a thin alpha disk, , so the former effect should initially dominate, but our disk is quite thick (), so the timescales on which these effects operate are not well separated. To understand the actual disk evolution, we must quantify the major heating and cooling effects.

Figures 1112, and 13 show the major entropy sources and sinks for different levels of initial magnetization. From the energy and lepton number source terms provided by the leakage code, a radiative entropy sink term can be computed (see Appendix for details). Cooling from advection into the black hole is straightforwardly measured by monitoring entropy flux at the inner boundary. Adiabatic expansion and nuclear reactions (in nuclear statistical equilibrium) do not affect entropy, while shocks, reconnection, and turbulent dissipation should only heat. Thus, the total heating rate should be


where is the time derivative of the total entropy. Unfortunately, it is difficult to separate the various possible heating sources, as they will all appear in the code via the stabilizing dissipation terms in our shock-capturing MHD scheme. We normalize each source term by the instantaneous total entropy of the disk, giving the source terms the quality of inverse timescales.

The entropy budget plots Figs. 1112, and 13 tell a clear story. At early times, there is strong heating in all cases from shocks as the disk, still nonlinearly perturbed from equilibrium, pulsates and axisymmetrizes. This heating ceases about 30 ms after merger as the disk settles. It is especially clear in the nonmagnetized case (Fig. 11) that this happens before the neutrino luminosity drops. The neutrino luminosity drops quickly thereafter, on a fraction of the initial thermal timescale, as radiation cools the disk enough to decrease itself. This rapid cooling stops when the thermal timescale has increased to about  ms.

It is worth mentioning that the unmagnetized simulations show this initial strong heating regardless of the grid and methods used. Indeed, it is seen even in the original simulation presented in Foucart et al. (2014). The exact amount of early-time heating does vary noticeably from one method to another. It is greatest for the multipatch runs, perhaps because of numerical perturbations caused by switching to a radically different grid (Fig. 1).

For the nonmagnetized case, the final state is neutrino cooling-dominated. Accretion has nearly stopped, and advective cooling is negligible. The heating rate is significantly lower than the neutrino cooling rate, although the average of the former is still around a third of the latter. Note that the heating rate does occasionally become negative, presumably a sign of numerical error in the difficult-to-follow thermal evolution of the gas as it becomes ever more supersonic. This negative heating could be removed by a stricter lower limit on the entropy (see the Appendix for details), but this would bias numerical error toward heating, which might have an undesirable cummulative effect.

For magnetized cases, the heating rate remains above the neutrino cooling rate. However, this effect is largely cancelled by the strong advective cooling that takes place as hot material accretes into the black hole. Although the component entropy sources and sinks are larger than in the nonmagnetized case, the thermal timescale in these cases also increases to . At late times, the neutrino luminosity decreases slightly faster in the most highly magnetized case, although in all cases a luminosity of around erg s will be maintained till the end of the evolution.

Figure 11: The evolution of the total heating rate , neutrino cooling rate and advection cooling rate ratios to the total entropy for the nonmagnetized case. Heating and neutrino cooling rates drop significantly around . The advection cooling is almost zero at the end of the simulation.
Figure 12: The evolution of the total heating rate , neutrino cooling rate and advection cooling rate ratios to the total entropy for the case with weaker seed field (). Like in the nonmagnetized case, the total heating and neutrino cooling rates decrease significantly comparing with the early time. The advection cooling rate is considerably higher due to the MHD effects.
Figure 13: The evolution of the total heating rate , neutrino cooling rate and advection cooling rate ratios to the total entropy for the standard, stronger seed field. The disk shows the same qualitative thermal behavior as in the weaker seed field case.

In Figure 14, we plot late-time convergence for representative global quantities. Convergence at these times is difficult to achieve because it requires resolving not only the fastest-growing MRI wavelength but also sufficient inertial range that the average transport effects of turbulence are accurately captured. In the highest-resolution study to date of a BHNS post-merger system, Kiuchi et al. Kiuchi et al. (2015) were unable to demonstrate convergence even with grid spacing a few times smaller than we can afford. Thus, it is not surprising that we also obtain no better than qualitative convergence, i.e. the overall behavior is similar at all resolutions. Like Kiuchi et al., we see a tendency toward more vigorous turbulent heating at higher resolutions.

Figure 14: Convergence study of the total neutrino luminosity and heating rate for nonmagnetized and strong seed field multipatch runs. The results are qualitatively convergent in both cases.

iv.4 Comparison with previous studies

Magnetized black hole-neutron star mergers have been carried out by other groups Chawla et al. (2010); Etienne et al. (2012); Paschalidis et al. (2015); Kiuchi et al. (2015). In particular, Etienne et al. Etienne et al. (2013) have also inserted a poloidal field into a post-merger BHNS disk. The highest-resolution MHD BHNS merger simulation is that of Kiuchi et al Kiuchi et al. (2015). These simulations, like ours, have more realistic initial disk profiles than analytic tori would provide. However, there are important differences in our treatment of the thermal evolution of the disk. The above-mentioned studies, since they did not employ a finite-temperature nuclear equation of state, did not include neutrino cooling, which is present in our simulations, so their disks were presumably too hot. On the other hand, by inserting a seed field only when it was safe to apply the Cowling approximation, our disk had cooling for 15 ms without one of the major heating sources, so our disk is likely over-cooled. Furthermore, the convergence studies in Kiuchi et al. (2015) find that heating in the inner disk is very resolution dependent, with insufficient resolution leading to underestimates of thermally driven outflows in particular. A truly realistic thermal state of the disk would presumably be somewhere in between these extremes. (Of course, a truly realistic treatment would also require neutrino transport, not just a leakage approximation.)

Both the Etienne et al poloidal seed study and the highest-resolution simulations of Kiuchi et al Kiuchi et al. (2015) find sustained Blandford-Znajek Poynting flux polar jets. Our MHD disk evolutions do produce some unbound outflow over the simulation period (, average  s) and a magnetically-dominated polar region, but the polar field does not organize itself into a radial Blandford-Znajek like structure. The above-mentioned differences in thermal treatment in our simulation and resolution effects may play a role here. However, we also find that the presence or absence of strong winds and Poynting flux outflows is sensitive to the choice of seed field. When we evolve with a less confined initial field, we do see stronger matter outflow and polar magnetic flux. One might worry that these effects could somehow be suppressed by too strong a seed field, but simulations of neutrino-cooled disks with analytic initial conditions all find that stronger fields (as high as ) yield stronger winds and stronger Blandford-Znajek luminosity Janiuk et al. (2013). Even these analytic disks with weaker initial magnetic field () find strong unbounded outflows after evolving the magnetized disk long enough Siegel and Metzger (2017), and this might also turn out to be the case for our more confined B-field evolutions if evolved longer.

V Conclusions

We have carried out simulations of a BHNS post-merger system with a realistic initial state provided by a numerical relativity merger simulation , including both neutrino emission effects and magnetic field evolution. The initial magnetic field is applied as large poloidal loops confined in the post-merger disk. Because our simulations include the major heating and cooling sources, we can study the contribution of each thermal driving process as the disk settles toward thermal equilibirum. Without a magnetic field, there is no such thermal equilibrium, so after an initial phase of shock heating, the disk enters a phase of long-term cooling by neutrinos. With a strong seed magnetic field, the final state after several initial thermal timescales is a rough balance between MHD-related heating and advective cooling, with neutrino cooling being a secondary effect, driving the entropy down over longer timescales. This is roughly consistent with the long term evolution of two dimensional neutrino cooled -viscosity disks reported by Fernandez et al. Fernández et al. (2015), where neutrino cooling is only important at early times. In both magnetized and nonmagnetized cases, the main reason for settling is not a precise achievement of equilibrium, but an increase in the thermal timescale (from  ms to  ms) as the initially-high neutrino luminosity drops.

The considered magnetized 3D BHNS post-merger configuration provided the opportunity to test multiple methods for evolving the relativistic MHD equations. These show reassuring consistency over the first  ms, but realistic long-term evolution requires careful treatment of the energy variable, especially in how one handles the problematic recovery of primitive variables. The multipatch methods employed in some of our simulations can easily be applied to more general grid configurations Pazos et al. (2009).

The initial study of magnetized 3D BHNS post-merger disk evolution presented in this paper is limited in many ways. Only one BHNS system and one magnetic seed field geometry were used. Neutrino effects might be different for an opaque disk (e.g., Deaton et al. (2013)), and magnetohydrodynamic effects are known to be seed field-dependent Beckwith et al. (2008). Our leakage scheme neglects neutrino absorption, which could smooth temperature profiles and launch winds. Existing neutrino transport codes (e.g., Foucart et al. (2015)) can in the future be used to capture these effects. Finally, it would be interesting to carry out a similar study on NSNS post-merger systems.

The authors thank Zachariah Etienne, Scott Noble, Vasileios Paschalidis, Jean-Pierre De Villiers, John Hawley, and Hotaka Shiokawa, for helpful discussions and advice over the course of this project. M.D. acknowledges support through NSF Grant PHY-1402916. F.H. acknowledges support from the Navajbai Ratan Tata Trust at IUCAA, India. F.F. acknowledges support from Einstein Postdoctoral Fellowship grant PF4-150122, awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. H.P. gratefully acknowledges support from the NSERC Canada. L.K. acknowledges support from NSF grants PHY-1306125 and AST-1333129 at Cornell, while the authors at Caltech acknowledge support from NSF Grants PHY-1404569, AST-1333520, NSF-1440083, and NSF CAREER Award PHY-1151197. Authors at both Cornell and Caltech also thank the Sherman Fairchild Foundation for their support. Computations were performed on the Caltech compute clusters Zwicky and Wheeler, funded by NSF MRI award No. PHY-0960291 and the Sherman Fairchild Foundation. Computations were also performed on the SDSC cluster Comet under NSF XSEDE allocation TG-PHY990007N.

Appendix A Numerical improvements

a.1 Formulation

The fundamental equations to be evolved are the same as in our earlier MHD work Muhlberger et al. (2014). We write the metric


The fluid at each grid point is described by its set of “primitive variables”: baryonic density , temperature , electron fraction , and spatial components of the covariant 4-velocity . From , , and , the equation of state supplies the gas pressure , specific enthalpy , and sound speed . From , we know the Lorentz factor . The stress tensor is


where is the Faraday tensor. We assume a perfectly conducting fluid, , which fixes the electric field. The variables actually evolved (aside from the magnetic field, whose evolution is described below) are the conservative variables: a density variable , the proton density , an energy density variable , and a momentum density variable . In the above, is the determinant of the spatial metric. We evolve using an HLLE approximate Riemann solver Harten et al. (1983). Conservative formulations have the advantage that numerical dissipation in shock or turbulent subscale structures is automatically conservative. They have the disadvantage of not evolving a separate variable for the internal energy or entropy. Such information must be recovered by root finding from the conservative variables after each timestep, which can be expensive and (especially if kinetic energy dominates over internal energy in ) inaccurate.

The magnetic field can be described via the components of its 2-form or its vector field , related as . In a conducting medium, field lines advect with the fluid: . Since , we can alternatively evolve the vector potential 1-form , where . A vector potential evolution will automatically satisfy but will require specifying a gauge.

Our Cartesian grid simulations suppress monopoles via a constrained transport scheme, which requires staggering or between gridpoints. For the multipatch simulations described below, this would be very inconvenient because the patch coordinate transformations would have to account for each component of the field being at a different location, so we instead code two well-known methods that control while keeping all variables centered at the same gridpoints. The first is a centered vector potential method, implemented as in Giacomazzo et al. (2011). We find that the generalized Lorentz gauge, introduced in Farris et al. (2012), provides the best stability. The evolution for and the scalar potential are given by


where is a specifiable constant of order the mass of the system. Lorentz-type gauges lead to luminal characteristic speeds, but fortunately the speeds used in the HLLE fluxes used in the evolution of (see Giacomazzo et al. (2011)) can still be set to the physical, MHD wave maximum speed. The signal speeds for HLLE fluxes in the evolution, on the other hand, are set to the null .

The second magnetic evolution scheme is a covariant hyperbolic divergence cleaning method Liebling et al. (2010); Penner (2011); Mösta et al. (2014), in which an auxiliary evolution variable is introduced to damp monopoles. The Maxwell equation is replaced by , where is the 4-metric, the Faraday tensor, the unit time vector, and a specifiable damping constant. In components


where we set . Eq. 7 is in conservative form and can be evolved using our usual HLLE scheme, while Eq. 8 is evolved via straightforward second-order centered finite differencing.

Both of these methods require added numerical dissipation. Thus, we add Kreiss-Oliger dissipation to the magnetic evolution equations. For multipatch simulations, this step is done while time derivatives are being computed in the local patch coordinate system of evolution variable components in these coordinates.


is for divergence cleaning and for the vector potential method. is a second-derivative operator, and is the grid spacing in the -th direction, both computed in local patch coordinates. is a function of space, which vanishes on boundary points but may be otherwise chosen according to the problem Mattsson et al. (2004).

a.2 Cubed-sphere Multipatch Grids

Several groups have already implemented dynamics on spherical surfaces Sadourny (1972); Ronchi et al. (1996), 3D hydrodynamics Zink et al. (2008); Korobkin et al. (2011, 2012); Reisswig et al. (2013), 3D MHD Koldoba et al. (2002); Fragile et al. (2009); Romanova et al. (2012), and Einstein’s equations Thornburg (2004); Schnetter et al. (2006); Reisswig et al. (2007); Pazos et al. (2009); Pollney et al. (2011) with multipatch methods and cubed-sphere-like grids. The basic idea is to divide the computational domain into patches, each of which has its own local coordinate system in which it is a uniform Cartesian mesh. In the global coordinate system, each patch is distorted, and six distorted cubes can be fit together to fill a volume with spherical inner and outer boundaries. Time derivative calculations for timesteps are computed within the local patch coordinates and then transformed to the global coordinate system. Multipatch methods easily generalize to any combination of distorted cubes. For example, the central hole can be filled with a cube (as done in a test problem below), or the cubed-sphere could be surrounded by non-distorted cubes.

This method can be contrasted with other popular ways of evolving grids around black holes. One is the use of spherical-polar coordinate grids. The second is the use of Cartesian grids, with removal of the black hole interior accomplished either by excising all gridpoints within a spherical region (leading to an irregular-shaped “legosphere”) or by removing the interior via a radial coordinate transformation (“puncture”) Etienne et al. (2010b). All previous SpEC black hole-neutron star simulations use Cartesian grids with legosphere excision. We have been unable to find a stable implementation of this method for magnetized flows into a black hole. This is not surprising, since Cartesian grid faces even inside the horizon will have characteristic fields flowing into the grid, making the evolution ill-posed without boundary conditions providing information about the excised interior. Both spherical-polar and multipatch grids can naturally excise spherical regions (which can be distorted by coordinate transformations to fit the horizon shape as needed) and have no incoming characteristics if placed inside the apparent horizon (and outside the Cauchy horizon) of a stationary black hole. Multipatch methods have an advantage over spherical-polar grids that they do not suffer from coordinate singularities and grid pileup near the poles, which can be an issue for high-resolution spherical-polar simulations Shiokawa et al. (2012). Spherical-polar grids, on the other hand, have two advantages. First, for nearly axisymmetric systems, one can have much lower resolution in longitude than in latitude, a freedom not present in multipatch grids. Second, communication between patches in multipatch grids is by ghost zone overlaps. Ghost zone gridpoints will not match gridpoints on the overlapping live patch, so they must be filled by interpolation. This introduces a new source of error which will generally not exactly respect conservation laws and may create magnetic monopoles, although it should converge away with resolution. Which method is best most likely depends on the problem.

Since our conservative evolution equations are generally covariant, it is straightforward to evolve them in the local patch coordinates, shifting to global coordinates for ghost zone synchronization. For our WENO5 reconstruction method, we need three ghost zone layers on patch interior boundaries. Because of our methods of “fixing” problematic points described below, synchronizing variables is not quite the same as just synchronizing their time derivatives, and we find the former to be needed for stability. For the divergence cleaning method, any monopoles generated by interpolation of in ghost zones are damped (by design of Eqs 7 and 8) and remain small. For the vector potential method, we synchronize computed from on patch faces, where information is lacking on one side to compute the curl. It is crucial here to synchronize only the outermost layer of points, not the full 3-layer ghost zone region, because the latter will introduce monopoles in the ghost zones and lead rapidly to an instability there.

Our cubed-sphere grids are largely the same as those of other groups. A minor alteration in the ghost zones is illustrated in Fig. 15 to eliminate the presence of overlap regions which are “live” for both grids (i.e. neither is synced with respect to the other). Our fears that “live overlaps” would be dangerous have not been borne out, but the new arrangement does seem to propagate shocks a bit better and show less deviation in rest mass (interpolated ghost zones do not allow strict mass conservation in either case), although it cannot be generalized to more general multipatch structures.

Figure 15: An illustration of the synchronization of ghost zone regions at internal patch boundaries. As is standard practice with uniform grids, one grid is extended the full ghost zone width (three points, in our case) beyond the interface, while the other grid extends two points. For a cubed spheres setup, ghost zone extensions must be chosen to guarantee sufficient overlaps on 3-patch edges. Above, open circles are ghost zone points; filled circles are live points. The points B and H overlap and mark the interface. First H is set to B (which does not require interpolation). Then H can be used in the interpolation to get C.

Figures 16, 17, and 18 show some standard MHD test problems applied to the multipatch MHD code. Fig. 16 is the first Riemann problem from Balsara (2001) and Giacomazzo and Rezzolla (2007), containing a left-going fast rarefaction wave, a left-going compound wave, a contact discontinuity, a right-going slow shock and a right-going fast rarefaction wave. To test relativistic terms, we set lapse and shift , yielding the expected slowdown and advection. A cubical patch is added to the center to fill the inner hole, while the planar symmetry is imposed on the outer boundary, setting functions in the outer points to their values at the closest point in the interior on a line in the symmetry plane. The waves travel through interpolated boundaries without incident.

In Fig. 17, we evolve a Bondi accretion problem (the same as in Muhlberger et al. (2014)) with a radial magnetic field and maximum of 2.5. Both of these tests are performed with the divergence cleaning code. As in our earlier paper Muhlberger et al. (2014), we find better behavior when we add Kreiss-Oliger dissipation to all variables, with . Errors saturate after a few of evolution, with second-order convergence demonstrated except at the sonic point and the inner boundary.

Finally, we evolve a constant angular momentum Fishbone-Moncrief torus  Fishbone and Moncrief (1976) around a rapidly spinning black hole. We set the dimensionless spin of the black hole to , the angular momentum parameter to (), average , and the equation of state to a Gamma law with , making the problem very similar to a standard scenario studied by the HARM code McKinney and Gammie (2004); Shiokawa et al. (2012). We use the same radial and angular coordinate maps as in these studies. Like in Fragile et al. (2009), we find that it is necessary to tilt the grid in order for the current sheet formed by winding of the seed field to break in a reasonable time and initiate turbulence. We agree with their observation that this is an artifact of symmetries in the setup and should not be a worry for general problems. Fig. 18 shows on the left a snapshot of the density at , on the right a representation of the grid with resolution quartered for clarity. The actual evolution grid used 120 radial points and 60 angular points across each of the six patches. For this problem, we found it advantageous to have higher dissipation in problematic regions (low density regions and the viscinity of the black hole) and low dissipation inside the torus, where we wanted to resolve the MRI with modest resolution. There we set () inside the disk for (), and we set for or . Results qualitatively match the literature, with mass flow into the horizon , electromagnetic energy flux out of the horizon , and the generation of unbound matter.

Figure 16: Magnetized Riemann problem evolved on both a cubed-sphere multipatch grid (with about 240 grid points across the diameter of the spherical computational domain), together with the results for the same problem evolved on a Cartesian one-dimensional grid which is able to utilize the planar symmetry. The interface between the inner cube and outer cubed spheres is at 0.27.
Figure 17: Convergence test for Bondi accretion with a radial magnetic field. Shown here is the error in the conservative variable at , by which time it has settled. The error plotted is the absolute change in . The relative change of is about at the lower resolution. Second-order convergence breaks down at the sonic radius at , as expected. The grid consists of 48 domains, with each of the six patches split in two on each of its axes.
Figure 18: Meridional snapshot at of a turbulent accretion torus. Shown on the left is the density (on a logarithmic scale covering the four decades up to the maximum). On the right is the grid, with resolution reduced by about a factor of 4 for clarity. The effects of the radial and angular maps are visible, as is the regularity of the poles. The unusually close radial lines are nonmatching ghost zones.

a.3 Coordinate maps

In its fluid module, SpEC assumes uniform grid spacing in the coordinates on which the grid is defined, so nonuniformity can be acheived by introducing coordinate transformations between these grid coordinates and the original, “physical” coordinates. All simulations use logarithmic radial maps [i.e. uniform spacing in ], concentrating grid near the black hole. For Cartesian simulations, this naturally introduces a puncture, set at the desired excision radius . It leads to enormous distortions on the edges of the cubical grid, but since we only evolve in a sphere contained by the cube, this causes no problems. For multipatch simulations, the exponential map preserves the ratio between radial and transverse grid spacings; both increase with distance from the center.

We also add maps to concentrate grid near the equator. For multipatch runs, we use the angular map common for MHD disk simulations  McKinney and Gammie (2004) with . For Cartesian runs, this angular map unacceptably distorts grid cells, leading to artifacts in the evolution, so we instead use a cubic scale map on the z axis ( with and ).

Finally, we have carried out multipatch simulations using a radial map (composed with the logarithmic map) to concentrate grid on a ring coinciding with the high-density region. The map has the form


where is the grid radius, the physical radius, controls the width of the zoomed region, while , , and are set so that and coincide at the inner and outer radii, and the appropriate zoom factor () is acheived at .

a.4 Primitive variable recovery

Sometimes, due to numerical error, the evolved conservative variables may not correspond to any physical . In this case, we can “fix” the conservative variables to make primitive variable recovery possible using the prescription described in Appendix A of Muhlberger et al. (2014) (straightforwardly altered to take into account the minimum of being less than one Deaton et al. (2013)). Unfortunately, this introduces glitches in supersonic flows such as those in thin disks, usually seen as gridpoints at which the temperature discontinuously jumps to the equation of state table minimum. Although this is initially a cooling effect, the glitches create artificial heating. For nonmagnetized disk simulations, this ultimately stalls the cooling of the disk after only a small decrease in total entropy.

We remove this problem by introducing an auxiliary entropy evolution variable , where is the entropy per baryon. The use of entropy variables to reset problematic gridpoints and ameliorate accuracy problems in the evolution of internal energy by conservative codes has already been tried by other groups Balsara and Spicer (1999); Ruffert et al. (2001); Koldoba et al. (2016).

In the absence of subgrid-scale energy dissipation (shocks, reconnection, turbulence), the entropy of a fluid in nuclear statistical equilibrium evolves by advection and neutrino emission only Shapiro and Teukolsky (1983); Ruffert et al. (2001).


where and are the net neutrino energy and lepton number emission rates per volume, respectively, is the nucleon mass, and are chemical potentials. Note that, since we have excluded only heating effects, the evolved gives a lower bound on the true entropy.

Roughly speaking, we now have two energy variables, and , which are made to be consistent with each other at the beginning of each timesetep. Each step, we execute the following procedure.

  1. Evolve using an HLL approximate Riemann solver. must be evolved with a monotonic reconstructor to avoid new extrema. We use a second-order monotonized centered (MC2) limiter van Leer (1977). The other variables can be evolved with higher-order reconstruction like WENO.

  2. Compute , , and from the appropriate divisions of the conservative variables.

  3. If not, attempt to solve using and the other conservative variables except using the gnewton method as implemented by the GSL Scientific Library et al ().

  4. If a root is found, use it to compute the entropy, . If , accept the root. The parameter but is otherwise freely specifiable. We use .

  5. If a root was not found, or if it violates the condition in step 4, first check to see if the point is in the force-free regime. If so, use force-free recovery of . (See Muhlberger et al. (2014) for details on this solver and the conditions for its use.)

  6. If the point does not meet the force-free conditions, attempt to solve for using and the other conservative variables except , again using GSL’s gnewton. If a root was found in step 4, use it as the initial guess for the root solve. (Using instead of speeds up convergence in some difficult points, but probably makes little difference in general.)

  7. If a root could not be found with multidimensional root finding, attempt again with and other variables except , this time with GSL’s 1D brent root finder. Here we regard as the variable, solving Eq. (A24) of Muhlberger et al. (2014), with solved via a separate 1D solve of the condition on each iteration. This 1+1D solving is much slower but more robust than the 2D solver.

  8. If this fails, attempt a 1D bracketing algorithm for which uses rather than . (See Appendix A of Muhlberger et al. (2014)). If this fails, terminate the evolution with an error.

  9. If an acceptable root was found, apply other “atmosphere” fixes to the primitive variables at low densities: limits to the temperature and Lorentz factor in these regions.

  10. Recompute all conservative variables from these final primitive variables. and are now again consistent.

A simple sanity check on our implementation of the source terms in Eq. 11 is to alter the above to force the code to always use the evolved in primitive variable recovery, in which case one observes the disk cooling on a timescale of the total thermal energy divided by neutrino luminosity.

In Fig. 1, we have already shown the difference this method makes to the entropy evolution of the nonmagnetized disk. Significantly, all discontinuous artifacts are gone when the new method is used. Because the magnetized disk does not reach such low entropies, the choice of methods makes little difference for those simulations.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description