Evolution of the Magnetized, NeutrinoCooled Accretion Disk in the Aftermath of a Black Hole Neutron Star Binary Merger
Abstract
Black holetorus systems from compact binary mergers are possible engines for gammaray bursts (GRBs). During the early evolution of the postmerger remnant, the state of the torus is determined by a combination of neutrino cooling and magneticallydriven heating processes, so realistic models must include both effects. In this paper, we study the postmerger evolution of a magnetized black holeneutron star binary system using the Spectral Einstein Code (SpEC) from an initial postmerger state provided by previous numerical relativity simulations. We use a finitetemperature 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 generalrelativistic magnetohydrodynamics (MHD), including the use of cubedsphere 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 MHDrelated 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.
pacs:
I Introduction
The cause of shorthard 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 holeneutron star (BHNS) and neutron starneutron 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 baryonpoor ultrarelativistic outflow of the sort needed to explain GRB properties. The accretion gas cools primarily by neutrino emission, and so such systems are classified as neutrinodominated 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 neutrinoantineutrino 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 postmerger 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 finitetemperature 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 energyintegrated twomoment M1 transport scheme Shibata et al. (2011); Foucart et al. (2015). Comparing to the earlier times of evolution we found that the postmerger accretion disks become cold, and more neutronrich 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 lowdensity regions. No significant neutrinodriven 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.
Longterm 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). Largescale fields threading the BH ergosphere enables extraction of the black hole spin energy into a Poynting fluxdominated 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 Poyntingdominated 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 longterm 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 longterm postmerger evolution of BHNS (and highmass 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 onedimensional (axisymmetric, vertically summed), equilibrium NDAF studies Popham et al. (1999); Di Matteo et al. (2002); Janiuk et al. (2004). Onedimensional NDAFs were evolved by Janiuk et al. Janiuk et al. (2007), who found disks can become viscothermally 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 neutrinoantineutrino energy release and latetime 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 firstprinciples, 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 postmerger 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, postmerger 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 postmerger 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 finitetemperature 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 cubedsphere 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 latetime 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 MHDrelated 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 neutrinorelated 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 LattimerSwesty 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 freeemission 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 FermiDirac 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 nucleonnucleon 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 M127S9 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
(1) 
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 highdensity 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 volumeweighted integrated gas pressure to the volumeweighted 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 postmerger magnetic fields. We focus on this extreme case first for two reasons. First, it allows us to resolve the rapidlygrowing 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 fastestgrowing 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 strongfield, highdensity 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 strongseed disk simulations Janiuk et al. (2013)) that a strong field maximizes MHDrelated effects.
Name  (m)  (m)  BH singularity  Energy evolution  

B0PL0  213  2880  580  Puncture  
B0PL1  266  2285  467  Puncture  
13PL0  213  2880  580  Puncture  13  
13PL1  266  2285  467  Puncture  13  
13PL2  332  1806  376  Puncture  13  
B0ML1  178  2360  573  Multipatch  
13ML1  178  2360  573  Multipatch  13  
B0MEntL0  138  3128  763  Multipatch  Entropy  
B0MEntL1  178  2360  573  Multipatch  Entropy  
B0MEntL1r  178  1675  573  Multipatch  Entropy  
13MEntL0  138  3128  763  Multipatch  Entropy  13 
13MEntL1  178  2360  573  Multipatch  Entropy  13 
13MEntL1r  178  1675  573  Multipatch  Entropy  13 
13MEntL2  231  1790  426  Multipatch  Entropy  13 
36MEntL1  178  2360  573  Multipatch  Entropy  36 
36MEntL2  231  1790  426  Multipatch  Entropy  36 
36PL1  266  2285  467  Puncture  36 
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 irregularshaped cubicsphere or “legosphere” excision region. Points within a stencil of were evolved with onesided differencing. This proved numerically unstable for MHD evolutions–an unsurprising result given the presence of incoming characteristic speeds on legosphere boundary faces.
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 nonvacuum 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 “cubedsphere” 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 highdensity 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 latetime 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 maximumdensity ring, increasing resolution there by a factor of 2.5. (See Appendix for details.) In all cases, the entropy curves, and especially the latetime 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 fastestgrowing 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 fastestgrowing 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 highdensity region, so we use grids with .
Although it is simple to check , it is not possible to disentangle MRIdriven 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 fastestgrowing 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 magnetizedmultipatch simulation (case 36MEntL2 in table 1) on the other hand, tracks the similar puncture simulation 36PL1 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 (B0MEntL1), and high resolution for strongly magnetized case (13MEntL2), and the puncture evolution methods for the weakly magnetized case (36PL1) 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 ringlike 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 fieldstrength 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.
iv.2 Neutrino emission and optical depth
Fig 9 shows the neutrino luminosity for electronflavor 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 highdensity 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 highdensity region.
One possible influence on would be a change in the neutrino optical depth. Figure 10 plots the energyaveraged 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 antineutrino 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 magneticallydriven 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 11, 12, 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
(2) 
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 shockcapturing 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. 11, 12, 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 earlytime 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 coolingdominated. 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 difficulttofollow 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.
In Figure 14, we plot latetime convergence for representative global quantities. Convergence at these times is difficult to achieve because it requires resolving not only the fastestgrowing MRI wavelength but also sufficient inertial range that the average transport effects of turbulence are accurately captured. In the highestresolution study to date of a BHNS postmerger 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.
iv.4 Comparison with previous studies
Magnetized black holeneutron 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 postmerger BHNS disk. The highestresolution 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 abovementioned studies, since they did not employ a finitetemperature 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 overcooled. 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 highestresolution simulations of Kiuchi et al Kiuchi et al. (2015) find sustained BlandfordZnajek Poynting flux polar jets. Our MHD disk evolutions do produce some unbound outflow over the simulation period (, average s) and a magneticallydominated polar region, but the polar field does not organize itself into a radial BlandfordZnajek like structure. The abovementioned 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 neutrinocooled disks with analytic initial conditions all find that stronger fields (as high as ) yield stronger winds and stronger BlandfordZnajek 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 Bfield evolutions if evolved longer.
V Conclusions
We have carried out simulations of a BHNS postmerger 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 postmerger 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 longterm cooling by neutrinos. With a strong seed magnetic field, the final state after several initial thermal timescales is a rough balance between MHDrelated 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 initiallyhigh neutrino luminosity drops.
The considered magnetized 3D BHNS postmerger configuration provided the opportunity to test multiple methods for evolving the relativistic MHD equations. These show reassuring consistency over the first ms, but realistic longterm 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 postmerger 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 fielddependent 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 postmerger systems.
Acknowledgements.
The authors thank Zachariah Etienne, Scott Noble, Vasileios Paschalidis, JeanPierre 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 PHY1402916. F.H. acknowledges support from the Navajbai Ratan Tata Trust at IUCAA, India. F.F. acknowledges support from Einstein Postdoctoral Fellowship grant PF4150122, awarded by the Chandra Xray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS803060. H.P. gratefully acknowledges support from the NSERC Canada. L.K. acknowledges support from NSF grants PHY1306125 and AST1333129 at Cornell, while the authors at Caltech acknowledge support from NSF Grants PHY1404569, AST1333520, NSF1440083, and NSF CAREER Award PHY1151197. 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. PHY0960291 and the Sherman Fairchild Foundation. Computations were also performed on the SDSC cluster Comet under NSF XSEDE allocation TGPHY990007N.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
(3) 
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 4velocity . 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
(4) 
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 2form or its vector field , related as . In a conducting medium, field lines advect with the fluid: . Since , we can alternatively evolve the vector potential 1form , 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 wellknown 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
(5)  
(6) 
where is a specifiable constant of order the mass of the system. Lorentztype 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 4metric, the Faraday tensor, the unit time vector, and a specifiable damping constant. In components
(7)  
(8) 
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 secondorder centered finite differencing.
Both of these methods require added numerical dissipation. Thus, we add KreissOliger 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.
(9) 
is for divergence cleaning and for the vector potential method. is a secondderivative 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 Cubedsphere 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 cubedspherelike 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 cubedsphere could be surrounded by nondistorted cubes.
This method can be contrasted with other popular ways of evolving grids around black holes. One is the use of sphericalpolar 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 irregularshaped “legosphere”) or by removing the interior via a radial coordinate transformation (“puncture”) Etienne et al. (2010b). All previous SpEC black holeneutron 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 illposed without boundary conditions providing information about the excised interior. Both sphericalpolar 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 sphericalpolar grids that they do not suffer from coordinate singularities and grid pileup near the poles, which can be an issue for highresolution sphericalpolar simulations Shiokawa et al. (2012). Sphericalpolar 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 3layer ghost zone region, because the latter will introduce monopoles in the ghost zones and lead rapidly to an instability there.
Our cubedsphere 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.
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 leftgoing fast rarefaction wave, a leftgoing compound wave, a contact discontinuity, a rightgoing slow shock and a rightgoing 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 KreissOliger dissipation to all variables, with . Errors saturate after a few of evolution, with secondorder convergence demonstrated except at the sonic point and the inner boundary.
Finally, we evolve a constant angular momentum FishboneMoncrief 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.
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 highdensity region. The map has the form
(10) 
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 subgridscale 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).
(11) 
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.

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

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

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 ().

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

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 forcefree regime. If so, use forcefree recovery of . (See Muhlberger et al. (2014) for details on this solver and the conditions for its use.)

If the point does not meet the forcefree 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.)

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.

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.

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.

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.
References
 Nakar (2007) E. Nakar, Phys. Rep. 442, 166 (2007), arXiv:astroph/0701748 .
 Berger (2014) E. Berger, Ann. Rev. Astron. Astrophys. 52, 43 (2014), arXiv:1311.2603 [astroph.HE] .
 Popham et al. (1999) R. Popham, S. E. Woosley, and C. Fryer, Astrophys. J. 518, 356 (1999), arXiv:astroph/9807028 [astroph] .
 Di Matteo et al. (2002) T. Di Matteo, R. Perna, and R. Narayan, Astrophys. J. 579, 706 (2002), arXiv:astroph/0207319 .
 Janiuk et al. (2004) A. Janiuk, R. Perna, T. Di Matteo, and B. Czerny, Mon. Not. Roy. Astron. Soc. 355, 950 (2004), arXiv:astroph/0406362 .
 Paczynski (1991) B. Paczynski, Acta Astronomica 41, 257 (1991).
 Jaroszynski (1996) M. Jaroszynski, Astron. Astrophys. 305, 839 (1996), arXiv:astroph/9506062 [astroph] .
 Richers et al. (2015) S. Richers, D. Kasen, E. O’Connor, R. Fernandez, and C. D. Ott, Astrophys. J. 813, 38 (2015), arXiv:1507.03606 [astroph.HE] .
 Just et al. (2016) O. Just, M. Obergaulinger, H. T. Janka, A. Bauswein, and N. Schwarz, Astrophys. J. 816, L30 (2016), arXiv:1510.04288 [astroph.HE] .
 Blandford and Znajek (1977) R. D. Blandford and R. L. Znajek, Mon. Not. Roy. Astr. Soc. 179, 433 (1977).
 Meszaros and Rees (1997) P. Meszaros and M. Rees, Astrophys.J. 482, L29 (1997), arXiv:astroph/9609065 [astroph] .
 Shakura and Sunyaev (1973) N. I. Shakura and R. A. Sunyaev, in X and GammaRay Astronomy, IAU Symposium, Vol. 55, edited by H. Bradt and R. Giacconi (1973) p. 155.
 Lee et al. (2005) W. H. Lee, E. RamirezRuiz, and D. Page, Astrophys. J. 632, 421 (2005), arXiv:astroph/0506121 .
 Dessart et al. (2009) L. Dessart, C. D. Ott, A. Burrows, S. Rosswog, and E. Livne, Astrophys. J. 690, 1681 (2009), arXiv:0806.4380 .
 Perego et al. (2014) A. Perego, S. Rosswog, R. M. Cabezón, O. Korobkin, R. Käppeli, A. Arcones, and M. Liebendörfer, Mon. Not. Roy. Astr. Soc. 443, 3134 (2014), arXiv:1405.6730 [astroph.HE] .
 Masada and Sano (2008) Y. Masada and T. Sano, Astrophys. J. 689, 1234 (2008), arXiv:0808.2338 .
 Guilet et al. (2015) J. Guilet, E. Mueller, and H.T. Janka, Mon. Not. Roy. Astron. Soc. 447, 3992 (2015), arXiv:1410.1874 [astroph.HE] .
 Foucart et al. (2015) F. Foucart, E. O’Connor, L. Roberts, M. D. Duez, R. Haas, L. E. Kidder, C. D. Ott, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi, Phys. Rev. D91, 124021 (2015), arXiv:1502.04146 [astroph.HE] .
 Kiuchi et al. (2015) K. Kiuchi, Y. Sekiguchi, K. Kyutoku, M. Shibata, K. Taniguchi, and T. Wada, Phys. Rev. D 92, 064034 (2015), arXiv:1506.06811 [astroph.HE] .
 Deaton et al. (2013) M. B. Deaton, M. D. Duez, F. Foucart, E. O’Connor, C. D. Ott, L. E. Kidder, C. D. Muhlberger, M. A. Scheel, and B. Szilágyi, Astrophys. J. 776, 47 (2013), arXiv:1304.3384 [astroph.HE] .
 Foucart et al. (2014) F. Foucart, M. B. Deaton, M. D. Duez, E. O’Connor, C. D. Ott, R. Haas, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 90, 024026 (2014), arXiv:1405.1121 [astroph.HE] .
 Ruffert et al. (1996) M. Ruffert, H.T. Janka, and G. Schaefer, Astron. Astrophys. 311, 532 (1996), astroph/9509006 .
 Rosswog and Liebendörfer (2003) S. Rosswog and M. Liebendörfer, Mon. Not. Roy. Astr. Soc. 342, 673 (2003), arXiv:astroph/0302301 [astroph] .
 O’Connor and Ott (2010) E. O’Connor and C. D. Ott, Class. Quantum Grav. 27, 114103 (2010), arXiv:0912.2393 [astroph.HE] .
 Shibata et al. (2011) M. Shibata, K. Kiuchi, Y. Sekiguchi, and Y. Suwa, Progress of Theoretical Physics 125, 1255 (2011), arXiv:1104.3937 [astroph.HE] .
 Balbus and Hawley (1998) S. A. Balbus and J. F. Hawley, Reviews of Modern Physics 70, 1 (1998).
 Hawley et al. (1995) J. F. Hawley, C. F. Gammie, and S. A. Balbus, Astrophys. J. 440, 742 (1995).
 De Villiers and Hawley (2003) J.P. De Villiers and J. F. Hawley, Astrophys.J. 592, 1060 (2003), arXiv:astroph/0303241 [astroph] .
 McKinney and Gammie (2004) J. C. McKinney and C. F. Gammie, Astrophys. J. 611, 977 (2004), arXiv:astroph/0404512 .
 Chawla et al. (2010) S. Chawla, M. Anderson, M. Besselman, L. Lehner, S. L. Liebling, et al., Phys. Rev. Lett. 105, 111101 (2010), arXiv:1006.2839 [grqc] .
 Etienne et al. (2012) Z. B. Etienne, Y. T. Liu, V. Paschalidis, and S. L. Shapiro, Phys. Rev. D 85, 064029 (2012), arXiv:1112.0568 [astroph.HE] .
 Etienne et al. (2013) Z. B. Etienne, Y. T. Liu, V. Paschalidis, and S. L. Shapiro, ArXiv eprints (2013), arXiv:1303.0837 [astroph.HE] .
 Paschalidis et al. (2015) V. Paschalidis, M. Ruiz, and S. L. Shapiro, Astrophys. J. 806, L14 (2015), arXiv:1410.7392 [astroph.HE] .
 Liu et al. (2008) Y. T. Liu, S. L. Shapiro, Z. B. Etienne, and K. Taniguchi, Phys. Rev. D78, 024012 (2008), arXiv:0803.4193 [astroph] .
 Anderson et al. (2008) M. Anderson et al., Phys. Rev. Lett. 100, 191101 (2008), arXiv:0801.4387 [grqc] .
 Giacomazzo et al. (2011) B. Giacomazzo, L. Rezzolla, and L. Baiotti, Phys. Rev. D 83, 044014 (2011), arXiv:1009.2468 [grqc] .
 Rezzolla et al. (2011) L. Rezzolla, B. Giacomazzo, L. Baiotti, J. Granot, C. Kouveliotou, et al., Astrophys.J. 732, L6 (2011), arXiv:1101.4298 [astroph.HE] .
 Kiuchi et al. (2014) K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and T. Wada, Phys. Rev. D 90, 041502 (2014), arXiv:1407.2660 [astroph.HE] .
 Kiuchi et al. (2015) K. Kiuchi, P. CerdáDurán, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Phys. Rev. D 92, 124034 (2015), arXiv:1509.09205 [astroph.HE] .
 Beckwith et al. (2008) K. Beckwith, J. F. Hawley, and J. H. Krolik, Astrophys. J. 678, 1180 (2008), arXiv:0709.3833 [astroph] .
 Wan (2017) M.B. Wan, Phys. Rev. D95, 104013 (2017), arXiv:1606.09090 [astroph.HE] .
 Sano et al. (2004) T. Sano, S.i. Inutsuka, N. J. Turner, and J. M. Stone, Astrophys. J. 605, 321 (2004), arXiv:astroph/0312480 [astroph] .
 Shi et al. (2016) J.M. Shi, J. M. Stone, and C. X. Huang, Mon. Not. Roy. Astron. Soc. 456, 2273 (2016), arXiv:1512.01106 [astroph.HE] .
 Neilsen et al. (2014) D. Neilsen, S. L. Liebling, M. Anderson, L. Lehner, E. O’Connor, and C. Palenzuela, Phys. Rev. D 89, 104029 (2014), arXiv:1403.3680 [grqc] .
 Palenzuela et al. (2015) C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. L. Caballero, E. O’Connor, and M. Anderson, Phys. Rev. D 92, 044045 (2015), arXiv:1505.01607 [grqc] .
 Janiuk et al. (2007) A. Janiuk, Y. Yuan, R. Perna, and T. Di Matteo, Astrophys. J. 664, 1011 (2007), arXiv:0704.1325 .
 Setiawan et al. (2004) S. Setiawan, M. Ruffert, and H.T. Janka, mnras 352, 753 (2004), astroph/0402481 .
 Fernández and Metzger (2013) R. Fernández and B. D. Metzger, Mon. Not. Roy. Astr. Soc. 435, 502 (2013), arXiv:1304.6720 [astroph.HE] .
 Fernández et al. (2015) R. Fernández, D. Kasen, B. D. Metzger, and E. Quataert, Mon. Not. Roy. Astr. Soc. 446, 750 (2015), arXiv:1409.4426 [astroph.HE] .
 Shibata et al. (2007) M. Shibata, Y. Sekiguchi, and R. Takahashi, Prog. Theor. Phys. 118, 257 (2007), arXiv:0709.1766 .
 Barkov and Baushev (2011) M. V. Barkov and A. N. Baushev, New Astronomy 16, 46 (2011), arXiv:0905.4440 [astroph.HE] .
 Shibata and Sekiguchi (2012) M. Shibata and Y. Sekiguchi, Prog.Theor.Phys. 127, 535 (2012), arXiv:1206.5911 [astroph.HE] .
 Janiuk et al. (2013) A. Janiuk, P. Mioduszewski, and M. Moscibrodzka, Astrophys.J. 776, 105 (2013), arXiv:1308.4823 [astroph.HE] .
 Siegel and Metzger (2017) D. M. Siegel and B. D. Metzger, ArXiv eprints (2017), arXiv:1705.05473 [astroph.HE] .
 Goodman and Xu (1994) J. Goodman and G. Xu, Astrophys. J. 432, 213 (1994).
 Obergaulinger et al. (2009) M. Obergaulinger, P. CerdáDurán, E. Müller, and M. A. Aloy, Astron. Astrophys. 498, 241 (2009).
 Hawley and Balbus (1992) J. F. Hawley and S. A. Balbus, Astrophys. J. 400, 595 (1992).
 Brandenburg (1996) A. Brandenburg, apjlet 465, L115 (1996).
 Lattimer and Swesty (1991) J. M. Lattimer and F. D. Swesty, Nucl. Phys. A535, 331 (1991).
 Perego et al. (2016) A. Perego, R. M. Cabezón, and R. Käppeli, ApJ Suppl. 223, 22 (2016), arXiv:1511.08519 [astroph.IM] .
 Noble et al. (2009) S. C. Noble, J. H. Krolik, and J. F. Hawley, Astrophys. J. 692, 411 (2009), arXiv:0808.3140 [astroph] .
 Piontek and Ostriker (2004) R. A. Piontek and E. C. Ostriker, Astrophys. J. 601, 905 (2004), arXiv:astroph/0310510 [astroph] .
 Pino and Mahajan (2008) J. Pino and S. M. Mahajan, Astrophys. J. 678, 1223 (2008), arXiv:0802.0274 [astroph] .
 Baiotti and Rezzolla (2006) L. Baiotti and L. Rezzolla, Phys. Rev. Lett. 97, 141101 (2006), arXiv:grqc/0608113 [grqc] .
 Faber et al. (2007) J. A. Faber, T. W. Baumgarte, Z. B. Etienne, S. L. Shapiro, and K. Taniguchi, Phys. Rev. D76, 104021 (2007), arXiv:0708.2436 [grqc] .
 Liu et al. (2009) Y. T. Liu, Z. B. Etienne, and S. L. Shapiro, Phys. Rev. D 80, 121503(R) (2009).
 Etienne et al. (2010a) Z. B. Etienne, Y. T. Liu, and S. L. Shapiro, Phys. Rev. D 82, 084031 (2010a), arXiv:1007.2848 [astroph.HE] .
 Muhlberger et al. (2014) C. D. Muhlberger, F. H. Nouri, M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, M. A. Scheel, B. Szilágyi, and S. A. Teukolsky, Phys. Rev. D 90, 104014 (2014), arXiv:1405.2144 [astroph.HE] .
 Thornburg (2004) J. Thornburg, Class. Quant. Grav. 21, 3665 (2004), arXiv:grqc/0404059 [grqc] .
 Schnetter et al. (2006) E. Schnetter, P. Diener, E. N. Dorband, and M. Tiglio, Class. Quant. Grav. 23, S553 (2006), arXiv:grqc/0602104 [grqc] .
 Reisswig et al. (2007) C. Reisswig, N. T. Bishop, C. W. Lai, J. Thornburg, and B. Szilágyi, Class. Quantum Grav. 24, S327 (2007), arXiv:grqc/0610019 .
 Pazos et al. (2009) E. Pazos, M. Tiglio, M. D. Duez, L. E. Kidder, and S. A. Teukolsky, Phys. Rev. D80, 024027 (2009), arXiv:0904.0493 [grqc] .
 Pollney et al. (2011) D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, and P. Diener, Phys. Rev. D 83, 044045 (2011).
 Zink et al. (2008) B. Zink, E. Schnetter, and M. Tiglio, Phys. Rev. D77, 103015 (2008), arXiv:0712.0353 [astroph] .
 Korobkin et al. (2011) O. Korobkin, E. B. Abdikamalov, E. Schnetter, N. Stergioulas, and B. Zink, Phys. Rev. D 83, 043007 (2011), arXiv:1011.3010 [astroph.HE] .
 Korobkin et al. (2012) O. Korobkin, S. Rosswog, A. Arcones, and C. Winteler, Mon. Not. Roy. Astr. Soc. 426, 1940 (2012), arXiv:1206.2379 [astroph.SR] .
 Reisswig et al. (2013) C. Reisswig, R. Haas, C. D. Ott, E. Abdikamalov, P. Mösta, D. Pollney, and E. Schnetter, Phys. Rev. D 87, 064023 (2013), arXiv:1212.1191 [astroph.HE] .
 Koldoba et al. (2002) A. V. Koldoba, M. M. Romanova, G. V. Ustyugova, and R. V. E. Lovelace, Astrophys. J. 576, L53 (2002), arXiv:astroph/0209598 [astroph] .
 Fragile et al. (2009) P. C. Fragile, C. C. Lindner, P. Anninos, and J. D. Salmonson, Astrophys. J. 691, 482 (2009), arXiv:0809.3819 [astroph] .
 Romanova et al. (2012) M. M. Romanova, G. V. Ustyugova, A. V. Koldoba, and R. V. E. Lovelace, Mon. Not. Roy. Astron. Soc. 421, 63 (2012), arXiv:1111.3068 [astroph.SR] .
 Shibata et al. (2006) M. Shibata, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys. Rev. D 74, 104026 (2006), astroph/0610840 .
 Squire and Bhattacharjee (2014) J. Squire and A. Bhattacharjee, Phys. Rev. Lett. 113, 025006 (2014), arXiv:1406.6582 [physics.plasmph] .
 Pessah and Psaltis (2005) M. E. Pessah and D. Psaltis, Astrophys. J. 628, 879 (2005), arXiv:astroph/0406071 [astroph] .
 Foucart et al. (2015) F. Foucart, E. O’Connor, L. Roberts, M. D. Duez, R. Haas, L. E. Kidder, C. D. Ott, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 91, 124021 (2015), arXiv:1502.04146 [astroph.HE] .
 Harten et al. (1983) A. Harten, P. D. Lax, and B. van Leer, SIAM Review 25, 35 (1983).
 Farris et al. (2012) B. D. Farris, R. Gold, V. Paschalidis, Z. B. Etienne, and S. L. Shapiro, Phys. Rev. Lett. 109, 221102 (2012), arXiv:1207.3354 [astroph.HE] .
 Liebling et al. (2010) S. L. Liebling, L. Lehner, D. Neilsen, and C. Palenzuela, Phys. Rev. D81, 124023 (2010), arXiv:1001.0575 [grqc] .
 Penner (2011) A. J. Penner, Mon. Not. Roy. Astr. Soc. 414, 1467 (2011), arXiv:1011.2976 [astroph.HE] .
 Mösta et al. (2014) P. Mösta, B. C. Mundim, J. A. Faber, R. Haas, S. C. Noble, T. Bode, F. Löffler, C. D. Ott, C. Reisswig, and E. Schnetter, Class. Quantum Grav. 31, 015005 (2014), arXiv:1304.5544 [grqc] .
 Mattsson et al. (2004) K. Mattsson, M. Svärd, and J. Nordström, J. Sci. Comput. 21, 57 (2004).
 Sadourny (1972) R. Sadourny, Monthly Weather Review 144, 136 (1972).
 Ronchi et al. (1996) C. Ronchi, R. Iacono, and P. S. Paolucci, J. Comput. Phys. 124, 93 (1996).
 Etienne et al. (2010b) Z. B. Etienne, Y. T. Liu, and S. L. Shapiro, Phys. Rev. D 82, 084031 (2010b), arXiv:1007.2848 .
 Shiokawa et al. (2012) H. Shiokawa, J. C. Dolence, C. F. Gammie, and S. C. Noble, Astrophys. J. 744, 187 (2012), arXiv:1111.0396 [astroph.HE] .
 Balsara (2001) D. Balsara, Astrophys. J. Suppl. Ser. 132, 83 (2001).
 Giacomazzo and Rezzolla (2007) B. Giacomazzo and L. Rezzolla, New frontiers in numerical relativity. Proceedings, International Meeting, NFNR 2006, Potsdam, Germany, July 1721, 2006, Class. Quant. Grav. 24, S235 (2007), arXiv:grqc/0701109 [grqc] .
 Fishbone and Moncrief (1976) L. G. Fishbone and V. Moncrief, Astrophys. J. 207, 962 (1976).
 Balsara and Spicer (1999) D. S. Balsara and D. Spicer, Journal of Computational Physics 148, 133 (1999).
 Ruffert et al. (2001) M. Ruffert, H. T. Ruffert, and H. T. Janka, Astron. Astrophys. 380, 544 (2001), arXiv:astroph/0106229 [astroph] .
 Koldoba et al. (2016) A. V. Koldoba, G. V. Ustyugova, P. S. Lii, M. L. Comins, S. Dyda, M. M. Romanova, and R. V. E. Lovelace, New Astronomy 45, 60 (2016), arXiv:1508.00003 .
 Shapiro and Teukolsky (1983) S. L. Shapiro and S. A. Teukolsky, Research supported by the National Science Foundation. New York, WileyInterscience, 1983, 663 p. (1983).
 van Leer (1977) B. van Leer, J. Comput. Phys. 23, 276 (1977).
 (103) M. G. et al, GNU Scientific Library Reference Manual, 3rd ed.