Long-term evolution of a merger-remnant neutron star in general relativistic magnetohydrodynamics I: Effect of magnetic winding
LLong-term evolution of a merger-remnant neutron star in general relativisticmagnetohydrodynamics I: Effect of magnetic winding
Masaru Shibata,
1, 2
Sho Fujibayashi, and Yuichiro Sekiguchi
3, 2 Max Planck Institute for Gravitational Physics (Albert Einstein Institute),Am M¨uhlenberg 1, Potsdam-Golm 14476, Germany Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, 606-8502, Japan Department of Physics, Toho University, Funabashi, Chiba 274-8510, Japan (Dated: February 3, 2021)Long-term ideal and resistive magnetohydrodynamics (MHD) simulations in full general relativityare performed for a massive neutron star formed as a remnant of binary neutron star mergers.Neutrino radiation transport effects are taken into account as in our previous papers. The simulationis performed in axial symmetry and without considering dynamo effects as a first step. In the idealMHD, the differential rotation of the remnant neutron star amplifies the magnetic-field strengthby the winding in the presence of a seed poloidal field until the electromagnetic energy reaches ∼
10% of the rotational kinetic energy, E kin , of the neutron star. The timescale until the maximumelectromagnetic energy is reached depends on the initial magnetic-field strength and it is ∼ ∼ G. After a significant amplification ofthe magnetic-field strength by the winding, the magnetic braking enforces the initially differentiallyrotating state approximately to a rigidly rotating state. In the presence of the resistivity, theamplification is continued only for the resistive timescale, and if the maximum electromagneticenergy reached is smaller than ∼
3% of E kin , the initial differential rotation state is approximatelypreserved. In the present context, the post-merger mass ejection is induced primarily by the neutrinoirradiation/heating and the magnetic winding effect plays only a minor role for the mass ejection. PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg
I. INTRODUCTION
Theoretical exploration for the evolution of the mergerremnants of binary neutron stars has become a hot topicin this decade [1], because such systems can shine as anelectromagnetic counterpart of gravitational waves emit-ted in their inspiral stage and the signals bring us a va-riety of information for the nature of the merger processand neutron stars that cannot be obtained only by thegravitational-wave observation, as the first observationof the binary neutron star merger shows [2, 3]. Thereare a variety of the possibilities for the remnants of thebinary neutron star mergers [1] and for the correspond-ing electromagnetic counterparts [4, 5]. Irrespective ofthe possibilities, the key for the strong electromagneticemission is the mass ejection from the remnant. Thus,one important aspect of the theoretical study is to clar-ify the properties of the ejecta such as total mass, typicalvelocity, and typical elements.The promising component of the remnant that inducesthe mass ejection is a disk (or torus) surrounding thecentral compact object, which is either a massive neu-tron star or a black hole. The reason for this is that thedisk is differentially rotating, approximately with the Ke-plerian rotational profile, and is likely to be magnetizedbecause the disk matter stems from neutron stars. In thissituation, the magnetorotational instability (MRI) [6] isactivated, and as a result, a turbulence would be devel-oped, enhancing the turbulence viscosity. The result-ing viscous heating and angular momentum transportinevitably enhance the activity of the disk, and even- tually, the mass ejection is induced [7–18]. In addition,the amplified magnetic field could eventually constitute aaligned poloidal magnetic field along the rotational axisand this could further enhance the mass ejection effi-ciency [10, 14, 15]. With these considerations, numericalsimulations have been extensively performed in the lastdecade to clarify the details of the mass ejection from theaccretion disks.Although many works have been done for exploringthe evolution of the accretion disks formed as remnantsof the neutron-star merger, there are only a few worksfor directly evolving the remnant massive neutron starfor a long timescale of (cid:38) (cid:36) ), where (cid:36) is the cylindrical radius: Reflecting the nearly irrota-tional velocity field of the pre-merger stage of two neu-tron stars, the inner region of the remnant neutron starsshould be differentially rotating with the angular velocitythat increases with the cylindrical radius ( d Ω /d(cid:36) > d Ω /d(cid:36) >
0, the MRI does not play anyrole.However, irrespective of the sign of d Ω /d(cid:36) , the wind- a r X i v : . [ a s t r o - ph . H E ] F e b ing of the magnetic-field lines occurs in the presence ofa seed poloidal magnetic field together with differentialrotation. In contrast to the MRI for which the growth ofthe magnetic-field strength exponentially proceeds, the(toroidal) magnetic-field strength increases only linearlywith time by the winding. Nevertheless, the growth rateis not still negligible in the presence of the rapid rotationbecause the amplification factor is approximately writ-ten as Ω t (e.g., Refs. [21–23]) and Ω ∼ / s (seeSec. III) for the merger remnants. Thus in 0 . ∼ times unlessstrong resistive processes are present. That is, even inthe presence of a seed poloidal magnetic field typically of10 –10 G (which may be a conservative value for themerger remnants [24]), the (toroidal) magnetic field couldbe amplified to the order of 10 –10 G in ∼ . c = 1 = G where c and G denote the speed of light and the gravitational constant, respectively (but c is often recovered to clarify the unitsin the following sections). Latin and Greek indices de-note the spacetime and space components, respectively.In Sec. II, we suppose to use Cartesian coordinates forthe spatial components whenever equations are written. II. BASIC EQUATIONS AND METHOD FORNUMERICAL COMPUTATIONSA. Brief summary
In this work, we numerically solve Einstein’s equation,ideal or resistive MHD equations, evolution equations forthe lepton fractions including the electron fraction, and(approximate) neutrino-radiation transfer equations. Ex-cept for the MHD parts (see the following subsectionsfor details), the numerical implementation is the same asthat in our latest viscous-hydrodynamics work [16, 17]:Einstein’s equation is solved using the original ver-sion of the Baumgarte-Shapiro-Shibata-Nakamura for-malism [28] together with the puncture formulation [29],Z4c constraint propagation prescription [30], and 5th-order Kreiss-Oliger dissipation. The axial symmetry forthe geometric variables is imposed using the cartoonmethod [31, 32] with the 4th-order accuracy in space.For evolving the lepton fractions, we take into accountelectron and positron capture, electron-positron pair an-nihilation, nucleon-nucleon bremstrahlung, and plasmondecay [11, 17]. We employ a tabulated equation of statebased on the DD2 equation of state [33] for a relativelyhigh-density part and the Timmes (Helmhoitz) equationof state for the low-density part [34]. In this tabulatedequation of state, thermodynamics quantities such as ε , P , and h are written as functions of ρ , Y e , and T where ε , P , h (= c + ε + P/ρ ), ρ , Y e and T are the specific in-ternal energy, pressure, specific enthalpy, rest-mass den-sity, electron fraction, and matter temperature, respec-tively. We choose the lowest rest-mass density to be0 . / cm in the table, and the atmosphere density for ρ ∗ := ρu t √− g in the hydrodynamics simulation is cho-sen to be 10 g / cm in the central region of r ≤
100 kmand it is decreased down to 1 g / cm with the dependenceof ∝ r − in the outer region (i.e., for the far region it is1 g / cm ). Here u µ and g denote the four velocity of thefluid and the determinant of the spacetime metric g µν , re-spectively. ρ ∗ obeys the continuity equation in the formof ∂ t ρ ∗ + ∂ i ( ρ ∗ v i ) = 0 , (2.1)where v i = u i /u t is the three velocity, dx i /dt . FromEq. (2.1), the conserved rest mass is defined by M ∗ := (cid:90) ρ ∗ d x. (2.2) B. Maxwell’s equations
First we write down the equations for the electromag-netic fields, which are derived from Maxwell’s equations: ∇ µ F µν = − πj ν , (2.3) ∇ µ ∗ F µν = 0 . (2.4)Here ∇ µ is the covariant derivative with the respect to g µν , F µν is the electromagnetic tensor, j µ is the currentfour vector, and ∗ F µν is the dual of F µν defined by ∗ F µν := 12 (cid:15) µναβ F αβ , (2.5)with (cid:15) µναβ the Levi-Civita tensor of (cid:15) txyz = √− g . Usingthe unit timelike vector normal to the spacelike hyper-surface, n µ , we define the electric and magnetic fields by E µ := F µν n ν , (2.6) B µ := 12 n α (cid:15) αµνβ F νβ , (2.7)and thus, the electromagnetic tensor and its dual arewritten as F µν = n µ E ν − n ν E µ + (cid:15) µνα B α , (2.8) ∗ F µν = (cid:15) µνα E α − n µ B ν + n ν B µ , (2.9)where (cid:15) ναβ = n µ (cid:15) µναβ . Since E µ n µ = B µ n µ = 0, E t = B t = 0. Note the definition of n µ = − α ∇ µ t with α thelapse function.Then Maxwell’s equations are written in terms of E µ and B µ as follows: For the constraint equations, D k E k = 4 πρ e , (2.10) D k B k = 0 , (2.11)and for the evolution equations, ∂ t E i − L β E i = αKE i − D k (cid:0) α(cid:15) kij B j (cid:1) − πα ¯ j i , (2.12) ∂ t B i − L β B i = αKB i + D k (cid:0) α(cid:15) kij E j (cid:1) , (2.13)where ρ e := − j a n a and ¯ j i := γ ik j k with γ ij the spatialmetric defined by γ µν := g µν + n µ n ν . D k denotes thecovariant derivative with respect to γ ij , β i the shift vec-tor, K the trace of the extrinsic curvature K ij , and L β denotes the Lie derivative with respect to β i . We numer-ically solve the evolution equations in the forms [27] ∂ t E i = − ∂ k (cid:0) β i E k − β k E i + α(cid:15) kij B j (cid:1) − π (cid:0) J i − Qβ i (cid:1) , (2.14) ∂ t B i = − ∂ k (cid:0) β i B k − β k B i − α(cid:15) kij E j (cid:1) , (2.15)where E i := √ γ E i , B i := √ γ B i , J i := √− g ¯ j i , and Q := √ γ ρ e with γ = det( γ ij ) (in curvilinear coordi-nates, the definition should be appropriately modified byexcluding the contribution of the coordinates in γ ). Wenote ( − g ) = α γ . In this notation, the constraint equa-tions are written in the simple forms as ∂ i E i = 4 πQ, (2.16) ∂ i B i = 0 . (2.17) To close the equations, we in general need the Ohm’slaw, for which we here write as j µ = ˜ ρ e u µ + σ c ( F µν u ν + α d ∗ F µν u ν ) , (2.18)where ˜ ρ e := − j µ u µ = w − ( ρ e − σ c E µ u µ + σ c α d B µ u µ )is the charge density observed in the frame comovingwith the fluid, and σ c is the conductivity with w := − n µ u µ = αu t . Note that the resistivity is defined by η := 1 / (4 πσ c ). The third term in the right-hand side ofEq. (2.18) denotes a dynamo term (in the simplest ver-sion) with α d being the so-called α -dynamo parameter(see, e.g., Ref. [35] for the dynamo theory and Ref. [36]for the relativistic formulation). The dynamo term is aphenomenological one and we do not consider it in thispaper (i.e., α d = 0) except for Appendix A 5. Note that F µν u ν and − ∗ F µν u ν (=: b µ ) denote the electric and mag-netic fields in the frame comoving with the fluid.In the ideal MHD case for which we suppose σ c → ∞ ,it is convenient to write the electromagnetic tensor interms of the magnetic field in the frame comoving withthe fluid, b µ , as F µν = (cid:15) µναβ u α b β , (2.19) ∗ F µν = b µ u ν − b ν u µ . (2.20)Then, it becomes trivial that the corresponding electricfield, F µν u ν , vanishes. Then, the condition of F µν u ν = 0with Eq. (2.8) gives the electric field, E i , as E i = − w (cid:15) ijk u j B k , (2.21)and thus, E µ u µ (= E i u i ) = 0.In the ideal MHD case, the current is not determinedfrom Eq. (2.18) due to the fact that σ c → ∞ . In-stead, it has to be determined from Eq. (2.12) for E i given by Eq. (2.21). Thus, only Eq. (2.13) or (2.15)becomes the evolution equation for the electromagneticfield, which determines the magnetic-field evolution. Us-ing Eq. (2.21), this equation is written in the well-knownform as ∂ t B i = ∂ k (cid:0) B k v i − B i v k (cid:1) . (2.22) C. Ideal MHD
In the ideal MHD, the energy-momentum tensor iswritten as T µν = (cid:18) ρh + b π (cid:19) u µ u ν + (cid:18) P + b π (cid:19) g µν − b µ b ν , (2.23)where b = b µ b µ . Here, we have the relations as αb t = B k u k , (2.24) wb i = B i + B k u k u i , (2.25) w b = B + ( B k u k ) , (2.26)with B = B k B k , and thus, the energy-momentum ten-sor is also written in terms of B k .For evolving the MHD equations, we define S := √ γ T µν n µ n ν and S i := −√ γ T µν n µ γ νi , which are writ-ten as S = ρ ∗ hw − √ γP + √ γ π (cid:20) B − w (cid:8) B + ( B k u k ) (cid:9)(cid:21) , (2.27) S i = ρ ∗ hu i + √ γ π (cid:0) B u i − B k u k B i (cid:1) . (2.28)Their evolution equations are ∂ t S + ∂ j (cid:104) S v j + √ γP tot ( v j + β j ) − √− g πw B k u k B j (cid:105) = √− gK ij S ij − S k D k α, (2.29) ∂ t S i + ∂ j (cid:104) S i v j + √− gP tot δ ji − √− g πw B j ( B i + u i B k u k ) (cid:105) = − S ∂ i α + S k ∂ i β k − √ γS jk ∂ i γ jk , (2.30)where P tot := P + b / π and S ij := γ µi γ νj T µν .Multiplying B i to Eq. (2.28), we obtain B i S i = ρ ∗ hB i u i . Thus, Eq. (2.28) is rewritten as S i + √ γ πρ ∗ h B k S k B i = (cid:18) ρ ∗ h + √ γ π B (cid:19) u i . (2.31)Using this, the normalization relation of u µ is written interms of γ ij , ρ ∗ , S i , and B i as γ ij (cid:18) S i + √ γ πρ ∗ h B k S k B i (cid:19) (cid:18) S j + √ γ πρ ∗ h B k S k B j (cid:19) × (cid:18) ρ ∗ h + √ γ π B (cid:19) − + 1 = w . (2.32)Hence, for given (evolved) values of γ ij , ρ ∗ , S i , S , and B i together with a given equation of state, this equationtogether with Eq. (2.27) constitute simultaneous equa-tions for h and w , and thus, they are used for the so-called primitive recovery procedure in the ideal MHD.Our numerical approach for the primitive recovery whenemploying tabulated equations of state is essentially thesame as that in the viscous hydrodynamics case [25]. D. MHD in general cases
In the general MHD case, not only the magnetic fieldbut also the electric field are explicitly included in theenergy-momentum tensor as T µν = ρhu µ u ν + P g µν + 14 π (cid:34) E + B γ µν + n µ n ν ) − E µ E ν − B µ B ν + ( n µ (cid:15) ναβ + n ν (cid:15) µαβ ) E α B β (cid:35) , (2.33) where E = E k E k . We note that in the ideal MHDcase, using Eq. (2.21) together with the fact that E i / (cid:112) E k E k , B i / (cid:112) B k B k , and γ ij u j / (cid:112) γ kl u k u l can con-stitute the orthonormal bases of the three space, theenergy-momentum tensor of Eq. (2.23) is derived in astraightforward manner.For the general MHD case, S and S i are written as S = ρ ∗ hw − √ γP + 18 π √ γ (cid:0) E + B (cid:1) , (2.34) S i = ρ ∗ hu i + √ γ π (cid:15) ijk E j B k , (2.35)and their evolution equations are ∂ t S + ∂ j (cid:34) S v j + √ γ (cid:18) P − E + B π (cid:19) ( v j + β j )+ 14 π √− g(cid:15) jkl E k B l (cid:35) = √− gK ij S ij − S k D k α, (2.36) ∂ t S i + ∂ j (cid:34) S i v j + √− g (cid:18) P + E + B π (cid:19) δ ji − √− g π (cid:16) E i E j + B i B j (cid:17) − √ γ π ( v j + β j ) (cid:15) ikl E k B l (cid:35) = − S ∂ i α + S k ∂ i β k − √− gS jk ∂ i γ jk , (2.37)where S ij = ρhu i u j + P γ ij + 14 π (cid:20) − E i E j − B i B j + 12 γ ij ( E + B ) (cid:21) . (2.38)Using Eq. (2.35), the normalization relation for u µ iswritten as γ il (cid:18) S i − √ γ π (cid:15) ijk E j B k (cid:19) (cid:18) S l − √ γ π (cid:15) ljk E j B k (cid:19) ( ρ ∗ h ) − +1 = w . (2.39)Thus, for the resistive MHD, this equation together withEq. (2.34) constitute the simultaneous equations for h and w for given (evolved) values of γ ij , ρ ∗ , S i , S , B i , and E i , and are used for the primitive recovery procedure.Practically, we solve the MHD equations (both in idealand resistive MHD) using the cylindrical coordinates( (cid:36), ϕ, z ) in this work. For this procedure, we write theequations for S , S y , ρ ∗ , and the conservation equationsfor the lepton fraction in a conservative form as in ourviscous hydrodynamics simulations [16, 17]. This is inparticular important to guarantee the conservation of therest mass and angular momentum in numerical compu-tation. E. Numerical methods for solving Maxwell’sequations
We solve the ideal MHD equations mostly in the samemethod as described in Ref. [37]: We evolve the hydrody-namics equations for S , S i , ρ ∗ , and the lepton fractionusing a high-resolution shock capturing scheme (a 3rd-order upwind scheme). In the new implementation, oneimprovement is made for solving the induction equationon the evolution of the magnetic field, B k . The previ-ous implementation solved the equations for B x and B z using a constraint transport scheme [38] and that for B y by an upwind scheme, which is the same as for evolv-ing S , S i , and ρ ∗ . However, we were aware of the factthat this method was so diffusive that the magnetic-fieldenergy was spuriously decreased with a short timescalewithin 100 ms with the typical grid resolution of the gridspacing ∼
200 m (we use the grid resolution of DD2-135M of Ref. [17] for most of the present simulationsand that of DD2-135L of Ref. [17] for 4 test simulations).Thus, for a long-term simulation with the duration morethan seconds, this scheme is practically useless (althougha sophisticated high-order scheme may solve this prob-lem). In the new implementation, we solve the induc-tion equation for B x and B z by a high-resolution upwindscheme, which is the same as for evolving S , S i , and ρ ∗ , and the induction equation for B y by the 4th-ordercentered finite differencing with the operation of the 5th-order Kreiss-Oliger dissipation. To approximately pre-serve the divergence-free condition of Eq. (2.17), a di-vergence cleaning is introduced (in the same manner asin the resistive MHD scheme: see below). We find thatin this method, the stable numerical evolution is feasi-ble and also spurious oscillation is suppressed with thereasonable choice for the coefficient of the Kreiss-Oligerterm. By this procedure, the diffusive evolution of themagnetic field is suppressed in a much better mannerthan in our previous implementation. The maximum er-ror size of the divergence-free condition, which is definedby | ( ∂ i B i ) ∆x/ B z | where ∆x is the grid spacing that cov-ers remnant neutron stars, remains to be ∼ − for theentire simulation time in the present numerical resolu-tion.In the resistive MHD, we also employ the same shockcapturing scheme as that employed in the ideal MHDfor the hydrodynamics part that evolve S , S i , ρ ∗ , andthe lepton fraction. On the other hand, Maxwell’s equa-tions are solved in a different procedure. In addition,we change the method of the time evolution depend-ing on the magnitude of the conductivity. In the fol-lowing, we describe the methods for the low- and high-conductivity cases separately, focusing only on the casethat the timescale of 1 / (4 πσ c ) is shorter than the timestep interval of the numerical simulation, ∆t : A partiallyimplicit scheme is employed for appropriately handlingthe equation of the electric field for the high values of σ c .For the low-conductivity case, we evolve the electro-magnetic equations using the method that is used for evolving geometrical variables and incorporating a di-vergence cleaning prescription. Specifically, we rewriteEqs. (2.14) and (2.15) in the forms: ∂ t E i = − ∂ k (cid:0) β i E k − β k E i + α(cid:15) kij B j (cid:1) − π (cid:0) J i − Qβ i (cid:1) , (2.40) ∂ t B i = − ∂ k (cid:0) β i B k − β k B i − α(cid:15) kij E j (cid:1) , − αγ / γ ij ∂ j φ B (2.41) ∂ t φ B = β k ∂ k φ B − ακφ B − αγ − / ∂ k B k , (2.42)where φ B is a new auxiliary variable associated with thedivergence cleaning, κ is a constant, and J i − Qβ i iswritten as J i − Qβ i = Qv i + σ c α d B k u k ( v i + β i )+ ασ c (cid:104) wA ij E j + (cid:15) ijk u j B k + α d (cid:0) − w B i + (cid:15) ijk u j E k (cid:1) (cid:105) , (2.43)with A ij = δ ij − w − ¯ u i u j and ¯ u i = γ ij u j . We noticethat in our scheme, we always evaluate Q by ∂ k E k / π .In numerical implementation, the transport terms inEqs. (2.40) and (2.41) are evaluated by the 4th-ordercentered finite differencing, and the transport term inEq. (2.42) is evaluated by the 4th-order upwind schemeas done for the geometric variables. Whenever we evolve E i , B i , and φ B , we operate the Kreiss-Oliger dissipationprocedure in the same manner as we do for the geometricvariables.To evolve the system, we basically use the 4th-orderRunge-Kutta scheme but for E i we partially employan implicit scheme. This is implemented in the fol-lowing manner: In Eq. (2.40), we move the term of4 πσ c αwA ij E j that appears in the current term to theleft-hand side and write the left-hand side of the evolu-tion equation in the form ∂ t E i + 4 πσ c αwA ij E j = 1 ∆t (cid:0) E i + 4 πσ c αwA ij E j ∆t − E i (cid:1) := 1 ∆t (cid:0) M ij E j − E i (cid:1) , (2.44)where E i denotes the values of E i in a previous time-stepand M ij = (1 + ζ ) δ ij − ζ ¯ u i u j /w with ζ = 4 πσ c αw∆t .Here, the inverse of M ij is simply written as( M − ) ij = 11 + ζ δ ij + ζ (1 + ζ )( w + ζ ) ¯ u i u j . (2.45)Thus, the updated values of E i is easily obtained in eachRunge-Kutta time step as E i = ( M − ) ij (cid:20) E j − ∆t F jT − π∆t (cid:110) Qv j − σ c ( − α d B k u k )( v j + β j )+ ασ c (cid:16) (cid:15) jkl u k B l + α d ( − w B j + (cid:15) jkl u k E l ) (cid:17)(cid:111)(cid:21) , (2.46)where F jT denotes the term associated with the transportterm, which is evaluated in the explicit manner. We sup-pose that α d is small (or zero) and thus the last term ofEq. (2.46) is not included for the implicit prescription,although it is straightforward to include it.As we illustrate in Appendix A, this implementationenables to solve several test-bed problems successfullyeven for the case that σ c is large such that ζ (cid:29)
1. Theelectromagnetic fields are also evolved stably even in theabsence of the upwind scheme.For the high-conductivity case close to the ideal MHDcase (specifically for ζ ≥
100 in our present implemen-tation), we modify this procedure. The reason for themodification is that the equation of the magnetic fieldapproximately reduces to the induction equation (with aweak diffusion term: see below) as ∂ t B i ≈ ∂ k ( v i B k − v k B i ) + O ( σ − ) . (2.47)Here in Eq. (2.47) we omit the divergence cleaning termfor simplicity, although we include it in the practicalsimulations. As a natural consequence of reducing ap-proximately to a transport equation, for the case of thelong resistive dissipation timescale, the centered finitedifferencing could introduce a long-term growing unsta-ble mode in numerical computation, in particular for theevolution of the poloidal magnetic field leading to theviolation of divergence-free condition of ∂ k B k = 0 (notethat this condition is written only by the poloidal-fieldcomponent in the axially symmetric case). Below, wefirst derive Eq. (2.47) from Eq. (2.41).For the case that ζ (= 4 πσ c αw∆t ) (cid:29)
1, the expres-sion of the solution by the implicit method of Eq. (2.46)should be mathematically modified by replacing ( M − ) ij by ( M − ) ij = 1 ζ δ ij + 1 w + ζ ¯ u i u j . (2.48)This is easily understood that the stiff term associatedwith σ c in the right-hand side of the evolution equationfor E i enforces the exponential decay of the initial electricfield. Then, Eq. (2.46) is rewritten into the form E i = ( M − ) ij (cid:16) E j − ∆t F j − π∆t ( Qv j + α d σ c d j ) (cid:17) − w (cid:15) jkl u k B l =: ˜ E i − w (cid:15) jkl u k B l , (2.49)where d j = B k u k ( v j + β j ) + α ( − w B j + (cid:15) jkl u k E l ) denotesthe term associated with the dynamo. Hence, the termsassociated with the transport for E i (see Eq. (2.41)) isrewritten as − β i B k + β k B i + α(cid:15) kij E j = v i B k − v k B i + α(cid:15) kij ˜ E j . (2.50)Thus, Eq. (2.47) is derived from Eq. (2.41) (here we sup-pose again that α d is small). As in the ideal MHD case, the transport part ofEq. (2.47) should be evaluated by using an upwindscheme for the numerical stability. On the other hand,the last term of Eq. (2.50) constitutes a diffusion term inEq. (2.41), which can be evaluated by the centered finitedifferencing. In our numerical implementation, the trans-port term is handled in the same manner as in the idealMHD case, while the diffusion term is handled by the4th-order centered finite differencing. In addition, thedivergence cleaning and 5th-order Kreiss-Oliger dissipa-tion are implemented as in the relatively low conductivitycase.Before closing this subsection, we briefly comment onour artificial prescription for handling the region in whichthe electromagnetic energy density is much larger thanthe rest-mass energy density, because such regions oftencannot be solved accurately in the MHD simulation. Inthe present paper, if the ratio of B / (4 πρc ) (in the idealMHD case) or ( E + B ) / (8 πρc ) (in the resistive MHDcase) is larger than 10 , we set u i = 0 artificially. In ad-dition, if w (= αu t ) exceeds 5 or ε exceeds 10 after theprimitive-recover process, we also artificially set u i = 0and ε is determined from the condition that the tempera-ture is kT = 0 . k is the Boltzmann constant.The reason for this is that in these high-electromagnetic-energy cases, the primitive recovery often fails. The re-gion with these extreme conditions appears for the casethat the magnetic winding proceeds significantly and asa result a strong outflow is driven from the neutron-starsurface toward a low-density atmosphere. In the absenceof these artificial prescriptions, the computation oftencrashed in our current implementation. Since we give upresolving the region with w ≥
5, we cannot explore avery-high Lorentz-factor jet in the present simulations.
F. Initial condition and relevant timescales
As in our series of the papers [11, 17, 25], the ini-tial condition for the matter field is supplied from theresult of a simulation for binary neutron star mergers.Specifically, we employ the DD2-135M model of Ref. [17]:a merger remnant of binary neutron stars with eachneutron-star mass 1 . M (cid:12) . As touched in Sec. II E, foursimulations are also performed using a lower-resolutionmodel, DD2-135L, to check the dependence of the resultson the grid resolution.We then superimpose electromagnetic fields for whichthe energy density is much smaller than the internal andkinetic energy density of the fluid. In this work, we ini-tially give a poloidal magnetic field written in the form B (cid:36) = − (cid:36) ∂ z A ϕ , (2.51) B z = 1 (cid:36) ∂ (cid:36) A ϕ , (2.52)where A ϕ = A (cid:36) max (cid:18) PP max − δ A , (cid:19) , (2.53)with δ A = 10 − and P max is the maximum pressure. A determines the initial field strength. With this setting,the magnetic field is present for a high-density regionwith ρ (cid:38) . g / cm . The toroidal magnetic field is setto be zero ( B y = 0) initially. The electric field is de-termined by the ideal MHD condition (2.21). We alsoperformed two ideal MHD simulations with δ A = 10 − (for which the settings agree approximately with Bh-R0-Y and Bh-R0-N in Table I). Irrespective of the values of δ A , the magnetic-field is strong only inside the remnantmassive neutron star for which the magnetic-field effectappears in the most remarkable manner in the presentcontext (i.e., in the axisymmetric simulation with no dy-namo effect). Indeed, it is found that the results onthe evolution of the electromagnetic energy and angu-lar velocity profile depend only weakly on the value of δ A . Thus in the following, we focus only on the resultsof δ A = 10 − .As Eq. (2.53) shows, we pay attention to the MHD ef-fect only for the high-density region, i.e., for the remnantneutron star, and do not pay attention to the MHD ef-fect in the disk. Because the angular velocity of the diskdecreases along the cylindrical-radius direction, the MRIwould play a crucial role for its evolution. For explor-ing the evolution, we need a long-term non-axisymmetricsimulation that can capture the effects of turbulent mo-tion and resulting effective viscosity induced. This is be-yond the scope of this paper.We also performed simulations initially with a purelytoroidal magnetic field of B y = A (cid:36)z max (cid:18) PP max − − , (cid:19) . (2.54)Here the dependence on the coordinates, (cid:36)z , stems fromthe required symmetry for B y in the axially and planely(with respect to the z = 0 plane) symmetric assumption.In axisymmetric simulations with this initial condition,the magnetic field should be simply preserved or decaywith the timescale determined by σ c (cf. Eq. (2.59) be-low). In Appendix B, we show that our implementationcan derive the expected numerical results.Table I lists the initial conditions (with the purelypoloidal magnetic field). ˆ B init and E B denote the maxi-mum value of B z γ − / (= B z γ − / ) and the electromag-netic energy at the initial state. Here, the electromag-netic energy is defined by E B = 18 π (cid:90) ( B + E ) √− g d x. (2.55)Note that E is much smaller than B . Since the volumeof the neutron star is of the order of 10 cm , the volumeaveraged magnetic-field strength is initially ∼ G for E B ∼ erg. The labels, Bv, Bh, Bm, and Bl, in the model name of Table I refer to the initial magnetic-fieldstrength (very high, high, medium, and low).Because the remnant neutron star which we employis differentially rotating, with this setting, the toroidalmagnetic-field strength is initially increased linearly withtime by the winding effect as (e.g., Refs. [21–23]) B T ≈ B (cid:36) d Ω d ln (cid:36) t, (2.56)where B (cid:36) is the initial local magnitude of the (cid:36) compo-nent of the magnetic field. This growth continues untilthe magnetic-field energy increases to ∼
10% of the ro-tational kinetic energy, E kin , of the neutron star. Thusthe winding timescale is defined approximately by τ wind := (cid:115) E kin E B, (Ω max − Ω ) − ≈ . (cid:18) E kin erg (cid:19) / (cid:18) E B, erg (cid:19) − / × (cid:18) Ω max − Ω / s (cid:19) − , (2.57)where Ω max and Ω denote the maximum value of theangular velocity and angular velocity at the center, re-spectively. E B, is the initial magnetic-field energy, and E kin is defined by E kin := 12 (cid:90) ρ ∗ hu ϕ v ϕ d x. (2.58)For the model employed in this paper, the initial value of E kin is ≈ . × erg. In the estimate of Eq. (2.57),we supposed that the region of the maximum angular ve-locity would be located in an outer region of the remnantmassive neutron star (see, e.g., Ref. [16]) and it governsthe winding and the growth of the electromagnetic en-ergy. Also, we assumed that a part of the poloidal mag-netic field for which the electromagnetic energy is ∼ E B, contributes to the winding.In the case of the finite value of σ c , the electromag-netic fields are dissipated. The dissipation timescale isapproximately written as τ dis ≈ π(cid:96) σ c c = 0 .
13 s (cid:18) (cid:96) (cid:19) (cid:16) σ c s − (cid:17) , (2.59)where (cid:96) denotes the typical variation scale of the mag-netic field. Because we do not know the typical size of σ c ,we pay attention to the cases for which τ dis is comparableto τ wind , i.e., σ c = 10 –10 s − .For the case of τ dis (cid:29) τ wind , obviously, the resis-tive dissipation does not play a role. By contrast, for τ dis (cid:28) τ wind , the magnetic field is diffused out before themagnetic winding significantly amplifies the magnetic-field strength. Also, the winding effect does not directlydetermine the final magnetic-field profile for τ dis (cid:28) τ wind . TABLE I. Initial condition and set-up for the numerical sim-ulation. We list the maximum value of ˆ B init = B z γ − / , themagnetic-field energy, and the value of σ − in units of mil-liseconds. The last column shows the on or off of the irradi-ation/heating and pair-annihilation of neutrinos. For all theinitial conditions, the total baryon mass is M ∗ = 2 . M (cid:12) ,the gravitational mass is M = 2 . M (cid:12) , the total rotationalkinetic energy is E kin ≈ . × erg, and the total angularmomentum is J = 4 . × g cm /s.Model ˆ B init (G) E B (erg) σ c (s − ) Pair & Heat.Bv-R0-Y 4 . × . × ∞ OnBh-R0-Y 2 . × . × ∞ OnBh-R0-N 2 . × . × ∞ OffBh-Rw-Y 2 . × . × OnBh-Rl-Y 2 . × . × OnBh-Rm-Y 2 . × . × OnBh-Rh-Y 2 . × . × OnBm-R0-Y 1 . × . × ∞ OnBm-R0-N 1 . × . × ∞ OffBm-Rl-Y 1 . × . × OnBm-Rm-Y 1 . × . × OnBm-Rh-Y 1 . × . × OnBl-Rl-Y 5 . × . × On This point is understood by the following model for theevolution of B y , for which the evolution equation is ap-proximately written as ∂ t B y ≈ (cid:36) B i ∂ i Ω + η ∆ B y , (2.60)where ∆ denotes the Laplacian operator (different from ∆ ). Here, we picked up only the terms related to thewinding and resistive dissipation. For this equation, weconsider a toy model in which the first term in the right-hand side is not time-varying as (cid:36) B i ∂ i Ω = F ( x i ) =: η ∆ F ( x i ). Then, it is easily found that for t (cid:28) τ dis , themagnetic winding approximately determines the solutionas B y ≈ tF ( x i ). On the other hand, for t (cid:38) τ dis , theasymptotic solution is B y ≈ − F ( x i ). Thus, the windinghistory is not reflected, although the order of magnitudeof F is written as ∼ Ω τ dis |B i | , which reflects the windingin t ≤ τ dis . We note that near the rotational axis, F ∝ (cid:36) , and hence, F ∝ (cid:36) . Thus, even if B y is positiveduring the winding, the final relaxed value for it can benegative. III. EVOLUTION OF A REMNANT NEUTRONSTARA. Summary of the evolution
Table I lists the models for which we performed simu-lations. The initial electromagnetic-field energy is variedamong 1 . × –9 . × erg. With this setting, themagnetic winding timescale is approximately 0.1–1 s. The simulations are performed for the ideal MHD mod-els (referred to as R0 in the model name) and for the re-sistive MHD models with σ c = 10 , 10 , 10 , and 10 s − (referred to as Rh, Rm, Rl, and Rw in the model name).For σ c = 10 s − , the dissipation timescale by the resis-tivity is ∼
100 s, which is much longer than the windingtimescale in our setting. Thus, this model is employedto check that the ideal MHD results are approximatelyreproduced in the resistive MHD implementation.For most of the models, neutrino effects including theirradiation/heating and pair annihilation are taken intoaccount as in our previous papers [17, 25]. To clarifythe importance of the neutrino effects on the evolutionof the massive neutron stars and on the mass ejection,for the ideal MHD case, we perform two simulations inwhich the effects of both the neutrino irradiation/heatingand neutrino pair annihilation are switched off (modelsBh-R0-N and Bm-R0-N: in the presence of the neutrinoeffects, the model is referred to with the label Y).As we touched in Secs. II E and II F, numerical sim-ulations were performed with two grid resolutions fortwo ideal MHD (Bh-R0-Y and Bh-R0-N) and two resis-tive MHD (Bh-Rw-Y and Bh-Rl-Y) cases. By comparingthe results of two different resolutions, we find that thenumerical results depend weakly on the grid resolution(see Appendix C). Specifically, for the lower grid resolu-tion, the amplification of the magnetic-field strength issuppressed, and thus, the maximum magnetic-field en-ergy becomes smaller, in particular in the presence ofthe neutrino irradiation/heating and pair-annihilation ef-fects. For the resistive MHD simulation, we also find thatthe dissipation is spuriously enhanced for the lower gridresolution. However, besides these quantitative effects,the evolution process of the neutron star is not modifiedqualitatively by changing the grid resolution.In the present context, the typical MHD process in theremnant neutron star of neutron-star mergers in the pres-ence of an initial seed poloidal magnetic field is governedby the magnetic winding effect, for which the growthtimescale of the toroidal magnetic-field energy is writ-ten approximately by Eq. (2.57). When the electro-magnetic energy exceeds a few % of the rotational ki-netic energy of the remnant neutron star (for this case E kin ≈ . × erg), the magnetic braking starts play-ing an important role. As a consequence, the angularmomentum in the neutron star is redistributed and thedifferentially rotating configuration of the angular veloc-ity is enforced to approach a (approximately) rigidly ro-tating one. Then, the angular velocity in the centralregion exceeds ∼ ∼ . -3 -2 -1 t E B ( e r g ) t (s) Bv-R0-YBh-R0-YBm-R0-YBh-R0-NBm-R0-NBh-Rw-Y -3 -2 -1 t E B ( e r g ) t (s) Bh-Rh-YBh-Rm-YBh-Rl-YBh-Rw-YBm-Rh-YBm-Rm-YBm-Rl-YBl-Rl-Y
FIG. 1. Evolution of the electromagnetic energy for all the ideal MHD models, Bv-R0-Y, Bh-R0-Y, Bh-R0-N, Bm-R0-Y,Bm-R0-N, and for a resistive model with a high conductivity, Bh-Rw-Y (left panel) and for all the resistive models employedin this paper (right panel). addition, the magnetic-field profile in the neutron star ismodified significantly by the mass outflow (and result-ing magnetic-flux escape) process. The detail on thesefindings is described in Sec. III A 1.In the presence of only low resistivity (i.e., a high valueof σ c ), the evolution process of the remnant neutron staris essentially the same as in the ideal MHD case. In thepresence of a high resistivity with which the dissipationtimescale of the magnetic field is shorter than the windingtimescale, not only the winding but also the magnetic-field dissipation plays an important role. For such cases,the magnetic winding does not occur significantly andthe magnetic-field effect on the angular velocity profile isminor (see Sec. III A 2 for details).
1. Ideal MHD
Figure 1 shows the electromagnetic energy as a func-tion of time. The left panel plots the results for all theideal MHD models together with one resistive model ofa high conductivity (Bh-Rw-Y) and the right panel plotsthose for all the resistive MHD models. In the idealMHD cases, the electromagnetic energy, E B , is initiallyincreased by the magnetic winding effect by several or-ders of magnitude. In this early stage, E B is approxi-mately proportional to t irrespective of the models (in-cluding the resistive MHD models except for the case of σ c = 10 s − for which the resistive dissipation timescaleis quite short; the order of 10 ms).For the ideal MHD case, after E B exceeds ∼
3% of E kin , the increases of E B is decelerated (i.e., the slope of E B becomes shallower than t ), indicating that the wind-ing effect is weaken and the magnetic braking becomesimportant. However, the increase of E B still continuesuntil it reaches ∼
10% of E kin , and then, the growthof the electromagnetic field is saturated. The satura- tion energy is slightly smaller for the smaller initial fieldstrength. This is primarily due to the fact that moremagnetic fluxes are escaped associated with the mass out-flow for longer-term evolution, and may be partly due tothe numerical dissipation, i.e., in the longer-term evolu-tion the damage by the numerical dissipation is more se-rious: We note that the simulations are always performedfor ∼ rotational period of the neutron star.The saturated energy is slightly higher in the absenceof the neutrino irradiation/heating and pair-annihilationeffects. Our interpretation for this is that in the ab-sence of these effects, the mass outflow from the neutronstar is suppressed (cf. Sec. III B), resulting escape of themagnetic flux (associated with the mass outflow) is alsosuppressed, and as a result, the amplification of the elec-tromagnetic field proceeds to a higher level.As we already mentioned briefly, the mass ejection isdriven primarily by the neutrino irradiation/heating (seeSec. III B for details). However, with the increase of theelectromagnetic energy toward the saturation, the massejection from the neutron star is slightly enhanced. Thisis due to the magnetic pressure primarily resulting fromthe significantly amplified toroidal magnetic field. Bythis enhancement of the outflow, a part of the magneticflux escapes from the neutron star, and as a result, thetotal electromagnetic energy decreases. However, the de-crease rate becomes eventually low ( E B eventually showsthe oscillatory behavior) and the electromagnetic energyin average appears to relax to ∼ erg, i.e., ∼
1% ofthe rotational kinetic energy of the neutron star. Wenote that the rotational kinetic energy does not changesignificantly by the magnetic braking effects because it isalways larger by more than one order of magnitude thanthe electromagnetic energy, and thus, the magnetic-fieldeffect is minor. Instead, the rotational kinetic energyslightly increases with time because the neutron star con-tracts due to the emission of neutrinos of total energy of0 x (km)0200040006000800010000 Ω (r a d / s ) Bh-R0-Y( σ c = ∞ ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s x (km)0200040006000800010000 Ω (r a d / s ) Bm-R0-Y( σ c = ∞ ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s x (km)0200040006000800010000 Ω (r a d / s ) Bm-Rl-Y( σ c = 10 s − ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s x (km)0200040006000800010000 Ω (r a d / s ) Bm-Rm-Y( σ c = 10 s − ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s x (km)0200040006000800010000 Ω (r a d / s ) Bm-Rh-Y( σ c = 10 s − ) t = 0 .
00 s t = 0 .
10 s t = 0 .
20 s t = 0 .
30 s t = 0 .
40 s t = 0 .
68 s x (km)0200040006000800010000 Ω (r a d / s ) Bh-Rw-Y( σ c = 10 s − ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s
FIG. 2. Evolution of the angular velocity profile for models Bh-R0-Y (top left), Bm-R0-Y (top right), Bm-Rl-Y (middle left),Bm-Rm-Y (middle right), Bm-Rh-Y (bottom left), and Bh-Rw-Y (bottom right). We note that due to the neutrino cooling,the remnant neutron star contracts and its central region spins up. For model Bm-Rh-Y, the increase of the angular velocityis predominantly due to this effect. ∼ erg, and hence, the neutron star slightly spins upin the absence of the electromagnetic effect, as shown inRef. [17].For the nearly ideal MHD case, i.e., for the resis-tive MHD case with a tiny resistivity (high conductivity, σ c = 10 s − ; model Bh-Rw-Y), the curve of E B is sim-ilar to the corresponding ideal MHD model (Bh-R0-Y)as found in the left panel of Fig. 1. This is natural be-cause for this model, the resistive dissipation timescaleis quite long ∼
100 s (cf. Eq. (2.59)). In the late stagewith t (cid:38) x (km)10 B P ( G ) Bm-R0-Y( σ c = ∞ ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s B T ( G ) Bm-R0-Y( σ c = ∞ )10 z (km) − − − t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s x (km)10 B P ( G ) Bm-R0-N( σ c = ∞ ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s B T ( G ) Bm-R0-N( σ c = ∞ )10 x (km) − − − t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s x (km)10 B P ( G ) Bm-Rm-Y( σ c = 10 s − ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s B T ( G ) Bm-Rm-Y( σ c = 10 s − )10 x (km) − − − t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s z (km)10 B P ( G ) Bm-R0-Y( σ c = ∞ ) t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s B T ( G ) Bm-R0-Y( σ c = ∞ )10 z (km) − − − t = 0 .
00 s t = 0 .
10 s t = 0 .
30 s t = 0 .
80 s t = 1 .
20 s t = 1 .
40 s t = 1 .
60 s t = 1 .
80 s
FIG. 3. Poloidal (left) and toroidal (right) magnetic-field strength along the x -direction with z = 1 km at selected time slices formodels Bm-R0-Y (top panel), Bm-R0-N (2nd top panel), and Bm-Rm-Y (3rd panel). The bottom panel shows the poloidal andtoroidal magnetic-field strengths along the z -direction with x = 1 km for model Rm-B0-Y. The field strengths of the poloidaland toroidal components are defined by B P = √B x B x + B z B z γ − / and B T = sign( B y ) (cid:112) B y B y γ − / ( ≈ B y γ − / ), respectively. t (cid:38) As we mentioned in Sec. I,initially, the angular velocity is an increase function of (cid:36) in the central region of the neutron star. This profileis modified and approaches gradually a rigidly rotatingstate for the region in which the strong magnetic fieldis present. Figure 2 illustrates that the angular velocityin the central region is increased by a factor of ∼ τ a = RB T / √ πρh ≈ .
11 s (cid:18) R
10 km (cid:19) (cid:18) B T G (cid:19) − (cid:18) ρ g / cm (cid:19) / , (3.1)where R , B T , and ρ denote the typical radius, toroidalmagnetic-field strength, and rest-mass density of the neu-tron star, and we set h ≈ τ a ,because the angular velocity approximately obeys a hy-perbolic equation [21]. However, in the present context,the poloidal magnetic flux is significantly escaped by theneutrino-driven mass outflow and a simple oscillation isnot seen. We also note that in the non-axisymmetric case,the turbulence and dynamo effects also would modify themagnetic-field profile significantly [23].The top and bottom panels of Fig. 3 show the poloidal(left) and toroidal (right) magnetic-field strength along x -direction with z = 1 km (top) and along the z -direction with x = 1 km (bottom) for model Rm-B0-Y. The poloidal and toroidal magnetic-field strengthsare defined, respectively, by √B x B x + B z B z γ − / andsign( B y ) (cid:112) B y B y γ − / ( ≈ B y γ − / ). It is found that thetoroidal-field strength (absolute value) is increased farbeyond 10 G inside the neutron star quite uniformly.By contrast, the maximum value of the poloidal fieldstrength does not increase but rather decreases in thelate phase. This decrease is due to the fact that associ-ated with the matter outflow, the magnetic flux is ejected We note that due to the neutrino cooling, the remnant neutronstar contracts and its central region spins up. The increase ofthe angular velocity is partly due to this effect. For the high-resistivity cases such as models Bh-Rh-Y and Bm-Rh-Y, theincrease of the angular velocity in the central region is causedmainly by this effect: This was confirmed by comparing the nu-merical result with no magnetic fields [17]. from the neutron star. By this process, instead, a strongpoloidal field is produced in the polar direction outsidethe neutron star (see the bottom panels of Fig. 3): Thetypical poloidal-field strength at the neutron-star pole is ∼ G and ∼ G at z ∼
100 km.Associated with the mass outflow from the neutron-star surface, magnetic-fluxes escape toward a large (cid:36) and large z direction. Because the angular velocity ofthe escaped matter decreases with the motion toward (cid:36) -direction, the generation of a negative toroidal magneticfield is proceeded along the field lines of the escaped mat-ter. This effect subsequently makes the toroidal magneticfield in the neutron star negative. For model Bm-R0-Y,the mass outflow from the neutron-star surface is appre-ciably caused by the neutrino-driven wind. As a result,the region with the negative toroidal field expands grad-ually with time (see the top-right panel of Fig. 3). Bycontrast, for model Bm-R0-N, the neutrino-driven windis absent and the mass outflow from the neutron-star sur-face is weak. Consequently, the region with the negativetoroidal field appears only in the vicinity of the neutronstar surface (see the second-top-right panel of Fig. 3).As indicated in Fig. 3, the toroidal magnetic-field en-ergy dominates over the poloidal magnetic-field one af-ter the significant winding: The poloidal-field energy is < .
1% of the toroidal field energy. However, this maybe an artifact in axisymmetric computations in which nodynamo effect is present [26]. In reality, the poloidal-field strength may become much higher than the presentresults.
2. Resistive MHD
For the resistive MHD models with low values of theconductivity (i.e., high resistivity), the magnetic-fieldstrength and the maximum value of E B reached by themagnetic winding depend on the values of σ c , and thus,are different from those for the ideal MHD case (see theright panel of Fig. 1 as well as Figs. 2 and 3). Specifi-cally, the maximum magnetic-field strength depends onthe ratio, τ dis /τ wind =: R τ : For smaller values of R τ (i.e., R τ (cid:46) E B are smaller because the resistive dissipa-tion of the electromagnetic energy proceeds faster thanthe magnetic winding. Thus, for the smaller values ofthe initial electromagnetic energy, these maximum val-ues are smaller for a given value of σ c . This is foundclearly from the results for σ c = 10 and 10 s − (com-pare the results of models Bh-Rm-Y and Bh-Rh-Y withBm-Rm-Y and Bm-Rh-Y, respectively). By contrast, for R τ (cid:38)
1, the resistive effect is not very remarkable: Forexample, the maximum value of E B always exceeds 1% ofthe rotational kinetic energy, E kin , because the magneticwinding sufficiently proceeds until the saturation of thegrowth is reached: See the results with σ c ≥ s − .We note that irrespective of the maximum value of E B reached, the toroidal magnetic-field energy dominates3 FIG. 4. Profiles for the density, temperature, specific entropy, electron fraction (left for these four plots), poloidal and toroidalmagnetic-field strength, the ratio of the rest-mass to magnetic energy density, and the ratio of the gas pressure to the magneticpressure (right for these four plots) after the saturation of the growth of the electromagnetic field for models Bm-R0-Y (toppanels) and Bm-Rm-Y (bottom panels). over the poloidal magnetic-field one in the end as in theideal MHD case. Again, this could be an artifact in ax-isymmetric computations with no dynamo effect in whichthe poloidal field cannot be amplified [26]. In the non-axisymmetric case, this is unlikely to be the case, andhence, the poloidal field may be also enhanced in thelater stage.For the case that the maximum value of E B is smallerthan ∼
3% of E kin , the effect of the magnetic brakingis weak, and hence, the differential rotation is modifiedweakly. As found from the middle-right and bottom-left panels of Fig. 2, the angular velocity near the ro-tational axis ( x = 0) does not increase appreciably for σ c ≤ s − (see also footnote 1). By contrast, theangular-velocity profile for σ c ≥ s − evolves in asimilar manner to that in the ideal MHD case (comparethe top-left and bottom-right panels or the top-right andmiddle-left panels of Fig. 2).The third row of Fig. 3 shows the evolution of themagnetic-field profile near the equatorial plane for aresistive model (Bm-Rm-Y). Figure 3 illustrates that the magnetic-field strength depends strongly on σ c : For σ c ≤ s − , the magnetic field is dissipated in 0 . (cid:46) . R τ is small, <
1, andhence, the electromagnetic energy caused by the windingis much smaller than the rotational kinetic energy andinternal energy of the neutron star. Hence, the gener-ated thermal energy plays only a minor role for the massoutflow and evolution of the neutron star (although thismay contribute a bit to enhancing the mass ejection). For R τ (cid:38)
1, the electromagnetic energy reached is ∼
10% of4the rotational kinetic energy. For such cases, however,the dissipation timescale is quite long; in our presentsetting, it is longer than ∼
3. Outcome after the winding
Figure 4 displays the profiles for the density, tempera-ture, specific entropy, electron fraction (left for these fourplots), poloidal and toroidal magnetic-field strength, theratio of the rest-mass to magnetic energy density, andthe ratio of the gas pressure to the magnetic pressure(right for these four plots) after the saturation of theelectromagnetic energy is reached for models Bm-R0-Y(ideal MHD model: top panels) and Bm-Rm-Y (moder-ately high resistivity model: bottom panels). These plotsshow that the system is composed of a massive neutronstar surrounded by a disk and an outflow driven from thepolar region of the neutron star toward the z -direction ir-respective of the presence or absence of the resistivity.However, by comparing the profiles for the two mod-els, two quantitative differences are found. First, themagnetic-field strength for Bm-R0-Y is stronger than forBm-Rm-Y as we already mentioned in Sec. III A 2 (seeFig. 3). This is clearly reflected outside the neutron starin Fig. 4. Second, the property of the outflow drivenfrom the polar surface of the neutron star depends on thestrength of an MHD effect: In the presence of the strongMHD pressure (i.e., the ideal MHD case), the matter den-sity in the polar region is smaller, and also, the velocity ofthe outflow is higher. This modifies the electron fractionof the matter in the polar region and ejecta: The elec-tron fraction is lower for model Bm-R0-Y for which theoutflow velocity is higher and the effect of the neutrinoirradiation to the ejecta is weaker (see also Sec. III B).
4. Remark: Comparison with previous work
In the present work, we find that the massive neu-tron star does not contract significantly and its structureis not modified essentially, after the magnetic windingand braking processes. This conclusion is in contrast toprevious work [39–41], which showed that differentiallyrotating neutron stars significantly contract, and if theyare hypermassive, they collapse to a black hole, after themagnetic winding effect and associated angular momen-tum transport by the magnetic braking take place. Themajor reason for the absence of the significant contrac-tion is that in the present work, a realistic angular veloc-ity profile is determined by the merger simulation (i.e., d Ω /d(cid:36) > . . . . . . . . . t (sec)012345 M e j e ( − M (cid:12) ) M ej , tot , B0-YBh-R0-YBm-R0-YBh-R0-NBm-R0-N FIG. 5. Ejecta mass as a function of time for the ideal MHDcases. Model B0-Y refers to the result in the absence of themagnetic-field effect but with the neutrino effects. locity decreases steeply with the cylindrical radius. For d Ω /d(cid:36) >
0, the centrifugal force in the central regionof the neutron star increases after the magnetic windingand subsequent magnetic braking proceed, whereas for d Ω /d(cid:36) < d Ω /d(cid:36) <
0, while for d Ω /d(cid:36) > B. Ejecta
As we reported in our previous papers [11, 17, 25],the mass ejection proceeds from the remnant of bi-nary neutron star mergers through the neutrino irradi-ation/heating even in the absence of any other effectslike the MHD and viscous effects. The rest mass ofthis neutrino-irradiation component is not very large as (cid:46) − M (cid:12) [17, 25]. In the context of the present paper,the mass ejection may be enhanced by the MHD effectthat results primarily from the amplified toroidal mag-netic field. In this subsection, we pay attention to thisenhancement. Because we finally find that the mass ejec-tion by the magnetic pressure enhanced by the magneticwinding in the neutron star is a minor effect, we here payattention to this topic focusing only on the ideal MHDresults.The ejecta component is determined using the samecriterion as in Refs. [16, 17]; we identify a matter com-ponent with | hu t | > h min as the ejecta. Here h min de-notes the minimum value of the specific enthalpy in theadopted equation-of-state table, which is ≈ . c . For5 . . . . . . . . . t (sec)0 . . . . . v e j e , p o l e ( c ) B0-YBh-R0-YBm-R0-YBh-R0-NBm-R0-N . . . . . . . . . t (sec)0 . . . . . Y e , p o l e B0-YBh-R0-YBm-R0-YBh-R0-NBm-R0-N
FIG. 6. Average velocity (left panel) and electron fraction (right panel) of the ejecta as functions of time for the ideal MHDsimulations and for one simulation without the MHD effect (B0-Y). the matter escaping from a sphere of r = r ext , we definethe ejection rates of the rest mass and energy (kineticenergy plus internal energy) at a given radius and timeby ˙ M eje := (cid:73) ρ √− gu i dS i , (3.2)˙ E eje := (cid:73) ρ ˆ e √− gu i dS i , (3.3)where ˆ e := hαu t − P/ ( ραu t ). The surface integral isperformed at r = r ext with dS i = δ ir r sin θdθdϕ forthe ejecta component. r ext is chosen to be 1000 km inthis work. ρ √− gu t (= ρ ∗ ) obeys the continuity equation for therest mass (see Eq. (2.1)), and thus, the time integra-tion for it is a conserved quantity. Also in the absenceof gravity and magnetic-field effects, ρ ˆ e √− gu t obeys theenergy conservation equation, and far from the centralregion, the sum of its time integration and its gravita-tional potential energy are approximately conserved (as-suming that the electromagnetic energy is much smallerthan the kinetic energy). Thus, the total rest mass andenergy (excluding the gravitational potential energy andelectromagnetic energy) of the ejecta (which escape awayfrom a sphere of r = r ext ) are calculated by M eje , esc ( t ) := (cid:90) t ˙ M eje dt, (3.4) E eje , esc ( t ) := (cid:90) t ˙ E eje dt. (3.5)Far from the central object, E eje , esc is approximatedby E eje , esc ≈ M eje , esc c + U + T kin + GM M eje , esc r ext , (3.6)where U and T kin are the values of the internal energyand kinetic energy of the ejecta at r ext → ∞ , respectively.The last term of Eq. (3.6) approximately denotes the contribution of the gravitational potential energy of thematter at r = r ext with M being the total gravitationalmass of the system [16]. Since the ratio of the internal en-ergy to the kinetic energy of the ejecta decreases with itsexpansion, we may approximate U/T kin ≈
0, and hence, E eje , esc by E eje , esc ≈ M eje , esc c + T kin + GM M eje , esc /r ext for the region far from the central object. We then definethe average velocity of the ejecta (for the component thatescapes from a sphere of r = r ext ) by v eje := (cid:115) E eje , esc − M eje , esc c − GM M eje , esc /r ext ) M eje , esc . (3.7)In the setting of the present paper, the mass outflowand resulting ejection occur only from the polar regionof the neutron star (cf. Fig. 4), because we do not takeinto account the realistic evolution of the disk surround-ing the neutron star, and mass outflow does not occurfrom it. Thus, in the present work, we perform the sur-face integral of Eqs. (3.2) and (3.3) only for θ ≤ ◦ where θ denotes the polar angle. Since the initial condi-tion of our simulations is prepared using the result of amerger simulation for binary neutron stars, the dynam-ical ejecta component is present from the beginning inthe computational domain [17]. However, the dynamicalejecta component is located primarily near the equatorialplane. Thus, with the restriction for the surface integralof θ ≤ ◦ , we can focus approximately on the post-merger mass ejection, which comes primarily from thepolar region of the neutron star in the present context.Figure 5 shows the rest mass of the post-merger ejectacomponent as a function of time for four ideal MHDsimulations. For comparison, we also plot the resultin the absence of the magnetic field (but with the neu-trino effects: model B0-Y) [17]. Figure 5 shows that themass ejection is driven primarily by the neutrino irradi-ation/heating and pair-annihilation heating toward thepolar region, because the rest mass of the ejecta in theabsence of these neutrino effects is less than half of thatin their presence and also less than that in the absence6of the magnetic-field effect (model B0-Y). By comparingthe results of models Bh-R0-Y and Rm-R0-Y with that ofmodel B0-Y, it is still found that the magnetic pressureenhanced by the magnetic winding increases the ejectamass. However, as obviously found from Fig. 5, the restmass of this ejecta component is of the order of 10 − M (cid:12) and by two orders of magnitude smaller than that ofthe viscosity-driven ejecta, which comes from disks (tori)around the neutron star [11, 17]. Therefore, the MHD ef-fect that stems from the winding in the neutron star doesnot contribute appreciably to increasing the ejecta mass.This is likely to be due to the facts that (i) the MHDeffect can strip the material only in a thin polar surfacelayer of the neutron star for which the total rest mass istiny and (ii) the gravitational potential near the neutronstar is so deep that the mass ejection cannot efficientlyoccur from its surface.The present work indicates that the mass ejection bythe MHD effect from the neutron star is much less signif-icant than the viscosity-driven mass ejection from disks(tori) surrounding the neutron star. These results sug-gest that the major source of the mass ejection from theremnant of binary neutron star mergers may be the disk,not the remnant neutron star. However, this speculationshould be explored in the future by the simulations inwhich other important effects such as dynamo effects aretaken into account.A noticeable MHD effect is found in the average veloc-ity and electron fraction of the ejecta (see Fig. 6; cf. Fig. 4as well). In the absence of the MHD effect (model B0-Y), the average velocity of the ejecta is at most 0 . c .However, in the presence of the strong MHD effect (mod-els Bh-R0-Y and Bm-R0-Y) the average velocity of theejecta is enhanced to be ∼ . c at t ∼ (cid:38) .
45, due to the irradiation bothby electron neutrinos and anti electron neutrinos, whilein the presence of the MHD effect, the electron fractionbecomes slightly lower, ∼ .
4. For example, for the idealMHD models Rh-B0-Y and Rm-B0-Y, the electron frac-tion in the polar region is as high as (cid:38) . . IV. SUMMARY
We performed ideal and resistive MHD simulations fora remnant neutron star of binary neutron star mergers ingeneral relativity with neutrino effects. As a first step,we paid attention to the effect of the magnetic windingfor the evolution of the remnant neutron star and re-sulting mass ejection. The initial matter profile for thesimulations was obtained from a merger simulation. Aseed poloidal magnetic field, for which the electromag-netic energy is much smaller than the rotational kineticand internal energy of the system, was initially superim-posed inside the remnant neutron star, and we focusedonly on the evolution of it; we did not pay attention tothe MHD evolution of the disk (torus) surrounding theneutron star.Because of the magnetic winding effect, the toroidalmagnetic field is generated and amplified in this set-ting. Since the toroidal magnetic-field strength increaseslinearly with time, the electromagnetic-field energy in-creases as ∝ t in the early growth stage. When theelectromagnetic-field energy exceeds ∼
3% of the rota-tional kinetic energy, the magnetic braking plays an im-portant role for the redistribution of the angular momen-tum inside the neutron star. For the ideal MHD case, thisalways occurs and, as a result of the angular momentumredistribution, the angular velocity, which is initially anincrease function of the cylindrical radius, is enforced toapproach a rigidly rotating state for the region in whichthe magnetic braking works well. The maximum elec-tromagnetic energy reached is ∼
10% of the rotationalkinetic energy of the neutron star, i.e., the maximumvalue of E B is ∼ erg. By the angular-momentum re-distribution, the rotational kinetic energy of the neutronstar is not significantly changed, while at the late stage,the electromagnetic energy relaxes to ∼ erg after themagnetic braking is exerted, because the magnetic fluxescapes from the neutron star associated with the massoutflow.For the resistive MHD simulation, the maximum elec-tromagnetic energy reached by the winding effect de-pends strongly on the magnitude of the conductivity;specifically, the maximum value is determined by the ra-tio, R τ = τ dis /τ wind . For R τ (cid:38)
1, the magnetic windingproceeds until the electromagnetic energy reaches ∼ R τ <
1, the resistive dissipationplays a role for suppressing the growth of the electro-magnetic energy. If the electromagnetic energy does notreach ∼
3% of the rotational kinetic energy, the magneticbraking effect is weak, and hence, the differentially rotat-ing state of the neutron star is preserved; the increase ofthe angular velocity near the central region is suppressed.7It is also found that the magnetic field amplified by thewinding in the neutron star does not contribute much toenhancing the mass ejection; for the mass ejection fromthe neutron star, the neutrino effects such as neutrino ir-radiation/heating and pair-annihilation heating are moreimportant than the MHD effect. Our interpretation forthis result is that (i) the MHD effect in the neutron starcan strip the material only in the thin polar surface re-gion of the neutron star for which the total rest mass istiny and (ii) the gravitational potential near the neutronstar is so deep that the mass outflow cannot be efficient.The present result suggests that the major source of themass ejection from the remnant of binary neutron starmergers may not be remnant neutron star but the disksurrounding it. However, this speculation should be ex-amined in more detail by the simulations in which otherimportant effects such as dynamo effects are taken intoaccount.In the present work, we performed axisymmetric sim-ulations. This implies that we might overlook importantMHD effects. One important missing effect is the turbu-lence induced by Parker [42] and Taylor instabilities [43](e.g., Ref. [44]). As we showed in this paper, the toroidalmagnetic field is inevitably amplified in the presence ofa seed poloidal magnetic field. It is well-known that inthe presence of the strong toroidal magnetic fields, (non-axisymmetric) Parker and Taylor instabilities can inducea turbulence, which could further induce the magnetic-field amplification through the dynamo effect. It is notclear at all what happens in such a situation. In addi-tion, it is not very clear what the final configuration ofthe magnetic field profile is after the onset of the turbu-lence (e.g., Ref. [45]). To explore these questions, weobviously need three-dimensional MHD simulations inthe subsequent work. An alternatively phenomenologi-cal approach is to incorporate a dynamo term to inducethe turbulence state [35, 36]. By this prescription, wedo not need three-dimensional simulations: An axisym-metric simulation would be reasonable to obtain someinsight on the effect of the turbulence. We plan to per-form simulations with a dynamo term in the subsequentwork.
ACKNOWLEDGMENTS
We thank Kenta Kiuchi and Federico Carrasco forhelpful discussions. MS and SF thank Yukawa Institutefor Theoretical Physics, Kyoto University for their hos-pitality during the first corona pandemic time in Ger-many, in which this project was started. This work wasin part supported by Grant-in-Aid for Scientific Research(Grant Nos. JP16H02183 and JP20H00158) of JapaneseMEXT/JSPS. Numerical computations were performedon Sakura and Cobra clusters at Max Planck Computingand Data Facility. -1-0.5 0 0.5 1 -0.4 -0.2 0 0.2 0.4 B y x t =10 t =1 t =10 (sol) FIG. 7. Numerical solution of B y at t = 1 (dotted curve)and 10 (solid curve) with the analytic solution of Eq. (A1) at t = 10 (dashed curve). Appendix A: Test-bed simulations for resistiveMHD implementation
We here present the results for a suit of test-bed sim-ulations performed with our resistive MHD implementa-tion. We employed several test-bed problems introducedin Ref. [36] (see also Refs. [46–49]): the problems of self-similar current sheet, resistive shock tube, and resistiverotor. We also performed test simulations for the prop-agation of the electromagnetic wave packet in the flatspacetime and for a dynamo closure.
1. Self-similar current sheet
This test-bed problem was first proposed in Ref. [46]and subsequently employed in many references [36, 47–49]. This is the test to check whether the magnetic fieldcorrectly diffuses out with the special-relativistic resistiveMHD implementation in the limit that the gas pressureis much higher than the magnetic pressure. In the onedimensional problem in which all the quantities dependonly on the x coordinate, an approximate (nearly exact)solution is written, e.g., in the following form: B y ( x, t ) = C erf (cid:18) √ σ c x √ t (cid:19) , (A1)where erf denotes the error function and C is a constantwith C / π much smaller than the gas pressure. Here weset C = 1 as the choice of the units for simplicity. Notethat our notation for the basic MHD equations is differentfrom the previous papers because of the presence of thefactor 4 π in front of j µ in Eq. (2.3).Following Ref. [36], we employ the initial condition at t = 1 as ρ = 1, p = 50, E i = 1, u i = 0, B x = B z = 0,and B y = B y ( x, t = 1). The Γ-law equation of state, p = (Γ − ρε with Γ = 4 / σ c is set to8 -1.5-1-0.5 0 -0.4 -0.2 0 0.2 0.4 B y x σ c =infty σ c =100.0 σ c =10.00 σ c =1.000 σ c =0.001 ρ FIG. 8. Numerical solutions of ρ (upper panel) and B y (lowerpanel) for the shock tube test with σ c = 10 , 10 , 10, 1, and10 − . The lower panel plots B y (not ˆ B y ). be 100 (the resistivity is 1 / (4 πσ c )). The computationaldomain is set up as x = [ − . .
6] which is coveredby 201 uniform grid points. The simulation is performeduntil t = 10. Figure 7 plots B y at t = 1 and 10. For t = 10, we also plot the functional form of B y describedin Eq. (A1). Figure 7 illustrates that our implementa-tion can reproduce the solution well. In this setting themaximum relative error is ≈ . × − .
2. Resistive shock tube
This problem was first proposed in Ref. [47]. The ini-tial condition for this problem proposed in Ref. [36] is( ρ, p, v x , v y , v z , ˆ B x , ˆ B y . ˆ B z )= (cid:40) (1 . , . , . , . , . , . , . , .
3) for x < , (1 . , . , − . , − . , . , . , − . , .
5) for x ≥ , where ˆ B i = B i / √ π (this renormalization is necessaryto align the unit with the previous work) and v i = u i /u t .The initial electric field is given by E i = − (cid:15) ijk v j B k . TheΓ-law equation of state with Γ = 5 / σ c is chosen to be 10 − , 1, 10, 10 , and 10 . The com-putation is performed from t = 0 to 0 .
55 covering thecomputational domain of x = [ − . .
5] by 401 uni-form grid points. The results are shown for ρ and B y in Fig. 8 as in Ref. [36]. A small unphysical oscillationcould be seen for the curves of ρ and B y in the absence ofthe Kreiss-Oliger dissipation for evolving Maxwell’s equa-tions. However, this dissipation term with a small coef-ficient of 1 /
640 cures the unphysical oscillation. We findthat the numerical solutions are quite similar to thosefound in Ref. [36].
3. Resistive rotor
This is a two-dimensional problem with the coordi-nates ( x, z ) and the initial condition proposed in Ref. [36]is given in the following manner. Inside a radius of r = √ x + z ≤ .
1, a uniformly rotating medium with ρ = 10 and the uniform angular velocity of Ω = 8 . v x = − z Ω and v z = x Ω). Outsidethe circle of r = 0 .
1, on the other hand, we set ρ = 1and Ω = 0. For the entire region, the pressure andthe magnetic field are initially uniform as p = 1 and( ˆ B x , ˆ B y , ˆ B z ) = (1 , , E i = − (cid:15) ijk v j B k . The Γ-law equation ofstate with Γ = 4 / t = 0 to 0 . x = [ − . .
5] and z = [ − . .
5] and is covered by a uniformgrid of 401 ×
401 points. This problem is numericallysolved for σ c = 10 (nearly ideal MHD case) and 1. Theresults for the profiles of the pressure and z -componentof the electric field are shown in Fig. 9. By compar-ing these results with Fig. 3 of Ref. [36], we find a good(qualitative/semi-quantitative) agreement.
4. Evolving electromagnetic fields
To check that our implementation for evolving the elec-tromagnetic fields works well in a practical computa-tional domain, we also performed several test-bed sim-ulations in three spatial dimensions varying the value of σ c from 0 to high values. For σ c = 0 and σ c → ∞ inthe flat spacetime of u i = 0, it is easy to derive analyticsolutions for the electromagnetic equations, and thus, wecompare the results of the simulations with these analyticsolutions.For σ c = 0, the electromagnetic fields obey wave equa-tions in vacuum, and general solutions for the basic equa-tions are easily derived. For example, an axisymmetricgeneral solution for the dipole radiation is written as B r = 1 r ∂∂r (cid:18) f ( r + t ) − f ( r − t ) r (cid:19) cos θ,B θ = − r ∂∂r (cid:20) r ∂∂r (cid:18) f ( r + t ) − f ( r − t ) r (cid:19)(cid:21) sin θ, (A2)and B ϕ = 0. For this case, E ϕ is only non-zero compo-nent of E i and it is derived straightforwardly from thefield equation. Here, f ( u ) is an arbitrary regular func-tion, and for f ( u ) = − ue − u / , the initial condition be-comes B x = xze − r / ,B y = yze − r / ,B z = (2 − x − y ) e − r / , (A3)and E i = 0. Note here that we choose the width of thewave packet (i.e., unity) as the unit of the length in thisproblem.9 -0.4 -0.2 0 0.2 0.4 x -0.4-0.2 0 0.2 0.4 z x -0.4-0.2 0 0.2 0.4 z -0.4 -0.2 0 0.2 0.4 x -0.4-0.2 0 0.2 0.4 z -0.6-0.4-0.2 0 0.2 0.4 0.6 -0.4 -0.2 0 0.2 0.4 x -0.4-0.2 0 0.2 0.4 z -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 FIG. 9. The results of the resistive rotor problem for σ c = 10 (left) and 1 (right). The upper and lower panels show p and E z , respectively. -0.5 0 0.5 1 1.5 2 0 2 4 6 8 10 t =0 t =5 t =10 B z x -0.5 0 0.5 1 1.5 2 0 2 4 6 8 10 B z x FIG. 10. Left: Propagation of electromagnetic waves ( B z ) in the flat spacetime for σ c = 0 at t = 0, 5, 10, and 15. At t = 15, the amplitude is much smaller than unity entirely. The solid and dashed curves show the numerical and exact solutions,respectively. Right: Dispersion of electromagnetic waves in the flat spacetime for 4 πσ c = 50 at t = 0, 50, and 100. The solidand dashed curves show the numerical and exact solutions, respectively (both curves approximately agree). For x <
0, thesolution has the reflection symmetry for both panels (and we do not plot it).
We then evolve this wave packet in the domain of x = [ −
10 : 10], y = [ −
10 : 10], and z = [0 : 10] withthe grid spacing of dx = 0 .
1. The reflection-symmetricboundary condition is imposed on the z = 0 plane andan outgoing wave boundary condition is imposed on the outer boundaries. The left panel of Fig. 10 plots the nu-merical profiles of B z along the x axis at t = 0, 5, 10,and 15 (solid curves) together with the analytic solutions(dashed curves). Note that at t = 15, the wave packet al-ready propagated away from the computational domain.0 -2 -1
0 2 4 6 8 10 B y m a x t λ =2 λ =1 λ =1/2 λ =1/3 λ =2/7 λ =1/4 λ =2/9 λ =1/4, λ z =1/2fastest growth FIG. 11. The maximum value of B y for the steady dynamotest for σ c = 10 and α d = 0 . λ x = 2, 1, 1 /
2, 1 / /
7, 1 /
4, 2 / a = 0 .
1, and b = 0 (no perturbation in the z direction) and with λ x = 1 / λ z = 1 / a = 0 .
1, and b =0 .
01. Note that the fastest growing mode has the wavelengthof 1 / / λ in theinset of the plot implies λ x . The relative error of the numerical solution to the ana-lytic one, e.g., at t = 5, is of order 10 − with this setting(except for the outer boundaries). This illustrates thatour implementation can well follow the propagation ofthe electromagnetic waves.For large values of σ c , on the other hand, the field equa-tions relax to parabolic ones. For the parabolic equationswith the initial condition of Eq. (A3), a solution for thedipole magnetic field is written as B r = 2 (cid:18) t πσ c (cid:19) − / exp (cid:18) − r t/πσ c (cid:19) cos θ,B θ = − r (cid:18) t πσ c (cid:19) − / (cid:18) − r t/ πσ c (cid:19) × exp (cid:18) − r t/πσ c (cid:19) cos θ,B ϕ = 0 . (A4)We again evolve this wave packet in the domain of x = [ −
10 : 10], y = [ −
10 : 10], and z = [0 : 10] with dx = 0 . πσ c = 50. The right panel of Fig. 10plots the profiles of B z along the x axis at t = 0, 50, and100 (solid curves) together with the analytic solutions(dashed curves). With the chosen grid resolution, themaximum relative error size is ≈
4% at t = 100 (exceptfor the region near the outer boundaries). The largesterror is always located near the region at which B z = 0,and for other regions, the accuracy is much better.
5. Steady dynamo
Following Ref. [36], we also performed this simple testto check that our implementation works well in the pres-ence of a dynamo term. In this test-bed problem, again,we pay attention only to solving the electromagnetic fieldsetting u i = 0.First we assume that B x = E x = 0, and B y , B z , E y ,and E z are functions of exp( iωt − ikx ) where ω and k arethe angular frequency and the wave number. Then, thefollowing dispersion relation is derived: ω − πiσ c ω − k ± πσ c α d k = 0 . (A5)For this equation, we find that the exponential growthmode is present for the case that k < πσ c α d , and theresultant expression of ω is ω = 2 πiσ c (cid:16) − (cid:112) − ( k/ πσ c ) + kα d /πσ c (cid:17) , (A6)and for the fastest growing mode with k = 2 πσ c α d , thisbecomes ω = − πiσ c (cid:18)(cid:113) α − (cid:19) . (A7)Because the basic equation is linear in E i and B i , itis straightforward to extend this analysis in the multidi-mensional case as long as we focus only on the trans-verse component. That is, for the case that we ini-tially prepare two independent modes proportional to,e.g., exp( iωt − ik x x ) and exp( iωt − ik z z ) where k x and k z are the wave numbers, the stability of the system isdetermined simply by analyzing the dispersion relationsin the x and z direction independently. Taking into ac-count this fact, we perform two-dimensional simulationssetting the region of x = [ − z = [ − B x = − b cos( k z z ) , (A8) B y = a sin( k x x ) + b sin( k z z ) , (A9) B z = − a cos( k x x ) , (A10)and E i is determined by the ideal MHD condition. Here, a and b are constants. In this setting, for the case that2 λ i σ c α d >
1, the initial seed field should grow exponen-tially with time, where λ i := 2 π/k i ( i is x or z ) denotesthe wave length. Note that our setting of the compu-tational domain can follow the waves only with λ i ≤ σ c = 10 and α d = 0 .
2. For this setting, the fastest growing mode has λ i = 1 / λ i = 1 / λ i > / B y for σ c = 10 and α d = 0 . λ x = 2, 1, 1 /
2, 1 / /
7, 1 /
4, 2 / a = 0 .
1, and b = 0 (no perturbation in the z direction) and with λ x = 1 / λ z = 1 / a = 0 .
1, and b = 0 .
01. The growth rates of the unstable modes arewell captured in each simulation. This figure illustratesthat our implementation can derive the predicted resultsirrespective of the prepared initial conditions.1
0 0.1 0.2 0.3 0.4 0.5 E B ( e r g ) t (s) σ c =10 (1/s) σ c =10 (1/s) σ c =10 (1/s) σ c =10 (1/s) τ diss =0.01 s τ diss =0.1 s τ diss =1 sideal MHD FIG. 12. Evolution of the magnetic-field energy for the rem-nant massive neutron star initially with a purely toroidal mag-netic field for σ c = 10 , 10 , 10 , and 10 s − (solid curves).The dashed lines denote ∝ exp( − t/τ diss ) with τ diss = 0 . .
1, and 1 . σ c = 10 , 10 , and 10 s − , which are con-sistent with Eq. (2.59). The result by the ideal MHD imple-mentation is also shown by the dotted curve, which agreesapproximately with the curve for σ c = 10 s − . Appendix B: Resistive evolution of the massiveneutron star with purely toroidal magnetic field
In this section, we present the results for the resistiveMHD evolution of the merger remnant neutron star withpurely toroidal magnetic fields and demonstrate that ourimplementation can follow the resistive dissipation of thetoroidal magnetic field successfully.Figure 12 plots the evolution of the electromagneticenergy for the remnant massive neutron star, for whicha purely toroidal magnetic field of Eq. (2.54) is super-imposed at t = 0, in the axisymmetric resistive MHDsimulations with σ c = 10 , 10 , 10 , and 10 s − (solidcurves). The result of an ideal MHD simulation is showntogether (dashed curve). Since the toroidal magneticfield is simply superimposed and thus the initial condi-tion is not exactly in an equilibrium state, the electro-magnetic energy initially decreases by 10–20% even inthe absence of the resistivity during the early relaxationstage. The subsequent long-term gradual decrease of theelectromagnetic energy for the ideal MHD case wouldbe partly due to the numerical dissipation or diffusion. Note again that this evolution process is valid only inthe assumption of the axial symmetry. In general cases,non-axisymmetric instability such as Parker and Taylorinstability [42, 43] could occur and the evolution processcould be significantly modified. In these simulations the neutrino cooling is incorporated and thedensity and pressure profiles are modified during the simulation.This effect also affects the evolution of the magnetic field config-uration. -3 -2 -1 E B ( e r g ) t (s) Bh-R0-YBh-R0-NBh-R0-YLBh-R0-NLBh-Rl-YBh-Rw-YBh-Rl-YLBh-Rw-YL
FIG. 13. Evolution of the electromagnetic energy for twoideal MHD models Bh-R0-Y, Bh-R0-N and for two resistiveMHD models Bh-Rw-Y and Bh-Rl-Y with two different gridresolutions. The solid and dashed curves show the results for ∆x = 200 m and ∆x = 250 m, respectively. For the resistive MHD case, the magnetic field de-creases exponentially with time after the early relaxationstage. The short-dashed slopes denote the relation of ∝ exp( − t/τ diss ) where τ diss is the decay timescale, andare plotted for approximately fitting the curves for theresistive MHD results. It is found that the dissipationtimescale of the magnetic field is ≈ .
01, 0 .
1, and 1 . σ c = 10 , 10 , and 10 s − , which are consistent with thetimescales estimated by Eq. (2.59). For σ c = 10 s − ,the dissipation timescale would be ∼ s, and muchlonger than the simulation time. Therefore, the resultfor this case agrees approximately with that in the idealMHD simulation: Our ideal and resistive MHD imple-mentations can derive approximately the identical re-sults. Appendix C: Grid resolution dependence
As a convergence test, we performed numerical simula-tions with two lower grid resolutions for two ideal MHD(Bh-R0-Y and Bh-R0-N) and two resistive MHD (Bh-Rw-Y and Bh-Rl-Y) cases. Figure 13 is the same asFig. 1 but for the comparison in the evolution of theelectromagnetic energy for two different resolutions. Thesolid and dashed curves show the results for the casesthat the neutron star is resolved with the grid spacing of ∆x = 200 m and ∆x = 250 m, respectively.As found from Fig. 13, the growth of the magnetic-field strength by the winding effect is followed with thelower grid resolution fairly well. However, the ampli-fication of the magnetic-field strength is suppressed forthe lower resolution, and thus, the maximum magnetic-field energy becomes smaller. This effect is not appre-ciable for the case that the effects of the neutrino ir-radiation/heating and pair-annihilation are switched off.2However, in the presence of these neutrino effects, thesuppression is quite large: The maximum magnetic-fieldstrength for the lower resolution is by a factor of ∼ [1] M. Shibata and K. Hotokezaka, Ann. Rev. Nucl. Part.Sci. , 41 (2019).[2] B. P. Abbott et al., Phys. Rev. Lett. , 161101 (2017).[3] LIGO Scientic and VIRGO Collaborations, Astrophys. J. , L12 (2017).[4] B. D. Metzger and E. Berger, Astrophys. J. , 48(2012).[5] K. Hotokezaka and T. Piran, Mon. Not. R. Astron. Soc. , 1430 (2015).[6] S. A. Balbus and J. F. Hawley, Rev. Mod. Phys. , 1(1998).[7] R. Fern´andez and B. D. Metzger, Mon. Not. Royal As-tron. Soc. , 502 (2013).[8] B. D. Metzger and R. Fern´andez, Mon. Not. Royal As-tron. Soc. , 3444 (2014).[9] O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely, and H.-Th. Janka Mon. Not. Royal Astron. Soc. , 541 (2015).[10] D. M. Siegel and B. D. Metzger, Phys. Rev. Lett. ,231102 (2017): Astrophys. J. , 52 (2018).[11] S. Fujibayashi, K. Kiuchi, N. Nishimura, Y. Sekiguchi,and M. Shibata, Astrophys. J. , 64 (2018).[12] R. Fern´andez, A. Tchekhovskoy, E. Quataert, F. Fou-cart, and D. Kasen, Mon. Not. R. Astron. Soc. , 3373(2019).[13] A. Janiuk, Astrophys. J. , 2 (2019).[14] I. M. Christie et al., Mon. Not. R. Astron. Soc. , 4811(2019).[15] J. M. Miller et al., Phys. Rev. D , 023008 (2019).[16] S. Fujibayashi, M. Shibata, S. Wanajo, K. Kiuchi, K.Kyutoku, and Y. Sekiguchi, Phys. Rev. D , 083029(2020).[17] S. Fujibayashi, S. Wanajo, K. Kiuchi, K. Kyutoku,Y. Sekiguchi, and M. Shibata, Astrophys. J. , 122(2020).[18] R. Fern´andez, F. Foucart, and J. Lippuner, submitted toMon. Not. R. Astron. Soc. (arXiv: 2005.14208).[19] P. M¨osta, D. Radice, R. Haas, E. Schnetter, and S.Bernuzzi, Astrophys. J. , L37 (2020).[20] M. Shibata, K. Taniguchi, and K. Ury¯u, Phys. Rev. D , 084021 (2005).[21] S. L. Shapiro, Astrophys. J. , 397 (2000).[22] M. Shibata, Y.-T. Liu. S. L. Shapiro, and B. C. Stephens,Phys. Rev. D , 104026 (2006).[23] L. Sun, M. Ruiz, and S. L. Shapiro, Phys. Rev. D ,064057 (2019).[24] K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, andT. Wada, Phys. Rev. D , 041502 (2014); K. Kiuchi,P. Cerda-Duran, K. Kyutoku, Y. Sekiguchi, and M. Shi-bata, Phys. Rev. D , 124034 (2015); K. Kiuchi et al.,Phys. Rev. D , 124039 (2018).[25] S. Fujibayashi, Y. Sekiguchi, K. Kiuchi, and M. Shibata, Astrophys. J. , 114 (2017).[26] T. G. Cowling, Mon. Not. R. Astron. Spoc. , 39(1933).[27] M. Shibata, Numerical Relativity (World Scientific, Sin-gapore, 2016).[28] M. Shibata and T. Nakamura, Phys. Rev. D ,5428(1995): T. W. Baumgarte and S. L. Shapiro, Phys.Rev. D , 024007(1998).[29] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-chower, Phys. Rev. Lett. , 111101 (2006): J. G. Baker,J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter,Phys. Rev. Lett. , 111102 (2006).[30] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W.Tichy, and B. Br¨ugmann, Phys. Rev. D , 084057(2013).[31] M. Alcubierre, et al., Int. J. Mod. Phys. D , 273 (2001).[32] M. Shibata, Prog. Theor. Phys. (2000), 325; M. Shi-bata, Phys. Rev. D (2003), 024033; M. Shibata andY. Sekiguchi, Prog. Theor. Phys. , 535 (2012).[33] S. Banik, M. Hempel, and D. Bandyophadyay, Astro-phys. J. Suppl. Ser. , 22 (2014).[34] F. X. Timmes and F. D. Swesty, Astrophys. J. Suppl. , 501 (2000).[35] A. Brandenburg and K. Subramanian, Phys. Rep. ,1 (2005).[36] N. Bucciantini and L. Del Zanna, Mon. Not. R. Astron.Soc. , 71 (2013).[37] M. Shibata and Y. Sekiguchi, Phys. Rev. D , 044014(2005).[38] C. R. Evans and J. F. Hawley, Astrophys. J. , 659(1988).[39] D. Duez, Y.-T. Liu, S. L. Shapiro, M. Shibata, and B. C.Stephens, Phys. Rev. Lett. , 031101 (2006).[40] M. Shibata, M. D. Duez, Y.-T. Liu, S. L. Shapiro, andB. C. Stephens, Phys. Rev. Lett. , 031102 (2006).[41] M. D. Duez, Y.-T. Liu, S. L. Shapiro, M. Shibata, andB. C. Stephens, Phys. Rev. D , 104515 (2006).[42] E. N. Parker, Astrophys. J., , 49 (1955).[43] R. J. Tayler, Mon. Not. R. Astron. Soc., , 365 (1973).[44] K. Kiuchi, S. Yoshida, and M. Shibata, Astron. and As-trophys. , A30 (2011).[45] V. Duez, J. Braithwaite, and S. Mathis, Astrophys. J., , L34 (2010).[46] S. S. Komissarov, Mon. Not. R. Astron. Soc. , 995(2007).[47] M. Dumbser and O. Zanotti, J. Comput. Phys., ,6991 (2009).[48] C. Palenzuela, L. Lehner, O. Reula, and L. Rezzolla,Mon. Not. R. Astron. Soc. , 1727 (2009).[49] Q. Qian, C. Fendt, S. Noble, and M. Bugli, Astrophys.J.834