The Role of Magnetic Field Geometry in the Evolution of Neutron Star Merger Accretion Discs
I.M. Christie, A. Lalakos, A. Tchekhovskoy, R. Fernández, F. Foucart, E. Quataert, D. Kasen
MMNRAS , 1–14 (2019) Preprint 5 July 2019 Compiled using MNRAS L A TEX style file v3.0
The Role of Magnetic Field Geometry in the Evolution of NeutronStar Merger Accretion Discs
I.M. Christie (cid:63) , A. Lalakos † , A. Tchekhovskoy , R. Fern´andez , F. Foucart ,E. Quataert , & D. Kasen , Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA Department of Physics, University of Alberta, Edmonton, AB T6G 2E1, Canada Department of Physics and Astronomy, University of New Hampshire, Durham, NH 03824, USA Departments of Physics & Astronomy, and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720, USA Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Received... / Accepted...
ABSTRACT
Neutron star mergers are unique laboratories of accretion, ejection, and r-process nucleosyn-thesis. We used 3D general relativistic magnetohydrodynamic simulations to study the role ofthe post-merger magnetic geometry in the evolution of merger remnant discs around station-ary Kerr black holes. Our simulations fully capture mass accretion, ejection, and jet produc-tion, owing to their exceptionally long duration exceeding 4 s. Poloidal post-merger magneticfield configurations produce jets with energies E jet ∼ (4 − × erg, isotropic equivalentenergies E iso ∼ (4 − × erg, opening angles θ jet ∼ − ◦ , and durations t j (cid:46) f ej ∼ −
40% of the post-mergerdisc mass, continuing out to times > E jet (cid:39) × erg, E iso ∼ erg, θ jet ∼ . − ◦ , and t j ∼ . (cid:38) f ej (cid:39) ff -axis observers. In comparison toGW 170817 / GRB 170817A, our simulations under-predict the mass of the blue relative to redcomponent by a factor of few. Including the dynamical ejecta and neutrino absorption mayreduce this tension.
Key words: accretion, accretion discs - stars: gamma-ray burst: general - stars: black holes -stars: jets
The recent detection of the neutron star (NS) merger GW170817in both gravitational and electromagnetic (EM) waves has markeda monumental achievement for multi-messenger astronomy (Ab-bott et al. 2017a,b). The first detected component of the EM sig-nal presented itself as a short duration γ -ray burst (GRB), detected ∼ . ∼ (cid:63) E-mail: [email protected] † E-mail: [email protected] tivistic jet launched by a compact object, either a NS or a black hole(BH) (Metzger & Berger 2012; Berger 2014). Moreover, computa-tional studies, using 3D general relativistic magnetohydrodynamic(GRMHD) simulations of BH accretion, in the context of NS merg-ers, have shown that a naturally forming, laterally non-uniform,structured jet can reproduce (Kathirgamaraju et al. 2018) the ob-served radio and X-ray afterglow emission from GW 170817 / GRB170817A (Alexander et al. 2018; Margutti et al. 2018).From this event, it is now widely believed that NS mergersare an important site for r-process nucleosynthesis in the Universe(Kasen et al. 2017; Cˆot´e et al. 2018; Hotokezaka et al. 2018).Confirmation of this fact comes from the photometric and spec-troscopic observations of the kilonova (Cowperthwaite et al. 2017;Chornock et al. 2017; Drout et al. 2017; Tanaka et al. 2017; Tan-vir et al. 2017), argued to be produced from mildly relativistic(i.e. speeds of v ∼ . c ), neutron-rich ejecta radioactively heated c (cid:13) a r X i v : . [ a s t r o - ph . H E ] J u l Christie et al.
Table 1.
Models considered and initial parameters. From left to right: Maximum magnetic field strength within the initial torus, simulation duration t max in seconds and in units of r g / c , plasma β of the initial torus, and simulation resolution in terms of the number of cells (in the radial, poloidal, and toroidaldirections), the e ff ective resolution near the midplane π/ ∆ θ , and cell aspect ratio near the midplane ∆ r : r ∆ θ : r ∆ φ . Note that the latter two are determined at aradial distance of 40 r g .Model Field Max Field Duration, t max Initial Total resolution, E ff ective θ -re- Cell aspect ratio,Name Geometry Strength (G) ( s ) (10 r g / c ) plasma (cid:104) β (cid:105) N r × N θ × N φ solution, π/ ∆ θ ∆ r : r ∆ θ : r ∆ φ BPS Poloidal 1 . × . . × ×
64 256 3 : 1 : 8BPW Poloidal 3 . × . × ×
128 640 8 : 1 : 10BT Toroidal 4 . × . . × ×
128 256 3 : 1 : 4 by r-process elements (Metzger et al. 2010; Roberts et al. 2011;Tanaka 2016; Metzger 2017). The kilonova transient was observedto transition from a blue optical component to an infrared one ina few days (Chornock et al. 2017, see references therin), whichis consistent with theory if considering both blue emission, fromlow-opacity light r-process elements, and red emission, from high-opacity heavy r-process elements.There are two main mechanisms responsible for mass ejec-tion in the kilonova. The first is through dynamical ejecta beingexpelled on ∼ ms timescales by tidal forces (Rosswog et al. 1999;Hotokezaka et al. 2013) or shock interactions (Oechslin et al. 2007;Sekiguchi et al. 2016a). The second mechanism involves outflowsfrom an accretion disc formed from bound merger material. Thisdisc can evolve on longer timescales (i.e. ∼
100 ms − ∼ . v (cid:38) . c ,above the upper limit found in Siegel & Metzger (2017, 2018).In NS mergers, a torus with a primarily toroidal field is ex-pected from the tidal disruption of one (or both) neutron stars andfrom flux freezing (e.g. Kiuchi et al. 2014). Traditionally, GRMHDsimulations of black hole accretion discs have used poloidal fluxloops to initialize the disc, which is known to launch relativisticjets and drive outflows (Blandford & Znajek 1977). Several stud-ies have found that while a purely toroidal seed magnetic field issu ffi cient for the MRI to operate in the disc, such systems pro-duce extremely weak jets (Beckwith et al. 2008; McKinney et al.2012). Recently, Liska et al. (2018) demonstrated, for a geomet-rically thick and radially extended accretion disc, that an initiallytoroidal magnetic field can generate a large-scale poloidal mag-netic flux through what appears to be an α − Ω dynamo (Mo ff att1978) and produce a very powerful jet of power comparable to (or Even though the magnetic field is expected to be toroidally dominated,the geometry may be more complex, containing small-scale orientationflips. even exceeding) the accretion power. In fact, if the results of Liskaet al. (2018) applied to less radially extended post-merger accretiondiscs, they would imply jets that are 4 − ff ect this configuration (e.g. purely poloidal andpurely toroidal geometries) has on the accretion rate, relativisticjets, and the large-scale outflows, including the implications forand connections with the observed kilonova of GW 170817 / GRB170817A. In Sec. 2, we briefly describe the simulation setup. InSec. 3, we present our results for the mass rate and energetics ofall outflows, including the relativistic jet. In Sec. 4, we discuss theconnection of our results with sGRB observations and the observedkilonova of GW 170817 / GRB 170817A while concluding in Sec. 5.
We performed simulations using
HARMPI , an enhanced version ofthe serial open-source code HARM (Gammie et al. 2003; Noble et al.2006), with the addition of several physical processes, e.g. neu-trino cooling and nuclear recombination (for more details, see F19).Throughout, we use spherical polar coordinates r , θ , φ in the Kerr-Schild foliation. For neutrino cooling, we adopt the emission ratesdescribed in Janka (2001) and suppress emission in optically thickregions by a factor of e − τ ν , where τ ν = ρ/ g cm − , (1)and ρ is the gas density. Simulations were initialized with a BH ofmass M BH = M (cid:12) , where M (cid:12) is the solar mass, and spin param-eter a = .
8, surrounded by a torus of mass 0 . M (cid:12) and con-stant initial electron fraction Y e = .
1. We employ an ideal gas lawequation of state (EOS) with a constant adiabatic index γ ad = / T is determined from the total pressure,with contributions from the radiation, electron, proton, and neutron https://github.com/atchekho/harmpi We note that these simulations, with a torus mass of 0 . M (cid:12) , were per-formed before the results of GW 170817 / GRB 170817 A were announcedand observational modeling was performed. Studies have since inferred aninitial torus mass of ∼ . M (cid:12) (Shibata et al. 2017).MNRAS000
1. We employ an ideal gas lawequation of state (EOS) with a constant adiabatic index γ ad = / T is determined from the total pressure,with contributions from the radiation, electron, proton, and neutron https://github.com/atchekho/harmpi We note that these simulations, with a torus mass of 0 . M (cid:12) , were per-formed before the results of GW 170817 / GRB 170817 A were announcedand observational modeling was performed. Studies have since inferred aninitial torus mass of ∼ . M (cid:12) (Shibata et al. 2017).MNRAS000 , 1–14 (2019) agnetic Fields and Neutron Star Merger Discs Table 2.
Summary of our results. From left to right: Cumulative jet energy E jet , cumulative isotropic-equivalent jet energy E iso , jet opening angle (cid:104) θ jet (cid:105) (averaged over both jets and up to 1 s), accreted mass M accr , ejected mass M ejec , ejected mass within the red kilonova component M ejec , red (with electronfraction Y e < .
25) and the blue component M ejec , blue ( Y e > . (cid:104) v r (cid:105) , the average radial speed within the red (cid:104) v r (cid:105) red and blue (cid:104) v r (cid:105) blue kilonova components, and the average electron fraction (cid:104) Y e (cid:105) of all ejecta. All mass values listed as percentages are normalized to the initialtorus mass (0 . M (cid:12) ) while speeds are normalized to the speed of light.Model E jet E iso (cid:104) θ jet (cid:105) M accr M ejec M ejec , red M ejec , blue (cid:104) v r (cid:105) (cid:104) v r (cid:105) red (cid:104) v r (cid:105) blue (cid:104) Y e (cid:105) Name (10 erg) (10 erg) ( ◦ ) (%) (10 − M (cid:12) ) (%) (10 − M (cid:12) ) (%) (10 − M (cid:12) ) (%) (10 − M (cid:12) )BPS 25 22 13 60 2 40 1 . . . .
18 0 .
17 0 . . . . . . .
99 27 0 .
89 3 0 . .
08 0 .
07 0 .
16 0 . . . . . .
89 25 0 .
83 2 0 .
066 0 .
05 0 .
05 0 .
08 0 . components: P = [1 + Y e ] ρ k Tm n + a rad T . (2)Here, a rad is the radiation constant and m n is the neutron mass. Theelectron fraction Y e is evolved according to the numerical proce-dures outlined in F19. We note that our choice for the adiabatic in-dex γ ad was selected by comparing with hydrodynamic simulationswhich use a physical EOS (see Appendix A1 of F19).We performed three simulations di ff ering only in the initialpost-merger magnetic field geometry within the torus. We consid-ered two models, one with a strong poloidal magnetic field config-uration (BPS, described in detail in F19) and one with a weak fieldconfiguration (BPW model). The initial conditions for both mod-els are described by a vector potential A φ ∝ r ρ , which is thenmodified to maximize the magnetic flux in the torus as describedin Tchekhovskoy et al. (2011). For each of the two poloidal con-figurations, we normalized the magnetic field strength such that thedensity-weighted ratio of gas to magnetic pressure within the disc, (cid:104) β (cid:105) ρ = (cid:82) ρ p gas d V (cid:82) ρ p mag d V , (3)is (cid:104) β (cid:105) ρ =
100 for BPS and 850 for BPW, respectively. Hered V = √− g d r d θ d φ is the volume element and g is the determi-nant of the metric. For BPS, the MRI is easily resolved at a mod-erate resolution throughout the torus and yet the magnetic field isnot too strong to violently distort the torus after being amplified bythe shear and the MRI. For BPW, the magnetic field is ∼ φ -direction as inBPS. We provide a summary of each configuration setup, includingthe adopted simulation resolution, in Table 1.The third and final configuration is a toroidal magnetic fieldmodel, denoted as model BT, with plasma β ≡ p gas / p mag = β value because: i) itwas feasible to resolve the MRI given the available computationalresources and ii) the magnetic pressure is low enough so it does notdisrupt the disc. In all simulations, our numerical grid extends fromjust inside the event horizon to ∼ r g in the radial direction andfrom 0 to π in the θ and φ -directions.We carried the simulations out to t max ∼ (3 − × r g / c (cid:39) − r g = GM BH / c is the gravitational radius of the BHand c is the speed of light. Along with the BPS model describedin F19, these are the longest run simulations to date, as measuredin the units of r g / c (e.g. longer than the 2 × r g / c duration inNarayan et al. 2012). This unusually long duration is necessary formass ejection to complete: the cumulative ejected mass dependenceon time flattens out at late times (see Fig. 6(b)). It is also necessaryto capture the jet activity that lasts several seconds (see Fig. 12). We provide a summary of our results in Table 2 and include videosof each simulation in Supplementary Information. Upon the start of the simulation, the disc shear leads to the devel-opment of the MRI, which amplifies the magnetic field and powersmagnetized turbulence in the disc. This drives accretion of gas ontothe black hole. As shown in Fig. 1(a), the mass accretion rate on theblack hole increases and peaks around 10 ms ( ∼ r g / c ). Themass accretion rate peaks slightly earlier for the strong poloidalcase and slightly later for weaker magnetic fields. Following thepeak, ˙ M accr decays in the form of a power-law whose slope is es-sentially independent of the post-merger field geometry. Interest-ingly, the power-law decay portion of ˙ M accr is roughly the same forall configurations, suggesting that the e ff ects of the magnetic fieldgeometry are not important qualitatively for the evolution of the ac-cretion disc past the initial burn-in period (see also Beckwith et al.2008). This decline in the accretion rate comes from the reductionin the mass of the disc, due to both accretion onto the BH and ejec-tion of gas in outflows.We can perform a more quantitative comparison by look-ing at the total amount of material accreted by the BH, M accr , asshown in Fig. 1(b) and Table 2. The amount of accreted mate-rial reaches an asymptotic value by ∼ M accr ∼
60% (0 . M (cid:12) ), followed by ∼
67% (0 . M (cid:12) ) for weak poloidal field model BPW, and ∼ . M (cid:12) ) for toroidal field model BT. Stronger poloidal magneticfields lead to stronger outflows, so there is less gas left to be con-sumed by the BH. Interestingly, the weaker poloidal magnetic fieldmodels accrete approximately the same amount of mass but do notreach the hydrodynamic limit (see F19). The simulated discs can eject energy in the form of outflowslaunched by the magnetic fields twisted by the rotation of the BH(Blandford & Znajek 1977; Komissarov 2001; Tchekhovskoy et al.2010b) or the accretion disc (Blandford & Payne 1982). Typically, https://goo.gl/ct7Htx : There are two sets of videos. The first con-tains two panels, with the left and right panels showing the logarithm ofdensity (in g cm − ) and the electron fraction Y e , respectively, in a verticalslice (see also Fig. 5). The second set displays the mass-weighted red (i.e. Y e < .
25 material) and blue (i.e. Y e < .
25 material) kilonova compo-nents and the jet (green) at a distance of r out = cm ≈ r g (see alsoFig. 10).MNRAS , 1–14 (2019) Christie et al.
Figure 1.
For a wide range of post-merger magnetic field geometries (seelegend), we find a very similar temporal trend in the rest mass accretionrate (panel a) on the BH such that at late times ( t (cid:38) × − s), there isnegligible di ff erence between each model. However, these slight di ff erencesin the temporal decline of ˙ M accr imprint themselves as a small variation inthe amount of material M accr (panel b) accreted on the BH, with a purelytoroidal configuration accreting the most material. The time at which theMRI fully develops coincides with the peak in ˙ M accr . A coloured version ofthis plot is available online. numerical simulations of BH accretion show a combination of thetwo: BH-powered relativistic jets surrounded by sub-relativisticdisc-powered winds (McKinney 2005; Hawley & Krolik 2006;Tchekhovskoy et al. 2011; Tchekhovskoy 2015). We compute thenet sum of these outflow powers through a surface integral P tot ( r ) = − (cid:9) T rt d A , (4)where d A = √− g d θ d φ is the area element, T µν = ( ρ + u + P + b ) u µ u ν + ( P + b ) δ µν − b µ b ν (5)is the stress-energy tensor, b µ is the magnetic four-vector, b = b µ b µ is twice the magnetic pressure, and u µ is the proper velocity (Gam-mie et al. 2003). To distinguish jets and winds, we make use ofthe specific energy flux, µ = − T rt / ( ρ u r ): the value of µ determinesthe maximum possible Lorentz factor an outflow would achieve ifall of its internal and magnetic energy were converted into kineticenergy. We refer to regions with µ ≥ µ < P jet and P wind re-spectively, evaluated at r out . Because relativistic jets are poweredby large-scale poloidal magnetic fields, it is perhaps not surprisingthat the strong poloidal flux model, BPS, forms powerful relativis-tic jets. In fact, Fig. 2(a) shows that the jet power ramps up shortlyafter the light crossing time, r out / c (cid:39) .
033 s, and flattens out at
Figure 2.
The jet power P jet (panel a, solid lines) and jet e ffi ciency η jet (panel b) at r out = cm ≈ r g strongly depend upon the post-mergermagnetic field geometry (see legend in panel b), producing a ∼ timesdi ff erence in P jet between the di ff erent geometries. For all geometries, η jet approaches η MAD ∼ P jet and η jet . During jet disruptions, the EM powerat r out (thin solid line) is therefore solely contained within the surround-ing disc winds. The power contained within the disc winds is ∼ few timesless powerful than the jet and follows a similar temporal trend as P jet . Acoloured version of this plot is available online. P jet (cid:39) × erg s − until ∼ . ffi ciency, η jet ≡ P jet / (cid:104) ˙ M accr c (cid:105) – the ratio of jetto accretion power – from 1% at t ∼ .
05 s to 100% at 0 . ffi ciency by two orders of magnitude impliesthat, unlike a typical expectation that jet power follows mass ac-cretion rate, there is no one-to-one connection between the massaccretion rate and jet power . To understand this, it is helpful tolook at the behavior of the large-scale poloidal magnetic flux thatpasses through the BH and powers the relativistic jets, Φ BH = . (cid:90) r = r H | B r | d A , (6)where the integral is over both hemispheres of the event horizon, r H = r g [1 + √ − a ], and the factor of 0 . P jet ∝ Φ , (7)the constancy of the jet power would imply the constancy of Φ BH MNRAS000
05 s to 100% at 0 . ffi ciency by two orders of magnitude impliesthat, unlike a typical expectation that jet power follows mass ac-cretion rate, there is no one-to-one connection between the massaccretion rate and jet power . To understand this, it is helpful tolook at the behavior of the large-scale poloidal magnetic flux thatpasses through the BH and powers the relativistic jets, Φ BH = . (cid:90) r = r H | B r | d A , (6)where the integral is over both hemispheres of the event horizon, r H = r g [1 + √ − a ], and the factor of 0 . P jet ∝ Φ , (7)the constancy of the jet power would imply the constancy of Φ BH MNRAS000 , 1–14 (2019) agnetic Fields and Neutron Star Merger Discs Figure 3.
The physical Φ BH (panel a) and normalized φ BH ≡ Φ BH / ( ˙ M accr r g c ) / (panel b) poloidal magnetic flux on the BH, poweringthe relativistic jets, varies significantly with the post-merger field geometry(see legend in panel a). For all configurations, φ BH eventually reaches orapproaches a critical value of φ MAD ∼
50, demonstrating that in the contextof short GRBs, large-scale poloidal magnetic flux on the BH can becomedynamically-important and lead to the development of a MAD (see eqn. 8).A coloured version of this plot is available online. (Tchekhovskoy & Giannios 2015). Fig. 3(a) shows that within ∼
10 ms after the start of the simulation, the central BH receivesmost of the large-scale magnetic flux available after the merger, af-ter which the BH flux indeed remains approximately constant (towithin a factor of 2) until t ∼ . P mag ∝ Φ ,scales linearly with the disc mass that in turn is proportional to massaccretion rate, ˙ M accr . As a result, in a MAD, the BH magnetic fluxis proportional to the square root of the mass accretion rate, Φ BH (cid:39)
50 ( ˙ M accr r g c ) / , (8)where we have included the proportionality factor (Tchekhovskoy Figure 4.
The time-dependence of normalized poloidal magnetic flux on thenorthern hemisphere of the BH ( φ BH , north , yellow) shows flips in the mag-netic polarity of jets (see also Fig. 5). Namely, the polarity on the northernhemisphere of the BH flips sign several times throughout the simulation,resulting in the production of current sheets that propagate along the jets.Note that the total absolute magnetic flux exceeds that through the northernhemisphere, implying the presence of non-equatorial current sheet(s) on theBH event horizon. A colour version of this plot is available online. φ BH = Φ BH ( ˙ M accr r g c ) / , (9)is (cid:39)
50. Equations (7) and (8) imply that in the MAD state, P jet (cid:39) . a ˙ M accr c , (10)where we again have provided the proportionality factor(Tchekhovskoy 2015). That this numerical factor for our rapidlyspinning BH with a = . ffi ciency η jet , as seen in Fig. 2(b). Thus, at late times, the jetpower approximately equals the accretion power, and both decayas a power-law in time with a slope of ∼ − . M accr shown in Fig. 1.Our weaker post-merger poloidal field BPW model leads to anorder of magnitude weaker jet than in the BPS model: P jet initiallypeaks at 3 × erg s − and eventually plateaus at (cid:39) erg s − .As in BPS, this flattening is due to the constancy of the large-scalepoloidal magnetic flux on the BH, as seen in Fig. 3(a) between t (cid:39) . M crit (cid:39) Φ / (50 r g c ), at which,the magnetic flux becomes dynamically-important, leading to theformation of a MAD. Beyond this point, the jet power follows themass accretion rate, P jet ∼ ˙ M accr c , as is typical of MADs and seenin Fig. 2(b). Models BPS and BPW demonstrate that compact ac-cretion discs (as typical for short GRBs) can naturally reach a MADstate (see Proga & Zhang 2006 for 2D analog). Because the jetpower tracks the mass accretion rate, we are conveniently providedwith an inside view of the late time accretion on the BH.Surprisingly, we find that the purely toroidal post-merger ge-ometry BT model also launches jets of substantial power, P jet (cid:39) erg s − at t (cid:39) . t ∼ . − ff , with typical P jet ∼ erg s − . Fig. 2(b) shows that the corresponding jet e ffi ciencygradually increases in time, eventually approaching 100%. We dis-cuss this in more detail in Sec. 3.3. MNRAS , 1–14 (2019)
Christie et al.
Figure 5.
A sequence of vertical slices through density in our toroidal post-merger field simulation BT showing flips in magnetic polarity (solid and dottedlines representing positive and negative polarity, respectively) of the jets (blue regions). The magnetic polarity switches from negative (left panel) to positive(right panel), causing the BH magnetic flux to reconnect away, dropping the power in the relativistic jets (see Fig. 2), and the surrounding disc winds to chokethe jets (middle panel; see also Fig. 2(a)). This is the first demonstration of a striped jet formation in a BH accretion simulation with the stripes naturallyemerging through a large-scale poloidal flux dynamo in the disc rather than introduced through an initial condition. A colour version of this plot is availableonline.
How is it possible that even in the absence of any poloidal magneticflux, model BT produces relativistic jets, which require large-scalepoloidal magnetic flux (Beckwith et al. 2008; McKinney & Bland-ford 2009; McKinney et al. 2012)? This has never been seen in nu-merical simulations of compact discs, such as those expected in abinary merger. The yellow line in Fig. 4 shows the time-dependenceof the magnetic flux through the northern hemisphere of the BH, φ BH , north , normalized by the mass accretion rate (see eqn. 9). Whileinitially starting at zero, the magnetic flux increases in magnitude.This suggests that that the initial purely toroidal post-merger mag-netic field undergoes a dynamo-like process that can generate large-scale poloidal magnetic flux. We can see the dynamo action inthe movies of model BT through the emergence of poloidal mag-netic loops above and below the equatorial plane. This behavior isconsistent with an α − Ω like poloidal flux dynamo (Mo ff att 1978;Liska et al. 2018). The dimensionless flux evolves slowly with time,approaching the critical MAD state (although φ BH ∼
50 is notstrictly reached in the duration of the simulation). The polarity ofthe emerging poloidal magnetic flux appears to switch at random,likely reflecting the randomness of the magnetized turbulence un-derlying the dynamo.To understand the connection between the magnetic flux andthe jet power, it is helpful to look at the absolute magnetic flux, φ BH (see eqn. 9), the square of which controls the jet power. We seethat the jet power varies significantly (see Fig. 2(a)), with the jetsshutting o ff at multiple times. The jets become suppressed due tothe winds of the surrounding disc choking and disrupting the jets,especially at the times when the magnetic flux vanishes and the jetsare weakened. Fig. 5 illustrates how one such flux flip happens.The jets seen in blue in the left panel (at t ≈ .
03 s), have a well-defined structure. In contrast, in the middle panel (at t ≈ .
16 s), thesurrounding winds disrupt the jets. The right panel shows that even-tually the jets manage to push through (at t ≈ .
25 s). Such mag-netic flux polarity flips occur frequently and do not appear to show https://goo.gl/ct7Htx any obvious periodicity, as seen in Fig. 4. The average duration be-tween the flips appears to increase with the increasing simulationtime. We do not see such sign flips in the BPS and BPW models,suggesting that at least on large scales, the initial post-merger mag-netic flux dominates in these models the flux, if any, produced bythe dynamo. Note that jet power shut-o ff s occur more frequentlythan the polarity flips, indicating that many factors (not just thestrength of the magnetic flux but also, e.g., mass-loading of the po-lar regions by the ambient gas) determine the success of relativisticjet formation. When the jets are shut o ff , the EM power at r out issolely contained within the surrounding disc winds. This power,displayed as the thin blue line in Fig. 2(a), can contribute substan-tially to the total EM power of the combined jet + wind regions.At late times, t (cid:38) φ BH approachesthe critical value of 50, and Fig. 2(b) shows that the dimensionlessjet power η jet approaches the critical value of 1. This suggests thateven absent poloidal post-merger magnetic flux, the system man-ages to generate its own poloidal flux and approach (but not quitereach) the MAD state. With a longer simulation, it is plausible thatwe would see a full MAD state develop in our BT model, allowingus to use the jet power as an observational window in the accretionon the BH (see eqn. 10).Why has no simulation seen the development of strong BHmagnetic fields in simulations of compact discs with purely toroidalinitial magnetic field? There are several possible explanations.First, BH magnetic flux becomes substantial (Fig. 4) and the jetsbecome noticeably strong relative to the accretion flow (Fig. 2) onlyat t (cid:38) t ∼ × r g / c ,much longer than the typical simulation duration of ∼ r g / c .Thus, previous simulations might not have been long enough toobserve this e ff ect. Second, in order to see the dynamo action, weneeded to use very high resolutions, 512 × ×
128 cells (see Ta-ble 1). We found that while a simulation at twice as small resolution(i.e.256 × ×
64 cells) marginally resolved the toroidal MRI, itdid not resolve the poloidal MRI, and did not show noticeable signsof large-scale poloidal magnetic flux dynamo.We note that in the context of radially-extended accretion
MNRAS000
MNRAS000 , 1–14 (2019) agnetic Fields and Neutron Star Merger Discs Figure 6.
The rest mass outflow rate’s ˙ M out (panel a), evaluated at r out = cm ≈ r g , strong dependence upon the post-merger field geometry(see legend) demonstrates that strong poloidal flux is required for launch-ing prompt mass outflows. The time in which the initial outflow reaches r out depends on the post-merger geometry, stemming from the developing andsaturation time of the MRI as well as the velocity of the outflows’ strongdependence on the post-merger field geometry (see Table 2). The large vari-ation in ˙ M out for times (cid:46) ∼
12% ( ∼ × − M (cid:12) )di ff erence in the amount of ejected material M ejec (panel b). At late times( t (cid:38) M out becomes insensitive to the post-merger geometry, displayinga temporal trend of ∝ t − . in all models, resulting in a flattening of M ejec .A coloured version of this plot is available online. discs, Liska et al. (2018) found the operation of large-scale poloidalflux dynamo and the formation of strong jets, with η jet (cid:38)
1. This iscomparable to the jet e ffi ciency we find, but only at very late times.In fact, it takes our simulations three times longer than those ofLiska et al. (2018) to reach η jet ∼
1. Why is this so? If the largeradial extent of the disc is a prerequisite for the dynamo to operatee ffi ciently and produce powerful jets, it would take our small disca substantial amount of time until it appreciably expands radially.Importantly, unlike Liska et al. (2018), our jets also show polarityflips. We discuss potential reasons for this di ff erence in Sec. 4.2. Mass outflows and their composition are particularly important asthey determine the luminosity, color, and duration of the kilonova(see Secs. 3.5 and 4.4). We quantify the ejecta by measuring themass outflow rate ˙ M out through a sphere of radius r out = cm ≈ r g . This is su ffi ciently far from the BH to avoid the interac-tions with the turbulent and “viscously” expanding accretion disc. As shown in Fig. 6(a), the outflows reach r out earliest for strong The disc eventually does expands out to such large radii, however, by thattime it has very low density and carries little mass. post-merger poloidal magnetic fields, model BPS, followed at latertimes by weak poloidal, BPW, and purely toroidal, BT, models.This time di ff erence results from not only the MRI reaching sat-uration earlier for stronger poloidal magnetic fields, but also fromstronger poloidal fields launching faster outflows, as seen in Ta-ble 2. Namely, the average radial velocity of the ejecta for theBPS model is (cid:104) v r (cid:105) ∼ . c , much higher than ∼ . c for BPWand ∼ . c for BT.The amount of ejected material also varies by model, as shownin Fig. 6(b). For instance, in the BPS model, the mass outflow raterapidly ramps up in a fraction of a second and plateaus at ˙ M out ∼ − M (cid:12) s − . In contrast, in weak poloidal, BPW, and toroidal, BT,models the outflows remain about an order of magnitude weakerand catch up to the BPS model only by the end of the first sec-ond. This implies that strong post-merger poloidal magnetic flux isconducive to launching prompt mass outflows. Interestingly, massoutflows are much more similar at late times, past the first sec-ond: the outflow rate in all 3 models largely decays as a power-law,˙ M out ∝ t − . , suggesting that qualitatively the ability of post-mergersystems to launch outflows at late times becomes insensitive to thepost-merger magnetic field geometry. The above di ff erences, pri-marily the prompt mass ejection, lead in the strong poloidal fieldBPS model to an overall largest mass ejection, carrying 40% of theinitial torus mass (0 . M (cid:12) ). This is a third more than the ∼ . M (cid:12) ) for BPW and ∼
27% (0 . M (cid:12) ) for BTmodels.To isolate the e ff ects of post-merger magnetic fields, F19 com-pared the strong poloidal field BPS model to an otherwise identi-cal hydrodynamic model. They found that the strong poloidal mag-netic fields in the BPS model ejected about twice as much mass asin the hydrodynamic model, primarily because the hydrodynamicmodel was missing the mass ejection during the first second afterthe merger. Even our torodial BT model ejects more mass (about athird more) than the hydrodynamic models of F19. In Secs. 3.2 and 3.4, we analyzed the energetics and mass of ouroutflows. However, these outflows are expected to consist of mate-rial with a range of compositions with spatially varying Y e values.Fig. 7 shows, at di ff erent times, the breakdown of ejecta mass M out into bins of Y e . Because the post-merger torus is initialized with anelectron fraction of Y e = .
1, it is not surprising that at early times (cid:46) . (i.e Y e ≤ . Y e ∼ .
1. However, the amount of materialcentered at Y e ∼ . ∼ . − r out for all The average radial velocities, determined by (cid:104) v r (cid:105) = (cid:8) (( ρ + u + P ) u r u r + Pg rr ) d A / (cid:8) ρ u r d A , are slightly higher than those found in F19 due to av-eraging over momentum rather than density. Nuclear reaction network calculations show that this critical value of Y e ∼ .
25 separates the point at which no lanthanides are formed (Lippuner &Roberts 2015). These elements are key for the opacity and hence, the colorof the kilonova (Kasen et al. 2013).MNRAS , 1–14 (2019)
Christie et al. Y e d M [ M ] BPS [0,0.1] s [0,0.5] s [0,1] s [0,2] s [0,4] s [0,9.2] s Y e d M [ M ] BPW [0,0.1] s [0,0.5] s [0,1] s [0,2] s [0,4] s [0,4.4] s Y e d M [ M ] BT [0,0.1] s [0,0.5] s [0,1] s [0,2] s [0,4] s [0,4.3] s Figure 7.
Histograms of the cumulative amount of ejected material passingthrough r out = cm ≈ r g (in addition, see Fig. 6) vs electron frac-tion Y e , for each post-merger geometry. Increasing the post-merger poloidalfield strength not only ejects more material at earlier times, but also spreadsthis material over a broader range of Y e values, producing a more extendedlanthanide-poor (i.e. Y e ≥ .
25) region, influenced by the increasing impor-tance of positron capture. A colour version of this plot is available online. post-merger geometries with a spread in Y e extending to larger val-ues (cid:38) .
25. These regions form close to the BH at times t (cid:46) . ff ectivelyincreasing Y e . We find that the mean electron fraction within allejecta is almost independent of the post-merger geometry, with (cid:104) Y e (cid:105) being 0 .
16 for BPS, 0 .
19 for BPW, and 0 .
18 for BT (see Table 2). -2.0 -1.5 -1.0 -0.5 0.0 log ( v r / c ) d M [ M ] Y e -2.0 -1.5 -1.0 -0.5 0.0 log ( v r / c ) d M [ M ] Y e BPSBPWBT
Figure 8.
Histograms of the cumulative amount of ejected material passingthrough r out vs the logarithm of the radial velocity (normalized to c andequally spaced in 60 bins ranging from − Y e ≤ .
25, top panel) and lanthanide-poor ( Y e > .
25, bottom panel) regions (seealso Fig. 6). Low Y e -material is characterized by slower radial velocitieswhich are spread over a large range of obtainable v r values whereas low Y e -material has faster velocities confined within a narrower range of v r . Asdisplayed in Fig. 6, weaker or more toroidal post-merger field geometrieseject less material and at lower velocities. A coloured version of this plot isavailable online. As the interpretation of the lanthanide-rich and poor regionshas direct applicability to the observed kilonova (see Sec. 4.4 fordetails), it is helpful to investigate the physical and geometricalproperties of the two regions. In Fig. 9, we show how the massoutflow rate ˙ M out , averaged over the φ -direction, of the lanthanide-rich and poor matter is spread over time t and angle θ . Within thelanthanide-rich regions, the amount of material traversing through r out follows a similar trend presented in Figs. 6 and 7, with a major-ity of the material reaching r out at earlier times in the BPS modeland later for the BPW and BT models. For all geometries, thelanthanide-rich material is concentrated in the regions between therelativistic jet and the equatorial plane, with the peak amount ofmaterial being associated with the peak in ˙ M out seen in Fig. 6(a).Although the post-merger field geometry governs the total amountof ejected material, it weakly influences the ejected fraction (nor-malized to the initial torus mass ) of lanthanide-rich gas, ∼
90% ofthe total ejected material (see Table 2). The material passes through r out with mildly relativistic speeds; namely v r ∼ . − . c , as dis-played in the top panel of Fig. 8. The time averaged radial velocityof the lanthanide-rich regions is (cid:104) v r (cid:105) red ∼ . c for the BPS model,much faster than ∼ . c for the BPW model and ∼ . c for theBT model. Even larger speeds are found in material moving withinthe jet, however, it is orders of magnitude less dense than the sur-rounding winds (see Fig. 5 for comparison).For the lanthanide-poor regions, the material begins to passthrough r out at ∼ . ∼ . ∼ MNRAS000
90% ofthe total ejected material (see Table 2). The material passes through r out with mildly relativistic speeds; namely v r ∼ . − . c , as dis-played in the top panel of Fig. 8. The time averaged radial velocityof the lanthanide-rich regions is (cid:104) v r (cid:105) red ∼ . c for the BPS model,much faster than ∼ . c for the BPW model and ∼ . c for theBT model. Even larger speeds are found in material moving withinthe jet, however, it is orders of magnitude less dense than the sur-rounding winds (see Fig. 5 for comparison).For the lanthanide-poor regions, the material begins to passthrough r out at ∼ . ∼ . ∼ MNRAS000 , 1–14 (2019) agnetic Fields and Neutron Star Merger Discs c o s M eject =37% Y e BPS M out [ M / s ] 0.0 0.5 1.0 1.5 2.0-101 c o s M eject =3% Y e BPS M out [ M / s ]0 1 2 3 4-101 c o s M eject =27% BPW M out [ M / s ] 0.0 0.5 1.0 1.5 2.0-101 c o s M eject =3% BPW M out [ M / s ]0 1 2 3 4 t [ s ]-101 c o s M eject =25% BT M out [ M / s ] 0.0 0.5 1.0 1.5 2.0 t [ s ]-101 c o s M eject =2% BT M out [ M / s ] Figure 9.
A space-time diagram of ˙ M out (integrated over the φ -direction), as seen on a sphere of radius r out = cm ≈ r g as a function of cos θ and time t for all post-merger field geometries (please see the colour bars). The three rows, from top to bottom, show the results for the BPS, BPW, and BT models,respectively. The left panels show lanthanide-rich material ( Y e ≤ . Y e ≥ . . M (cid:12) . We find that low- Y e material ( Y e ≤ .
25) is ejected in large polar regions, extending from the equatorial plane of the disc to regions closeto the jet. Higher- Y e material ( Y e ≥ . https://goo.gl/ct7Htx ). A coloured version of this plot is available online. important. The amount of ejected material contained within thelanthanide-poor region is ∼
3% (i.e. ∼ − M (cid:12) ) within all mod-els, independent of the post-merger geometry (see Table 2). Thematerial passing through r out is contained within a much narrowerangular region than the lanthanide-rich material, with polar widthof ∆ θ ∼ ◦ − ◦ . We illustrate this in Fig. 10, by plotting ona sphere of radius r out for each post-merger geometry a temporalsnapshot (at t ∼ . , in addition to the relativistic jets (green).At early times, when the lanthanide-poor ejecta initially crosses r out , it is concentrated near the poles. For later times, the lanthanide-poor gas appears to emerge at larger polar angles (closer to theequatorial plane). This gas, initially obscured by the lanthanide-rich material, passes through r out at relativistic speeds larger thanthe lanthanide-rich material, namely v r (cid:38) . c with velocitiesreaching up to c , as presented in the bottom panel of Fig. 8. Thetime averaged radial velocity of the lanthanide-poor material is (cid:104) v r (cid:105) blue ∼ . c for the BPS model, ∼ . c for the BPW model, and ∼ . c for the BT model, eventually punching through and over- For videos displaying the full time evolution of Fig. 10, see: https://goo.gl/ct7Htx . This obscuration occurs very early in the disc evolution and should notbe mistaken with obscuration occurring within the late-time (i.e. t (cid:29) taking the lanthanide-rich material. This ejection is mildly asym-metric relative to the equatorial plane for the BPS and BT models,with more material being ejected through the southern hemisphere.This anisotropy manifests itself as a ∼ × − M (cid:12) di ff erence in thelanthanide-poor region of our BPS model and a ∼ . × − M (cid:12) di ff erence in the lathanide-rich region of the BT model. We have found that for a range of magnetic field geometries,post-merger accretion discs can naturally develop dynamically-important BH magnetic flux that turns them into MADs, providingus with an inside view of the mass accretion rate on the BH fromthe jet power. Namely, at early times, the jet power is set by theamount of large-scale poloidal magnetic flux present in the accre-tion flow, whereas at late times it is set by the mass accretion rate.This outcome is insensitive to the post-merger magnetic field ge-ometry which suggests that it is a robust phenomenon, even thoughit has not been found previously. Until now, smaller torii extendingout to (cid:46) r g embedded with a single poloidal field loop, have beenfound to lead to Standard And Normal Evolution (SANE, Narayanet al. 2012) discs, in which the gas pressure dominates the magnetic MNRAS , 1–14 (2019) Christie et al.
Figure 10.
Temporal snapshots ( t ≈ . r out = cm ≈ r g . The colorintensity of the blue and red kilonova components represent the mass-weighted rest-mass outflow rate ˙ M out (see Figs. 6 and 9) while the jet (green) intensityis weighted by its power (see Fig. 2). The blue component is confined within a much narrower region (polar width ∆ θ ∼ ◦ − ◦ ) than the red component.Initially obscured by the red component, the blue component moves faster, punching through and eventually overtakes the red component. The relativistic jetwithin each model is tightly collimated by the surrounding disc winds, obtaining opening angles θ jet (cid:46) ◦ (see also Fig. 11). For videos displaying the fulltime evolution, see https://goo.gl/ct7Htx . A coloured version of this plot is available online. pressure in the disc, (e.g., Gammie et al. 2003; McKinney 2005;Hawley & Krolik 2006). Indeed to obtain MADs, researchers havepreviously opted for very large accretion discs, much larger thanthose typically expected in the context of compact binary mergers:the large size of the discs allowed them to contain a large enoughpoloidal magnetic flux to flood the BH (Tchekhovskoy 2015; Haw-ley et al. 2015). For instance, to obtain MADs, Tchekhovskoy et al.(2011) considered discs of a large extent ∼ × r g while McKin-ney et al. (2012) simulated even larger, unbound discs that extendedout to infinity. Narayan et al. (2012) obtained similar results whilefocusing on discs that extended out to 10 r g .While such large discs are naturally expected in active galacticnuclei (AGN), binary mergers lead to much smaller discs. How cana small disc in a binary merger turn MAD? Instead of starting witha large amount of magnetic flux, a disc can start with little flux andevolve to the point when very little of the initial gas is left. At thislate time, what was initially a weak and subdominant magnetic fluxcan become a dynamically-important one. In fact, this is a naturalway of producing MADs in any system whose mass accretion ratedecreases over time, such as tidal disruption events (Tchekhovskoyet al. 2014) and core-collapse gamma-ray bursts (Tchekhovskoy& Giannios 2015). In this work, we demonstrated that MADs cannaturally develop in 3D numerical simulations of initially small ac-cretion discs, as typical for binary mergers (see Proga & Zhang2006 for a 2D analog). However, in order for a MAD state to occur,an unusually long evolution time (by GRMHD simulation durationstandards) is required, which might explain why this e ff ect has notbeen previously seen.In our strong poloidal magnetic field simulation, model BPS,the magnetic fields became dynamically-important around t MAD ≈ . φ BH = φ MAD (cid:39)
50 (see Fig. 3(b)). In aweaker poloidal magnetic field model, BPW, this happens at a fewtimes later time, t MAD ≈ φ BH ≈
35, more than half way tothe critical φ MAD value, by the end of the simulation. This suggests that given a longer duration (e.g. ∼
10 s), the BH magnetic flux canbecome dynamically-important even for this purely toroidal post-merger magnetic flux geometry. We discuss how this can happen inSec. 4.2.
Poloidal magnetic flux is a crucial prerequisite for jet formation.Indeed, it is the winding of the poloidal magnetic field by the BH(Blandford & Znajek 1977) or the inner regions of the disc (Bland-ford & Payne 1982) that is typically associated with magnetically-powered outflows. However, the shear between two merging neu-tron stars is expected to amplify the toroidal magnetic field com-ponent, naturally leading to a toroidally-dominated field geometry(the field direction might undergo polarity flips on small scales dueto the Kelvin-Helmholtz instability). In the absence of su ffi cientlystrong poloidal magnetic fields, how do binary mergers manage toproduce jets at all?Our toroidal field simulation, model BT, might shed lighton this long-standing problem. We find that as the simu-lation progresses, the accretion flow spontaneously developspoloidal magnetic field loops above and below the equatorialplane (see the movies in the Supplementary Information and athttps: // goo.gl / ct7Htx). This behaviour is similar to that of an α − Ω poloidal flux dynamo (Mo ff att 1978), in which buoyancy and Cori-olis forces work together to convert toroidal flux into poloidal mag-netic flux loops. Because of the chaotic nature of the dynamo, themagnetic flux polarity varies randomly from one loop to another.Once the newly formed magnetic flux loops reach the BH, theypower jets of alternating magnetic flux polarity, or striped jets.At large distances in the jet, magnetic reconnection in the currentsheets separating the regions of opposite polarity can provide nat-ural dissipation sites responsible for high-energy GRB jet emis-sion (Giannios & Uzdensky 2018). Note that the total poloidal fluxon the BH event horizon, determined by eqn. (6), is greater thanthe sum of the magnetic fluxes from its individual northern andsouthern hemispheric components. This non-zero di ff erence be-tween Φ BH and | Φ BH , north | + | Φ BH , south | emerges due to the presenceof the current sheets near the BH horizon.Recently, Liska et al. (2018) found that radially-extended ac- MNRAS000
Poloidal magnetic flux is a crucial prerequisite for jet formation.Indeed, it is the winding of the poloidal magnetic field by the BH(Blandford & Znajek 1977) or the inner regions of the disc (Bland-ford & Payne 1982) that is typically associated with magnetically-powered outflows. However, the shear between two merging neu-tron stars is expected to amplify the toroidal magnetic field com-ponent, naturally leading to a toroidally-dominated field geometry(the field direction might undergo polarity flips on small scales dueto the Kelvin-Helmholtz instability). In the absence of su ffi cientlystrong poloidal magnetic fields, how do binary mergers manage toproduce jets at all?Our toroidal field simulation, model BT, might shed lighton this long-standing problem. We find that as the simu-lation progresses, the accretion flow spontaneously developspoloidal magnetic field loops above and below the equatorialplane (see the movies in the Supplementary Information and athttps: // goo.gl / ct7Htx). This behaviour is similar to that of an α − Ω poloidal flux dynamo (Mo ff att 1978), in which buoyancy and Cori-olis forces work together to convert toroidal flux into poloidal mag-netic flux loops. Because of the chaotic nature of the dynamo, themagnetic flux polarity varies randomly from one loop to another.Once the newly formed magnetic flux loops reach the BH, theypower jets of alternating magnetic flux polarity, or striped jets.At large distances in the jet, magnetic reconnection in the currentsheets separating the regions of opposite polarity can provide nat-ural dissipation sites responsible for high-energy GRB jet emis-sion (Giannios & Uzdensky 2018). Note that the total poloidal fluxon the BH event horizon, determined by eqn. (6), is greater thanthe sum of the magnetic fluxes from its individual northern andsouthern hemispheric components. This non-zero di ff erence be-tween Φ BH and | Φ BH , north | + | Φ BH , south | emerges due to the presenceof the current sheets near the BH horizon.Recently, Liska et al. (2018) found that radially-extended ac- MNRAS000 , 1–14 (2019) agnetic Fields and Neutron Star Merger Discs Figure 11.
The opening angles of both jets (see the legend for explana-tion of line types) substantially di ff er between purely toroidal and purelypoloidal post-merger magnetic geometries (see legend), with the former,model BT, having smaller opening angles. In the BT model, the surround-ing disc winds disrupt the jet (i.e. θ jet = ◦ , see also Fig. 2), forcing thejet to punch through the material, producing a tightly collimated angle. Therange of opening angles in our simulations, 5 − ◦ is roughly consistentwith the range of inferred opening angles in short GRBs (Fong et al. 2015).A coloured version of this plot is available online. cretion discs initially threaded with a purely toroidal magnetic fieldcan produce powerful jets. Similar to our work, they see the signsof α − Ω dynamo and the formation in the accretion disc of poloidalfield loops of alternating polarity. However, they find that most ofthe loops become ejected in an outflow, and a single magnetic loopends up dominating the jet energetics and magnetic flux polarity.Why are our results di ff erent than theirs? There can be two possi-ble factors that can suppress outflows in our work and encouragethe disc to retain the poloidal flux loops instead of ejecting them.First, our smaller disc is more tightly bound and is therefore lessconducive to outflows. Second, our neutrino cooling makes the disceven more tightly bound, resulting in additional outflow suppres-sion. Both of these e ff ects might encourage the disc to retain mostof the freshly generated poloidal magnetic flux loops and encouragethem to accrete on the BH, leading to alternating BH magnetic fluxand striped jets. In future work, we will investigate the robustnessof this phenomenon and its relationship to the size of the accretiondisc, the presence of radiative or neutrino cooling, and the sensitiv-ity of the results to the numerical resolution. X-ray and optical afterglow observations of GRBs exhibiting asteepening in their temporal decline are often characterized by jetbreaks (Soderberg et al. 2006; Nicuesa Guelbenzu et al. 2011; Fonget al. 2012, 2014). The time associated with these breaks can un-cover characteristic properties of the jet, such as the jet openingangle θ jet (Sari et al. 1999; Frail et al. 2001). Although only a fewshort GRBs exhibit such a break, Fong et al. (2015) were able to es-timate their opening angles to span a wide range from ∼ (cid:104) θ jet (cid:105) ≈ ◦ ± ◦ .These observations raise an important question: what is the jetcollimating agent in short GRBs? In long-duration GRBs, whichcan also be tightly collimated (into opening angles as small as afew degrees, Frail et al. 2001; Cenko et al. 2010), a natural collimat-ing agent is the radially-extended stellar envelope of size (cid:38) r g (Tchekhovskoy et al. 2010a). It does not appear plausible, however, Figure 12.
Not only does the post-merger field geometry drastically e ff ectthe jet power (panel a), but it also results in a large di ff erence in the cumu-lative jet energy (panel b) and its isotropic equivalent (panel c). For the BPSmodel, a majority of the energy passes through r out within the first ∼ ∼ that a compact post-merger remnant disc extending out to ∼ r g can manage to collimate short GRB jets into the smallest observedshort GRB opening angles of ∼ ◦ . What is the way out of this co-nundrum? Our simulations reveal that magnetized turbulence leadsto angular momentum transport and viscous-like spreading of thedisc from the initial size of ∼ r g to (cid:38) r g . While this is an orderof magnitude smaller than the size of the stellar envelope in longGRBs, outflows launched from such an extended disc could sub-stantially collimate the jets. Does this lead to su ffi ciently small jetopening angles, θ jet , that span the range of short GRB observations?Figure 11 shows that at early times ( t ∼ . θ jet reaches its peak, the largest being θ jet ∼ ◦ for the BPS model, followed by ∼ ◦ for the BPW model, and ∼ ◦ for the BT model. These di ff erences in the opening angle MNRAS , 1–14 (2019) Christie et al. could result from: i) the lack of su ffi cient material surrounding thejet at early times, required to tightly collimate it, as the jet mate-rial launched around the polar axis moves faster than the surround-ing disc winds , and / or ii) the power contained within the jet is (cid:38)
50 times larger for the BPS model than the other models (seeFig. 2). For late times, t (cid:38) . r out at larger polar angles, pro-ducing a tight collimation of the jets. For BPS and BPW models,the time-average opening angles (averaged over their the jet activityperiod, i.e. t (cid:46) (cid:104) θ jet (cid:105) ∼ ◦ and ∼ . ◦ , respectively (see Table 2). For the purely toroidal BTmodel, the intermittence of θ jet follows that of P jet (see Fig. 2): thesurrounding disc winds disrupt the jets leading to θ jet = ◦ . As thejets reform near the BH, they have to drill through the disrupted ma-terial, resulting in a tighter collimation with a time-average value(averaged over 1 s) of (cid:104) θ jet (cid:105) ∼ . ◦ .The investigation of the jet opening angle’s dependence on thepost-merger geometry has important implications on the inferredenergies of the resulting afterglow. It is di ffi cult to measure directlythe intrinsic jet energy, E jet . More easily accessible is its isotropicequivalent energy, E iso , which is related to E jet by the beaming fac-tor f b ≡ − (cos θ jet , north + cos θ jet , south ) /
2, such that E iso = E jet / f b . Al-though the inferred value of E iso from GRB afterglows is model de-pendent (e.g. assumed particle spectrum, radiative e ffi ciency, den-sity of external ambient medium), a compilation of 38 short GRBsreports a median value of E iso ∼ × erg while the distributionspreads over a range of ∼ × erg to ∼ erg (Fong et al.2015). Assuming the median jet opening angle reported above (i.e. ∼ ◦ , Fong et al. 2015), this corresponds to a characteristic in-ferred intrinsic jet energy of E jet ∼ erg.We show in Figs. 12(b) and (c) the cumulative jet energyand its isotropic equivalent, as measured at a sphere of radius r out = cm ≈ r g for all three post-merger magnetic ge-ometries. For the BPS model, the large spike in P jet within thefirst ∼ E jet quickly reaching an asymptotic value of ∼ . × erg. At late times (cid:38) ∼ r out , beyond which it slowly reachesan asymptotic value of E iso ∼ . × erg. For weaker post-merger poloidal field geometries, such as our BPW model, the jetpower remains roughly constant throughout time and does not be-gin to decrease until (cid:38) E jet and E iso , this corresponds toa continuous increase in time, eventually leveling o ff at values of E jet ∼ . × erg and E iso ∼ × erg, respectively, at ∼ E jet ∼ × erg and E iso ∼ . × erg, respectively. However, it is important to notethe late time (i.e. (cid:38) . ff erence between E jet and E iso . Thereis a continuous increase in the former up until the end of the sim-ulation due to the late time increase in the jet power. However, atthis time, θ jet also increases significantly above its time averagedvalue (see Fig. 11). The increase in both quantities results in a late We note that in our simulations, we do not consider neutrino / anti-neutrino annihilation. Such e ff ects can deposit enough energy into the polarregions to drive mildly relativistic outflows, clearing the poles of baryons(Fujibayashi et al. 2017; Foucart et al. 2018). The distinction between the jet and disc winds is made by performing acut on the specific energy flux: µ = − T rt / ( ρ u r ) ≥ time decrease in the isotropic equivalent of the jet power presentingitself as a leveling o ff of E iso .Note that in our simulations the sub-relativistic winds sur-rounding the jets can carry a fraction of the jet energy in all models:typically, E wind ∼ × erg for BPS, ∼ erg for BPW, and ∼ × erg for BT models. This energy will eventually becomevisible at the forward shock, where the ejecta runs into the ambientmedium and produces the afterglow emission. Months-long radioand X-ray afterglow seen from GW170817 is thought to be pow-ered by the relativistic jet (Margutti et al. 2018; Kathirgamarajuet al. 2018; Alexander et al. 2018). A detailed analysis of the inferred properties of GW 170817 / GRB170817A has been reported in (Kasen et al. 2017; Kasliwal et al.2017; Kilpatrick et al. 2017). The total amount of ejected materialis estimated within the range of ∼ . − . M (cid:12) . To match the op-tical and infrared observations (Arcavi et al. 2017; Cowperthwaiteet al. 2017; Drout et al. 2017), Kasen et al. (2017) modeled thered (i.e. Y e < .
25) and blue (i.e. Y e > .
25) kilonova componentswith ejected masses M red ≈ . M (cid:12) and M blue ≈ . M (cid:12) and ve-locities of v red ≈ . c and v blue ≈ . c , respectively. Their analysissuggests that the mechanism for mass ejection is predominantly viaoutflows from a remnant accretion disc. Although our simulationsare consistent with this interpretation, there are several di ff erencesfound, for each post-merger field geometry, when making a com-parison to observations, as discussed below.From our results, we find that the total ejected mass (0 . M (cid:12) for BPS, 0 . M (cid:12) for BPW, and 0 . M (cid:12) for BT; see Fig. 6 andTable 2) is lower than the values listed above. To obtain valuesconsistent with observational modeling, we would require an ini-tial torus mass (cid:38) . M (cid:12) , if we were to simply rescale our results.As shown in the left column of Fig. 9, a majority of the ejectedmass has Y e ≤ .
25. As compared with observations, all configu-rations underpredict the ejected mass of the red component by (cid:46) ∼ . M (cid:12) , the amount of material within the red compo-nent (for all post-merger geometries) would be consistent with theinferred values from observational modeling of GW 170817 / GRB170817A. However, the average radial velocity of the red compo-nent varies significantly with the post-merger geometry, with theBT model obtaining smaller velocities than what is inferred (seeTable 2).For the blue kilonova component, all three models underpre-dict the inferred mass by several orders of magnitude ( ∼ − M (cid:12) from all models as compared to 0 . M (cid:12) from observational mod-eling). If our results were simply rescaled with the torus mass,we would require (cid:38) . M (cid:12) , much larger than what is expectedpost-merger. Moreover, it is useful to note that there is a di ff er-ence in the geometrical interpretation of the blue kilonova com-ponent. When modeling the emission from the blue component,Kasen et al. (2017) took the blue material to be painted over aspherical region within polar angles 0 ◦ ≤ θ ≤ ◦ . As discussedin Sec. 3.5 (see Figs. 9 and 10) for all three models, the blue mate-rial at r out = cm ≈ r g is confined within narrower regionsof polar width ∆ θ ∼ ◦ − ◦ . It could be, however, that our simu-lations may not have reached a free-expansion phase. If so, it couldbe likely that geometry of the blue region would be modified beforebeing observed as a kilonova.The aforementioned analysis therefore suggest that an initiallytoroidal field, for our idealized initial conditions, struggles to re- MNRAS000
25. As compared with observations, all configu-rations underpredict the ejected mass of the red component by (cid:46) ∼ . M (cid:12) , the amount of material within the red compo-nent (for all post-merger geometries) would be consistent with theinferred values from observational modeling of GW 170817 / GRB170817A. However, the average radial velocity of the red compo-nent varies significantly with the post-merger geometry, with theBT model obtaining smaller velocities than what is inferred (seeTable 2).For the blue kilonova component, all three models underpre-dict the inferred mass by several orders of magnitude ( ∼ − M (cid:12) from all models as compared to 0 . M (cid:12) from observational mod-eling). If our results were simply rescaled with the torus mass,we would require (cid:38) . M (cid:12) , much larger than what is expectedpost-merger. Moreover, it is useful to note that there is a di ff er-ence in the geometrical interpretation of the blue kilonova com-ponent. When modeling the emission from the blue component,Kasen et al. (2017) took the blue material to be painted over aspherical region within polar angles 0 ◦ ≤ θ ≤ ◦ . As discussedin Sec. 3.5 (see Figs. 9 and 10) for all three models, the blue mate-rial at r out = cm ≈ r g is confined within narrower regionsof polar width ∆ θ ∼ ◦ − ◦ . It could be, however, that our simu-lations may not have reached a free-expansion phase. If so, it couldbe likely that geometry of the blue region would be modified beforebeing observed as a kilonova.The aforementioned analysis therefore suggest that an initiallytoroidal field, for our idealized initial conditions, struggles to re- MNRAS000 , 1–14 (2019) agnetic Fields and Neutron Star Merger Discs produce the kilonova properties (e.g. mass and velocities) inferredfrom observational modeling. It should be noted, however, that sim-ply rescaling the mass of the initial torus would not be entirelyaccurate as discs with larger masses are more opaque to neutri-nos, requiring a more elaborate treatment of neutrino cooling thanwhat was used here (F19). Moreover for GW 170817, Shibata et al.(2017) estimated the range of potential torus masses to lie within0 . − . M (cid:12) , putting a constraint on our analysis. In order toprovide a more complete model, the inclusion of the post-mergerdynamical ejecta and a better neutrino transport scheme would benecessary, which could reduce the tension in the amount of massejected within the blue component. Here, we explored the role of the post-merger magnetic field config-uration on the long term disc evolution in the context of NS merg-ers. Beginning with either a purely poloidal or purely toroidal mag-netic field within the torus (see Table 1 for model initial setup), wefind the formation of a relativistic jet whose total power, energy,and opening angle are consistent with typical values inferred fromGRBs (see Figs. 2b, 11, and 12). For all three post-merger magneticfield configurations, we find that the jet power eventually reaches alevel of the order of ˙ M accr c , signifying that the disc has reached aMAD state (see eqn. 10 and Fig. 3b) and that at late times, the jetpower directly tracks the mass accretion rate on the BH. At earliertimes, the jet power is roughly constant and not strongly correlatedwith ˙ M BH , reflecting the large scale poloidal magnetic flux contentin the accretion flow.For the purely toroidal post-merger magnetic field configura-tion, we find the formation of a jet with energetics consistent withGRBs (Fig. 12). A dynamo-like process in the accretion disc leadsto the formation of alternating magnetic flux, shown in Fig. 4, thatpowers striped jets. If this result holds at higher resolution, the pro-duction of current sheets and their reconnection in the jets couldpower the prompt emission in GRBs (e.g. Spruit et al. 2001; Gian-nios & Spruit 2006; Beniamini et al. 2018).Concurrent with the launching of a jet, mass outflows are ex-pelled as disc winds. The driving mechanism of the winds is notfully understood but is most likely a combination of thermal andmagnetic e ff ects as well as small contributions from α -particle re-combination (Metzger et al. 2008; Siegel & Metzger 2018). The to-tal amount of ejected material contained within winds (see Fig. 6) isfound to be smaller than what is inferred from observational mod-eling of GW 170817 / GRB 170817A. However, our simulations as-sumed an initial torus mass of 0 . M (cid:12) . The initial torus mass, de-termined by the masses of the binary components and the assumedequation of state, is expected to be larger than what was assumedhere (Shibata et al. 2017).It is important to note several limitations of our analysis. Thefirst is our choice for the initial electron fraction of Y e = . / NS systems (Foucartet al. 2015, 2017), which have reported an electron fraction rangingfrom ∼ . − .
2. The second is the lack of neutrino absorption andtransport within our models. Recent results of GRMHD simula-tions combined with a full Monte Carlo neutrino transport methodhave shown that ∼
20% of the early time outflows are blue (Milleret al. 2019). As such, these approximations can lead to underesti-mating Y e at early times. Therefore, the fraction of the lanthanide- poor ejecta presented here (see Table 2) can be considered as alower bound.Our choice for the initial torus mass of 0 . M (cid:12) wasmade such that the torus remains optically thin and our ap-proximations for neutrino cooling are reasonable. Moreover, oursimulations results were completed before the observations ofGW170817 / GRB170817A. The results of using a larger torus, andits comparison with the observed kilonova are currently being ex-plored and will be presented elsewhere. We note, however, that formore massive discs, a more complete treatment of neutrino trans-port will be necessary (as in, e.g. Miller et al. 2019).Although the total amount of material contained within thewinds varies with the initial magnetic field configuration, shownin Fig. 6, its composition remains roughly fixed, with ∼ ∼
10% beinglanthanide-poor (see Figs. 7 and 9). A simple rescaling of our re-sults to an initial torus mass of 0 . M (cid:12) leads to lanthanide-richejection consistent with the red kilonova component inferred forGW170817, and with mass ejection falling short by a factor of ∼ ACKNOWLEDGEMENTS
IMC thanks Dr. K. Alexander, Dr. W. Fong, and Dr. B. Met-zger for their supportive discussions. RF acknowledges supportfrom the Natural Sciences and Engineering Research Council ofCanada (NSERC) through Discovery Grant RGPIN-2017-04286,and from the Faculty of Science at the University of Alberta. Thiswork was supported in part by a Simons Investigator award fromthe Simons Foundation (EQ) and the Gordon and Betty MooreFoundation through Grant GBMF5076. This was was supportedby NASA through grant 80NSSC18K0565 (FF, AT). This researchused resources of the National Energy Research Scientific Comput-ing Center (NERSC), which is supported by the O ffi ce of Scienceof the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The software used in this work was in part devel-oped by the DOE NNSA-ASC OASCR Flash Center at the Univer-sity of Chicago. Computations were performed at Carver, Hopper,and Edison (repositories m1186, m2058, m2401, and the scavengerqueue). REFERENCES
Abbott B. P., et al., 2017a, ApJ, 848, L12Abbott B. P., et al., 2017b, ApJ, 848, L13Alexander K. D., et al., 2018, ApJ, 863, L18Arcavi I., et al., 2017, Nature, 551, 64Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Beckwith K., Hawley J. F., Krolik J. H., 2008, ApJ, 678, 1180Beniamini P., Barniol Duran R., Giannios D., 2018, MNRAS, 476, 1785Berger E., 2014, ARA&A, 52, 43Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Blinnikov S. I., Novikov I. D., Perevodchikova T. V., Polnarev A. G., 1984,Soviet Astronomy Letters, 10, 177MNRAS , 1–14 (2019) Christie et al.
Cenko S. B., et al., 2010, ApJ, 711, 641Chornock R., et al., 2017, ApJ, 848, L19Cˆot´e B., et al., 2018, ApJ, 855, 99Cowperthwaite P. S., et al., 2017, ApJ, 848, L17Drout M. R., et al., 2017, Science, 358, 1570Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126Fern´andez R., Metzger B. D., 2013, MNRAS, 435, 502Fern´andez R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019,MNRAS, 482, 3373Fong W., et al., 2012, ApJ, 756, 189Fong W., et al., 2014, ApJ, 780, 118Fong W., Berger E., Margutti R., Zauderer B. A., 2015, ApJ, 815, 102Foucart F., et al., 2015, Phys. Rev. D, 91, 124021Foucart F., et al., 2016a, Phys. Rev. D, 93, 044019Foucart F., O’Connor E., Roberts L., Kidder L. E., Pfei ff er H. P., ScheelM. A., 2016b, Phys. Rev. D, 94, 123016Foucart F., et al., 2017, Classical and Quantum Gravity, 34, 044002Foucart F., Duez M. D., Kidder L. E., Nguyen R., Pfei ff er H. P., ScheelM. A., 2018, Phys. Rev. D, 98, 063007Frail D. A., et al., 2001, ApJ, 562, L55Fujibayashi S., Sekiguchi Y., Kiuchi K., Shibata M., 2017, ApJ, 846, 114Gammie C. F., McKinney J. C., T´oth G., 2003, ApJ, 589, 444Giannios D., Spruit H. C., 2006, A&A, 450, 887Giannios D., Uzdensky D. A., 2018, preprint, ( arXiv:1805.09343 )Goldstein A., et al., 2017, ApJ, 848, L14Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103Hawley J. F., Fendt C., Hardcastle M., Nokhrina E., Tchekhovskoy A.,2015, Space Sci. Rev., 191, 441Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., ShibataM., Taniguchi K., 2013, Phys. Rev. D, 87, 024001Hotokezaka K., Beniamini P., Piran T., 2018, International Journal of Mod-ern Physics D, 27, 1842005Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042Janka H. T., 2001, A&A, 368, 527Just O., Bauswein A., Ardevol Pulpillo R., Goriely S., Janka H.-T., 2015,MNRAS, 448, 541Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Na-ture, 551, 80Kasliwal M. M., et al., 2017, Science, 358, 1559Kathirgamaraju A., Tchekhovskoy A., Giannios D., Barniol Duran R.,2018, preprint, ( arXiv:1809.05099 )Kilpatrick C. D., et al., 2017, Science, 358, 1583Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., Wada T., 2014, Phys.Rev. D, 90, 041502Komissarov S. S., 2001, MNRAS, 326, L41Lippuner J., Roberts L. F., 2015, ApJ, 815, 82Liska M. T. P., Tchekhovskoy A., Quataert E., 2018, preprint,( arXiv:1809.04608 )Margutti R., et al., 2018, ApJ, 856, L18McKinney J. C., 2005, ApJ, 630, L5McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423,3083Metzger B. D., 2017, Living Reviews in Relativity, 20, 3Metzger B. D., Berger E., 2012, ApJ, 746, 48Metzger B. D., Piro A. L., Quataert E., 2008, MNRAS, 390, 781Metzger B. D., et al., 2010, MNRAS, 406, 2650Miller J. M., et al., 2019, arXiv e-prints, p. arXiv:1905.07477Misner C. W., Thorne K. S., Wheeler J. A., 1973, GravitationMo ff att H. K., 1978, Magnetic field generation in electrically conductingfluidsNarayan R., Paczynski B., Piran T., 1992, ApJ, 395, L83Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69Narayan R., Sa¸dowski A., Penna R. F., Kulkarni A. K., 2012, MNRAS, 426,3241Nicuesa Guelbenzu A., et al., 2011, A&A, 531, L6 Noble S. C., Gammie C. F., McKinney J. C., Del Zanna L., 2006, ApJ, 641,626Oechslin R., Janka H.-T., Marek A., 2007, A&A, 467, 395Paczynski B., 1986, ApJ, 308, L43Proga D., Zhang B., 2006, MNRAS, 370, L61Roberts L. F., Kasen D., Lee W. H., Ramirez-Ruiz E., 2011, ApJ, 736, L21Rosswog S., Liebend¨orfer M., Thielemann F.-K., Davies M. B., Benz W.,Piran T., 1999, A&A, 341, 499Sari R., Piran T., Halpern J. P., 1999, ApJ, 519, L17Savchenko V., et al., 2017, ApJ, 848, L15Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., Taniguchi K., 2016a,Phys. Rev. D, 93, 124046Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., Taniguchi K., 2016b,Phys. Rev. D, 93, 124046Shibata M., Fujibayashi S., Hotokezaka K., Kiuchi K., Kyutoku K.,Sekiguchi Y., Tanaka M., 2017, Phys. Rev. D, 96, 123012Siegel D. M., Metzger B. D., 2017, Physical Review Letters, 119, 231102Siegel D. M., Metzger B. D., 2018, ApJ, 858, 52Soderberg A. M., et al., 2006, ApJ, 650, 261Spruit H. C., Daigne F., Drenkhahn G., 2001, A&A, 369, 694Tanaka M., 2016, Advances in Astronomy, 2016, 634197Tanaka M., et al., 2017, PASJ, 69, 102Tanvir N. R., et al., 2017, ApJ, 848, L27Tchekhovskoy A., 2015, in Contopoulos I., Gabuzda D., Kylafis N., eds,Astrophysics and Space Science Library Vol. 414, The Formation andDisruption of Black Hole Jets. p. 45, doi:10.1007 /000