Long-term evolution of a merger-remnant neutron star in general relativistic magnetohydrodynamics I: Effect of magnetic winding

Long-term ideal and resistive magnetohydrodynamics (MHD) simulations in full general relativity are 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 simulation is performed in axial symmetry and without considering dynamo effects as a first step. In the ideal MHD, the differential rotation of the remnant neutron star amplifies the magnetic-field strength by the winding in the presence of a seed poloidal field until the electromagnetic energy reaches $\sim 10\%$ of the rotational kinetic energy, $E_{\rm kin}$, of the neutron star. The timescale until the maximum electromagnetic energy is reached depends on the initial magnetic-field strength and it is $\sim 1$ s for the case that the initial maximum magnetic-field strength is $\sim 10^{15}$ G. After a significant amplification of the magnetic-field strength by the winding, the magnetic braking enforces the initially differentially rotating state approximately to a rigidly rotating state. In the presence of the resistivity, the amplification is continued only for the resistive timescale, and if the maximum electromagnetic energy reached is smaller than $\sim 3\%$ of $E_{\rm kin}$, the initial differential rotation state is approximately preserved. In the present context, the post-merger mass ejection is induced primarily by the neutrino irradiation/heating and the magnetic winding effect plays only a minor role for the mass ejection.


I. INTRODUCTION
Theoretical exploration for the evolution of the merger remnants of binary neutron stars has become a hot topic in this decade [1], because such systems can shine as an electromagnetic counterpart of gravitational waves emitted in their inspiral stage and the signals bring us a variety of information for the nature of the merger process and neutron stars that cannot be obtained only by the gravitational-wave observation, as the first observation of the binary neutron star merger shows [2,3]. There are a variety of the possibilities for the remnants of the binary neutron star mergers [1] and for the corresponding electromagnetic counterparts [4,5]. Irrespective of the possibilities, the key for the strong electromagnetic emission is the mass ejection from the remnant. Thus, one important aspect of the theoretical study is to clarify the properties of the ejecta such as total mass, typical velocity, and typical elements.
The promising component of the remnant that induces the mass ejection is a disk (or torus) surrounding the central compact object, which is either a massive neutron star or a black hole. The reason for this is that the disk is differentially rotating, approximately with the Keplerian rotational profile, and is likely to be magnetized because the disk matter stems from neutron stars. In this situation, the magnetorotational instability (MRI) [6] is activated, and as a result, a turbulence would be developed, enhancing the turbulence viscosity. The resulting viscous heating and angular momentum transport inevitably enhance the activity of the disk, and even-tually, the mass ejection is induced [7][8][9][10][11][12][13][14][15][16][17][18]. In addition, the amplified magnetic field could eventually constitute a aligned poloidal magnetic field along the rotational axis and this could further enhance the mass ejection efficiency [10,14,15]. With these considerations, numerical simulations have been extensively performed in the last decade to clarify the details of the mass ejection from the accretion disks.
Although many works have been done for exploring the evolution of the accretion disks formed as remnants of the neutron-star merger, there are only a few works for directly evolving the remnant massive neutron star for a long timescale of 1 s [11,17] (more specifically, no long-term magnetohydrodynamics (MHD) simulation has been done: but see [19] for a short-term simulation). A unique property of the remnant massive neutron stars is found in their angular velocity profile, Ω( ), where is the cylindrical radius: Reflecting the nearly irrotational velocity field of the pre-merger stage of two neutron stars, the inner region of the remnant neutron stars should be differentially rotating with the angular velocity that increases with the cylindrical radius (dΩ/d > 0; e.g., Ref. [20]). This situation is in contrast to that for the disk rotating around a central object, which approximately has a Keplerian profile, i.e., the angular velocity decreases with the increase of the cylindrical radius, and thus, unstable to the MRI. For the inner region of the remnant neutron stars (in particular for the central part of it) for which dΩ/d > 0, the MRI does not play any role.
However, irrespective of the sign of dΩ/d , the wind-arXiv:2102.01346v1 [astro-ph.HE] 2 Feb 2021 ing of the magnetic-field lines occurs in the presence of a seed poloidal magnetic field together with differential rotation. In contrast to the MRI for which the growth of the magnetic-field strength exponentially proceeds, the (toroidal) magnetic-field strength increases only linearly with time by the winding. Nevertheless, the growth rate is not still negligible in the presence of the rapid rotation because the amplification factor is approximately written as Ωt (e.g., Refs. [21][22][23]) and Ω ∼ 5000 rad/s (see Sec. III) for the merger remnants. Thus in 0.2 s, the magnetic field can be amplified by ∼ 10 3 times unless strong resistive processes are present. That is, even in the presence of a seed poloidal magnetic field typically of 10 13 -10 14 G (which may be a conservative value for the merger remnants [24]), the (toroidal) magnetic field could be amplified to the order of 10 16 -10 17 G in ∼ 0.2 s after the merger. The resulting magnetic pressure could be a substantial fraction of the matter pressure, and hence, the MHD effect could play an important role for the late time evolution of the remnant neutron star.
Motivated by this consideration, we perform general relativistic MHD (GRMHD) simulations for a massive neutron star formed as a remnant of equal-mass binary neutron star mergers: As in our series of papers [11,17,25], we employ a numerical result of a simulation for binary neutron star mergers as the initial condition and perform an axisymmetric simulation. We solve Einstein's equation, GRMHD equations, and neutrinoradiation hydrodynamics equations altogether. Both the ideal and resistive MHD simulations are performed. We here note that the origin of the resistivity of the neutron stars is not certain. The purpose here is to simply show that in the presence of a hypothetical resistive MHD process, the evolution of the remnant massive neutron stars can be significantly modified.
In the assumption of axial symmetry, the poloidal magnetic field is not amplified even when the toroidalfield strength is significantly enhanced, due to the antidynamo property [26]. As a first step toward more detailed exploration of this topic, we do not consider any dynamo effects for the growth of the magnetic-field strength in this paper. Thus, we focus only on the growth of the toroidal magnetic field and its effect on the evolution of the remnant neutron star. The effect of the dynamo is planned to be studied in the subsequent work.
The paper is organized as follows: In Sec. II, we summarize the basic equations employed in the present numerical simulation paying particular attention to the GRMHD equations and our method for solving resistive MHD equations. Section III presents numerical results of the simulations paying particular attention to the growth of the magnetic-field energy by the winding and the resulting mass ejection. Section IV is devoted to a summary. In Appendix A, a number of the results of test simulations for our newly developed resistive MHD implementation are presented. Throughout this paper, we use the geometrical units of 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 units in the following sections). Latin and Greek indices denote the spacetime and space components, respectively. In Sec. II, we suppose to use Cartesian coordinates for the spatial components whenever equations are written.

II. BASIC EQUATIONS AND METHOD FOR NUMERICAL COMPUTATIONS
A. Brief summary In this work, we numerically solve Einstein's equation, ideal or resistive MHD equations, evolution equations for the lepton fractions including the electron fraction, and (approximate) neutrino-radiation transfer equations. Except for the MHD parts (see the following subsections for details), the numerical implementation is the same as that in our latest viscous-hydrodynamics work [16,17]: Einstein's equation is solved using the original version of the Baumgarte-Shapiro-Shibata-Nakamura formalism [28] together with the puncture formulation [29], Z4c constraint propagation prescription [30], and 5thorder Kreiss-Oliger dissipation. The axial symmetry for the geometric variables is imposed using the cartoon method [31,32] with the 4th-order accuracy in space. For evolving the lepton fractions, we take into account electron and positron capture, electron-positron pair annihilation, nucleon-nucleon bremstrahlung, and plasmon decay [11,17]. We employ a tabulated equation of state based on the DD2 equation of state [33] for a relatively high-density part and the Timmes (Helmhoitz) equation of state for the low-density part [34]. In this tabulated equation of state, thermodynamics quantities such as ε, P , and h are written as functions of ρ, Y e , and T where ε, P , h(= c 2 + ε + P/ρ), ρ, Y e and T are the specific internal energy, pressure, specific enthalpy, rest-mass density, electron fraction, and matter temperature, respectively. We choose the lowest rest-mass density to be 0.1 g/cm 3 in the table, and the atmosphere density for ρ * := ρu t √ −g in the hydrodynamics simulation is chosen to be 10 3 g/cm 3 in the central region of r ≤ 100 km and it is decreased down to 1 g/cm 3 with the dependence of ∝ r −3 in the outer region (i.e., for the far region it is 1 g/cm 3 ). Here u µ and g denote the four velocity of the fluid and the determinant of the spacetime metric g µν , respectively. ρ * obeys the continuity equation in the form of First we write down the equations for the electromagnetic fields, which are derived from Maxwell's equations: Here ∇ µ is the covariant derivative with the respect to g µν , F µν is the electromagnetic tensor, j µ is the current four vector, and * F µν is the dual of F µν defined by * with µναβ the Levi-Civita tensor of txyz = √ −g. Using the unit timelike vector normal to the spacelike hypersurface, n µ , we define the electric and magnetic fields by and thus, the electromagnetic tensor and its dual are written as Then Maxwell's equations are written in terms of E µ and B µ as follows: For the constraint equations, (2.10) D k B k = 0, (2.11) and for the evolution equations, where ρ e := −j a n a andj i := γ i k j k with γ ij the spatial metric defined by γ µν := g µν + n µ n ν . D k denotes the covariant derivative with respect to γ ij , β i the shift vector, K the trace of the extrinsic curvature K ij , and L β denotes the Lie derivative with respect to β i . We numerically solve the evolution equations in the forms [27] 14) where nates, the definition should be appropriately modified by excluding the contribution of the coordinates in γ). We note (−g) = α 2 γ. In this notation, the constraint equations are written in the simple forms as (2.17) To close the equations, we in general need the Ohm's law, for which we here write as is the charge density observed in the frame comoving with 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 of Eq. (2.18) denotes a dynamo term (in the simplest version) 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 a phenomenological one and we do not consider it in this paper (i.e., α d = 0) except for Appendix A 5. Note that F µν u ν and − * F µν u ν (=: b µ ) denote the electric and magnetic 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 in terms of the magnetic field in the frame comoving with the fluid, b µ , as Then, it becomes trivial that the corresponding electric field, F µν u ν , vanishes. Then, the condition of F µν u ν = 0 with Eq. (2.8) gives the electric field, E i , as and thus, E µ u µ (= E i u i ) = 0.
In the ideal MHD case, the current is not determined from Eq. (2.18) due to the fact that σ c → ∞. Instead, 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 electromagnetic field, which determines the magnetic-field evolution. Using Eq. (2.21), this equation is written in the well-known form as In the ideal MHD, the energy-momentum tensor is written as Here, we have the relations as with B 2 = B k B k , and thus, the energy-momentum tensor is also written in terms of B k . For evolving the MHD equations, we define S 0 := √ γ T µν n µ n ν and S i := − √ γ T µν n µ γ νi , which are written as Their evolution equations are where P tot := P + b 2 /8π 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 Using this, the normalization relation of u µ is written in terms of γ ij , ρ * , S i , and B i as Hence, for given (evolved) values of γ ij , ρ * , S i , S 0 , and B i together with a given equation of state, this equation together with Eq. (2.27) constitute simultaneous equations for h and w, and thus, they are used for the socalled primitive recovery procedure in the ideal MHD. Our numerical approach for the primitive recovery when employing tabulated equations of state is essentially the same as that in the viscous hydrodynamics case [25].

D. MHD in general cases
In the general MHD case, not only the magnetic field but also the electric field are explicitly included in the energy-momentum tensor as T µν = ρhu µ u ν + P g µν where E 2 = E k E k . We note that in the ideal MHD case, using Eq. (2.21) together with the fact that E i / E k E k , B i / B k B k , and γ ij u j / γ kl u k u l can constitute the orthonormal bases of the three space, the energy-momentum tensor of Eq. (2.23) is derived in a straightforward manner. For the general MHD case, S 0 and S i are written as and their evolution equations are where S ij = ρhu i u j + P γ ij Using Eq. (2.35), the normalization relation for u µ is written as Thus, for the resistive MHD, this equation together with Eq. (2.34) constitute the simultaneous equations for h and w for given (evolved) values of γ ij , ρ * , S i , S 0 , B i , and E i , and are used for the primitive recovery procedure. Practically, we solve the MHD equations (both in ideal and resistive MHD) using the cylindrical coordinates ( , ϕ, z) in this work. For this procedure, we write the equations for S 0 , S y , ρ * , and the conservation equations for the lepton fraction in a conservative form as in our viscous hydrodynamics simulations [16,17]. This is in particular important to guarantee the conservation of the rest mass and angular momentum in numerical computation.

E. Numerical methods for solving Maxwell's equations
We solve the ideal MHD equations mostly in the same method as described in Ref. [37]: We evolve the hydrodynamics equations for S 0 , S i , ρ * , and the lepton fraction using a high-resolution shock capturing scheme (a 3rdorder upwind scheme). In the new implementation, one improvement is made for solving the induction equation on the evolution of the magnetic field, B k . The previous 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 evolving S 0 , S i , and ρ * . However, we were aware of the fact that this method was so diffusive that the magnetic-field energy was spuriously decreased with a short timescale within 100 ms with the typical grid resolution of the grid spacing ∼ 200 m (we use the grid resolution of DD2-135M of Ref. [17] for most of the present simulations and that of DD2-135L of Ref. [17] for 4 test simulations). Thus, for a long-term simulation with the duration more than seconds, this scheme is practically useless (although a sophisticated high-order scheme may solve this problem). In the new implementation, we solve the induction equation for B x and B z by a high-resolution upwind scheme, which is the same as for evolving S 0 , S i , and ρ * , and the induction equation for B y by the 4th-order centered finite differencing with the operation of the 5thorder Kreiss-Oliger dissipation. To approximately preserve the divergence-free condition of Eq. (2.17), a divergence cleaning is introduced (in the same manner as in the resistive MHD scheme: see below). We find that in this method, the stable numerical evolution is feasible and also spurious oscillation is suppressed with the reasonable choice for the coefficient of the Kreiss-Oliger term. By this procedure, the diffusive evolution of the magnetic field is suppressed in a much better manner than in our previous implementation. The maximum error size of the divergence-free condition, which is defined by |(∂ i B i )∆x/B z | where ∆x is the grid spacing that covers remnant neutron stars, remains to be ∼ 10 −5 for the entire simulation time in the present numerical resolution.
In the resistive MHD, we also employ the same shock capturing scheme as that employed in the ideal MHD for the hydrodynamics part that evolve S 0 , S i , ρ * , and the lepton fraction. On the other hand, Maxwell's equations are solved in a different procedure. In addition, we change the method of the time evolution depending on the magnitude of the conductivity. In the following, we describe the methods for the low-and highconductivity cases separately, focusing only on the case that the timescale of 1/(4πσ c ) is shorter than the time step interval of the numerical simulation, ∆t: A partially implicit scheme is employed for appropriately handling the equation of the electric field for the high values of σ c .
For the low-conductivity case, we evolve the electromagnetic equations using the method that is used for evolving geometrical variables and incorporating a divergence cleaning prescription. Specifically, we rewrite Eqs. (2.14) and (2.15) in the forms: where φ B is a new auxiliary variable associated with the divergence cleaning, κ is a constant, and J i − Qβ i is written as We notice that in our scheme, we always evaluate Q by ∂ k E k /4π.
In numerical implementation, the transport terms in Eqs. (2.40) and (2.41) are evaluated by the 4th-order centered finite differencing, and the transport term in Eq. (2.42) is evaluated by the 4th-order upwind scheme as done for the geometric variables. Whenever we evolve E i , B i , and φ B , we operate the Kreiss-Oliger dissipation procedure in the same manner as we do for the geometric variables.
To evolve the system, we basically use the 4th-order Runge-Kutta scheme but for E i we partially employ an implicit scheme. This is implemented in the following manner: In Eq. (2.40), we move the term of 4πσ c αwA i j E j that appears in the current term to the left-hand side and write the left-hand side of the evolution equation in the form where E i 0 denotes the values of E i in a previous time-step and M i j = (1 + ζ)δ i j − ζū i u j /w 2 with ζ = 4πσ c αw∆t. Here, the inverse of M i j is simply written as Thus, the updated values of E i is easily obtained in each Runge-Kutta time step as where F j T denotes the term associated with the transport term, which is evaluated in the explicit manner. We suppose that α d is small (or zero) and thus the last term of Eq. (2.46) is not included for the implicit prescription, although it is straightforward to include it.
As we illustrate in Appendix A, this implementation enables to solve several test-bed problems successfully even for the case that σ c is large such that ζ 1. The electromagnetic fields are also evolved stably even in the absence of the upwind scheme.
For the high-conductivity case close to the ideal MHD case (specifically for ζ ≥ 100 in our present implementation), we modify this procedure. The reason for the modification is that the equation of the magnetic field approximately reduces to the induction equation (with a weak diffusion term: see below) as Here in Eq. (2.47) we omit the divergence cleaning term for simplicity, although we include it in the practical simulations. As a natural consequence of reducing approximately to a transport equation, for the case of the long resistive dissipation timescale, the centered finite differencing could introduce a long-term growing unstable mode in numerical computation, in particular for the evolution of the poloidal magnetic field leading to the violation of divergence-free condition of ∂ k B k = 0 (note that this condition is written only by the poloidal-field component in the axially symmetric case). Below, we first derive Eq.
This is easily understood that the stiff term associated with σ c in the right-hand side of the evolution equation for E i enforces the exponential decay of the initial electric field. Then, Eq. (2.46) is rewritten into the form denotes the term associated with the dynamo. Hence, the terms associated with the transport for E i (see Eq. (2.41)) is rewritten as Thus, Eq. (2.47) is derived from Eq. (2.41) (here we suppose again that α d is small).
As in the ideal MHD case, the transport part of Eq. (2.47) should be evaluated by using an upwind scheme for the numerical stability. On the other hand, the last term of Eq. (2.50) constitutes a diffusion term in Eq. (2.41), which can be evaluated by the centered finite differencing. In our numerical implementation, the transport term is handled in the same manner as in the ideal MHD case, while the diffusion term is handled by the 4th-order centered finite differencing. In addition, the divergence cleaning and 5th-order Kreiss-Oliger dissipation are implemented as in the relatively low conductivity case.
Before closing this subsection, we briefly comment on our artificial prescription for handling the region in which the electromagnetic energy density is much larger than the rest-mass energy density, because such regions often cannot be solved accurately in the MHD simulation. In the present paper, if the ratio of B 2 /(4πρc 2 ) (in the ideal MHD case) or (E 2 + B 2 )/(8πρc 2 ) (in the resistive MHD case) is larger than 10 3 , we set u i = 0 artificially. In addition, if w(= αu t ) exceeds 5 or ε exceeds 10 3 after the primitive-recover process, we also artificially set u i = 0 and ε is determined from the condition that the temperature is kT = 0.1 MeV where k is the Boltzmann constant. The reason for this is that in these high-electromagneticenergy cases, the primitive recovery often fails. The region with these extreme conditions appears for the case that the magnetic winding proceeds significantly and as a result a strong outflow is driven from the neutron-star surface toward a low-density atmosphere. In the absence of these artificial prescriptions, the computation often crashed in our current implementation. Since we give up resolving the region with w ≥ 5, we cannot explore a very-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 initial condition for the matter field is supplied from the result 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 each neutron-star mass 1.35M . As touched in Sec. II E, four simulations are also performed using a lower-resolution model, DD2-135L, to check the dependence of the results on the grid resolution.
We then superimpose electromagnetic fields for which the energy density is much smaller than the internal and kinetic energy density of the fluid. In this work, we initially give a poloidal magnetic field written in the form with δ A = 10 −5 and P max is the maximum pressure. A 0 determines the initial field strength. With this setting, the magnetic field is present for a high-density region with ρ 10 11.5 g/cm 3 . The toroidal magnetic field is set to be zero (B y = 0) initially. The electric field is determined by the ideal MHD condition (2.21). We also performed two ideal MHD simulations with δ A = 10 −2 (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 remnant massive neutron star for which the magnetic-field effect appears in the most remarkable manner in the present context (i.e., in the axisymmetric simulation with no dynamo effect). Indeed, it is found that the results on the evolution of the electromagnetic energy and angular velocity profile depend only weakly on the value of δ A . Thus in the following, we focus only on the results of δ A = 10 −5 . As Eq. (2.53) shows, we pay attention to the MHD effect only for the high-density region, i.e., for the remnant neutron star, and do not pay attention to the MHD effect in the disk. Because the angular velocity of the disk decreases along the cylindrical-radius direction, the MRI would play a crucial role for its evolution. For exploring the evolution, we need a long-term non-axisymmetric simulation that can capture the effects of turbulent motion and resulting effective viscosity induced. This is beyond the scope of this paper.
We also performed simulations initially with a purely toroidal magnetic field of Here the dependence on the coordinates, z, stems from the 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 decay with the timescale determined by σ c (cf. Eq. (2.59) below). In Appendix B, we show that our implementation can derive the expected numerical results. Table I lists the initial conditions (with the purely poloidal magnetic field).B init and E B denote the maximum value of B z γ −1/3 (= B z γ −1/6 ) and the electromagnetic energy at the initial state. Here, the electromagnetic energy is defined by Note that E 2 is much smaller than B 2 . Since the volume of the neutron star is of the order of 10 19 cm 3 , the volume averaged magnetic-field strength is initially ∼ 10 15 G for E B ∼ 10 48 erg. The labels, Bv, Bh, Bm, and Bl, in the model name of Table I refer to the initial magnetic-field strength (very high, high, medium, and low). Because the remnant neutron star which we employ is differentially rotating, with this setting, the toroidal magnetic-field strength is initially increased linearly with time by the winding effect as (e.g., Refs. [21][22][23]) where B 0 is the initial local magnitude of the component of the magnetic field. This growth continues until the magnetic-field energy increases to ∼ 10% of the rotational kinetic energy, E kin , of the neutron star. Thus the winding timescale is defined approximately by where Ω max and Ω 0 denote the maximum value of the angular velocity and angular velocity at the center, respectively. E B,0 is the initial magnetic-field energy, and E kin is defined by For the model employed in this paper, the initial value of E kin is ≈ 1.16 × 10 53 erg. In the estimate of Eq. (2.57), we supposed that the region of the maximum angular velocity would be located in an outer region of the remnant massive neutron star (see, e.g., Ref. [16]) and it governs the winding and the growth of the electromagnetic energy. Also, we assumed that a part of the poloidal magnetic field for which the electromagnetic energy is ∼ 10% of E B,0 contributes to the winding. In the case of the finite value of σ c , the electromagnetic fields are dissipated. The dissipation timescale is approximately written as where denotes the typical variation scale of the magnetic field. Because we do not know the typical size of σ c , we pay attention to the cases for which τ dis is comparable to τ wind , i.e., σ c = 10 7 -10 9 s −1 .
For the case of τ dis τ wind , obviously, the resistive dissipation does not play a role. By contrast, for τ dis τ wind , the magnetic field is diffused out before the magnetic winding significantly amplifies the magneticfield strength. Also, the winding effect does not directly determine the final magnetic-field profile for τ dis τ wind . I. Initial condition and set-up for the numerical simulation. We list the maximum value ofBinit = B z γ −1/3 , the magnetic-field energy, and the value of σ −1 c in units of milliseconds. The last column shows the on or off of the irradiation/heating and pair-annihilation of neutrinos. For all the initial conditions, the total baryon mass is M * = 2.95M , the gravitational mass is M = 2.64M , the total rotational kinetic energy is E kin ≈ 1.16 × 10 53 erg, and the total angular momentum is J = 4.65 × 10 49 g cm 2 /s. 46 10 9 On Bm-Rm-Y 1.1 × 10 15 5.9 × 10 46 10 8 On Bm-Rh-Y 1.1 × 10 15 5.9 × 10 46 10 7 On Bl-Rl-Y 5.6 × 10 14 1.5 × 10 46 10 9 On This point is understood by the following model for the evolution of B y , for which the evolution equation is approximately written as where ∆ denotes the Laplacian operator (different from ∆). Here, we picked up only the terms related to the winding and resistive dissipation. For this equation, we consider a toy model in which the first term in the righthand side is not time-varying as B i ∂ i Ω = F (x i ) =: η∆F 2 (x i ). Then, it is easily found that for t τ dis , the magnetic winding approximately determines the solution as B y ≈ tF (x i ). On the other hand, for t τ dis , the asymptotic solution is B y ≈ −F 2 (x i ). Thus, the winding history is not reflected, although the order of magnitude of F 2 is written as ∼ Ωτ dis |B i |, which reflects the winding in t ≤ τ dis . We note that near the rotational axis, F ∝ 3 , and hence, F 2 ∝ 5 . Thus, even if B y is positive during the winding, the final relaxed value for it can be negative.

III. EVOLUTION OF A REMNANT NEUTRON STAR
A. Summary of the evolution Table I lists the models for which we performed simulations. The initial electromagnetic-field energy is varied among 1.5 × 10 46 -9.4 × 10 47 erg. With this setting, the magnetic winding timescale is approximately 0.1-1 s.
The simulations are performed for the ideal MHD models (referred to as R0 in the model name) and for the resistive MHD models with σ c = 10 7 , 10 8 , 10 9 , and 10 11 s −1 (referred to as Rh, Rm, Rl, and Rw in the model name). For σ c = 10 11 s −1 , the dissipation timescale by the resistivity is ∼ 100 s, which is much longer than the winding timescale in our setting. Thus, this model is employed to check that the ideal MHD results are approximately reproduced in the resistive MHD implementation.
For most of the models, neutrino effects including the irradiation/heating and pair annihilation are taken into account as in our previous papers [17,25]. To clarify the importance of the neutrino effects on the evolution of the massive neutron stars and on the mass ejection, for the ideal MHD case, we perform two simulations in which the effects of both the neutrino irradiation/heating and neutrino pair annihilation are switched off (models Bh-R0-N and Bm-R0-N: in the presence of the neutrino effects, the model is referred to with the label Y).
As we touched in Secs. II E and II F, numerical simulations were performed with two 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. By comparing the results of two different resolutions, we find that the numerical results depend weakly on the grid resolution (see Appendix C). Specifically, for the lower grid resolution, the amplification of the magnetic-field strength is suppressed, and thus, the maximum magnetic-field energy becomes smaller, in particular in the presence of the neutrino irradiation/heating and pair-annihilation effects. For the resistive MHD simulation, we also find that the dissipation is spuriously enhanced for the lower grid resolution. However, besides these quantitative effects, the evolution process of the neutron star is not modified qualitatively by changing the grid resolution.
In the present context, the typical MHD process in the remnant neutron star of neutron-star mergers in the presence of an initial seed poloidal magnetic field is governed by the magnetic winding effect, for which the growth timescale of the toroidal magnetic-field energy is written approximately by Eq. (2.57). When the electromagnetic energy exceeds a few % of the rotational kinetic energy of the remnant neutron star (for this case E kin ≈ 1.16×10 53 erg), the magnetic braking starts playing an important role. As a consequence, the angular momentum in the neutron star is redistributed and the differentially rotating configuration of the angular velocity is enforced to approach a (approximately) rigidly rotating one. Then, the angular velocity in the central region exceeds ∼ 5000 rad/s (i.e., the rotational period is ∼ 1.2 ms). After the toroidal magnetic-field strength is significantly increased by the winding, the mass outflow from the neutron-star surface is enhanced by the increased magnetic pressure (although the mass ejection in the present context is driven primarily by the neutrino irradiation/heating: see Sec. III B). After the magnetic braking works significantly, the increase of the angular velocity in the central region is decelerated. In  addition, the magnetic-field profile in the neutron star is modified significantly by the mass outflow (and resulting magnetic-flux escape) process. The detail on these findings is described in Sec. III A 1.
In the presence of only low resistivity (i.e., a high value of σ c ), the evolution process of the remnant neutron star is essentially the same as in the ideal MHD case. In the presence of a high resistivity with which the dissipation timescale of the magnetic field is shorter than the winding timescale, not only the winding but also the magneticfield dissipation plays an important role. For such cases, the magnetic winding does not occur significantly and the magnetic-field effect on the angular velocity profile is minor (see Sec. III A 2 for details). Figure 1 shows the electromagnetic energy as a function of time. The left panel plots the results for all the ideal MHD models together with one resistive model of a high conductivity (Bh-Rw-Y) and the right panel plots those for all the resistive MHD models. In the ideal MHD cases, the electromagnetic energy, E B , is initially increased by the magnetic winding effect by several orders of magnitude. In this early stage, E B is approximately proportional to t 2 irrespective of the models (including the resistive MHD models except for the case of σ c = 10 7 s −1 for which the resistive dissipation timescale is quite short; the order of 10 ms).

Ideal MHD
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 2 ), indicating that the winding effect is weaken and the magnetic braking becomes important. However, the increase of E B still continues until it reaches ∼ 10% of E kin , and then, the growth of the electromagnetic field is saturated. The satura-tion energy is slightly smaller for the smaller initial field strength. This is primarily due to the fact that more magnetic fluxes are escaped associated with the mass outflow for longer-term evolution, and may be partly due to the numerical dissipation, i.e., in the longer-term evolution the damage by the numerical dissipation is more serious: We note that the simulations are always performed for ∼ 10 3 rotational period of the neutron star.
The saturated energy is slightly higher in the absence of the neutrino irradiation/heating and pair-annihilation effects. Our interpretation for this is that in the absence of these effects, the mass outflow from the neutron star is suppressed (cf. Sec. III B), resulting escape of the magnetic flux (associated with the mass outflow) is also suppressed, and as a result, the amplification of the electromagnetic field proceeds to a higher level.
As we already mentioned briefly, the mass ejection is driven primarily by the neutrino irradiation/heating (see Sec. III B for details). However, with the increase of the electromagnetic energy toward the saturation, the mass ejection from the neutron star is slightly enhanced. This is due to the magnetic pressure primarily resulting from the significantly amplified toroidal magnetic field. By this enhancement of the outflow, a part of the magnetic flux escapes from the neutron star, and as a result, the total electromagnetic energy decreases. However, the decrease rate becomes eventually low (E B eventually shows the oscillatory behavior) and the electromagnetic energy in average appears to relax to ∼ 10 51 erg, i.e., ∼ 1% of the rotational kinetic energy of the neutron star. We note that the rotational kinetic energy does not change significantly by the magnetic braking effects because it is always larger by more than one order of magnitude than the electromagnetic energy, and thus, the magnetic-field effect is minor. Instead, the rotational kinetic energy slightly increases with time because the neutron star contracts due to the emission of neutrinos of total energy of ∼ 10 53 erg, and hence, the neutron star slightly spins up in the absence of the electromagnetic effect, as shown in Ref. [17].
For the nearly ideal MHD case, i.e., for the resistive MHD case with a tiny resistivity (high conductivity, σ c = 10 11 s −1 ; model Bh-Rw-Y), the curve of E B is similar to the corresponding ideal MHD model (Bh-R0-Y) as found in the left panel of Fig. 1. This is natural because for this model, the resistive dissipation timescale is quite long ∼ 100 s (cf. Eq. (2.59)). In the late stage with t 1 s, the disagreement between the results for models Bh-R0-Y and Bh-Rw-Y is not negligible. Our interpretation for this is that the formulation and the finite-differencing scheme are not completely the same for the ideal and resistive MHD computations; e.g., in the resistive MHD we solve the equations for the electric field but in the ideal MHD we do not. In particular, for regions in which the electromagnetic energy density is much higher than the rest-mass energy density, we have to introduce an artificial treatment to avoid the crush of the computation, and for such regions the prescription becomes different in the two MHD implementations. Hence, in the long-term simulation, the numerical error accumulated by such an artificial treatment could cause this difference (accordingly, for t 1 s, the numerical results are not very accurate).
Associated with the magnetic braking by the amplified toroidal magnetic fields, the angular velocity profile of the neutron star is always modified for the ideal MHD case (see top two panels of Fig. 2). 1 As we mentioned in Sec. I, initially, the angular velocity is an increase function of in the central region of the neutron star. This profile is modified and approaches gradually a rigidly rotating state for the region in which the strong magnetic field is present. Figure 2 illustrates that the angular velocity in the central region is increased by a factor of ∼ 2 by this effect. Thus, the degree of the differential rotation is reduced by the winding. The timescale for this process is determined by the Alfvén timescale [21,22] where R, B T , and ρ denote the typical radius, toroidal magnetic-field strength, and rest-mass density of the neutron star, and we set h ≈ 1 for simplicity. Unless the magnetic-field profile is significantly modified, after the magnetic braking occurs, the angular velocity should basically show oscillatory behavior with the timescale of τ a , because the angular velocity approximately obeys a hyperbolic equation [21]. However, in the present context, the poloidal magnetic flux is significantly escaped by the neutrino-driven mass outflow and a simple oscillation is not seen. We also note that in the non-axisymmetric case, the turbulence and dynamo effects also would modify the magnetic-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 zdirection with x = 1 km (bottom) for model Rm-B0-Y. The poloidal and toroidal magnetic-field strengths are defined, respectively, by √ B x B x + B z B z γ −1/2 and sign(B y ) B y B y γ −1/2 (≈ B y γ −1/3 ). It is found that the toroidal-field strength (absolute value) is increased far beyond 10 15 G inside the neutron star quite uniformly. By contrast, the maximum value of the poloidal field strength does not increase but rather decreases in the late phase. This decrease is due to the fact that associated with the matter outflow, the magnetic flux is ejected from the neutron star. By this process, instead, a strong poloidal field is produced in the polar direction outside the neutron star (see the bottom panels of Fig. 3): The typical poloidal-field strength at the neutron-star pole is ∼ 10 14 G and ∼ 10 13 G at z ∼ 100 km. Associated with the mass outflow from the neutronstar surface, magnetic-fluxes escape toward a large and large z direction. Because the angular velocity of the escaped matter decreases with the motion towarddirection, the generation of a negative toroidal magnetic field is proceeded along the field lines of the escaped matter. This effect subsequently makes the toroidal magnetic field in the neutron star negative. For model Bm-R0-Y, the mass outflow from the neutron-star surface is appreciably caused by the neutrino-driven wind. As a result, the region with the negative toroidal field expands gradually with time (see the top-right panel of Fig. 3). By contrast, for model Bm-R0-N, the neutrino-driven wind is absent and the mass outflow from the neutron-star surface is weak. Consequently, the region with the negative toroidal field appears only in the vicinity of the neutron star surface (see the second-top-right panel of Fig. 3).
As indicated in Fig. 3, the toroidal magnetic-field energy dominates over the poloidal magnetic-field one after the significant winding: The poloidal-field energy is < 0.1% of the toroidal field energy. However, this may be an artifact in axisymmetric computations in which no dynamo effect is present [26]. In reality, the poloidalfield strength may become much higher than the present results.

Resistive MHD
For the resistive MHD models with low values of the conductivity (i.e., high resistivity), the magnetic-field strength and the maximum value of E B reached by the magnetic winding depend on the values of σ c , and thus, are different from those for the ideal MHD case (see the right panel of Fig. 1 as well as Figs. 2 and 3). Specifically, the maximum magnetic-field strength depends on the ratio, τ dis /τ wind =: R τ : For smaller values of R τ (i.e., R τ 1), the maximum values of the magnetic-field strength and E B are smaller because the resistive dissipation of the electromagnetic energy proceeds faster than the magnetic winding. Thus, for the smaller values of the initial electromagnetic energy, these maximum values are smaller for a given value of σ c . This is found clearly from the results for σ c = 10 7 and 10 8 s −1 (compare the results of models Bh-Rm-Y and Bh-Rh-Y with Bm-Rm-Y and Bm-Rh-Y, respectively). By contrast, for R τ 1, the resistive effect is not very remarkable: For example, the maximum value of E B always exceeds 1% of the rotational kinetic energy, E kin , because the magnetic winding sufficiently proceeds until the saturation of the growth is reached: See the results with σ c ≥ 10 9 s −1 .
We note that irrespective of the maximum value of E B reached, the toroidal magnetic-field energy dominates FIG. 4. Profiles for the density, temperature, specific entropy, electron fraction (left for these four plots), poloidal and toroidal magnetic-field strength, the ratio of the rest-mass to magnetic energy density, and the ratio of the gas pressure to the magnetic pressure (right for these four plots) after the saturation of the growth of the electromagnetic field for models Bm-R0-Y (top panels) and Bm-Rm-Y (bottom panels). over the poloidal magnetic-field one in the end as in the ideal MHD case. Again, this could be an artifact in axisymmetric computations with no dynamo effect in which the poloidal field cannot be amplified [26]. In the nonaxisymmetric case, this is unlikely to be the case, and hence, the poloidal field may be also enhanced in the later stage.
For the case that the maximum value of E B is smaller than ∼ 3% of E kin , the effect of the magnetic braking is weak, and hence, the differential rotation is modified weakly. As found from the middle-right and bottomleft panels of Fig. 2, the angular velocity near the rotational axis (x = 0) does not increase appreciably for σ c ≤ 10 8 s −1 (see also footnote 1). By contrast, the angular-velocity profile for σ c ≥ 10 9 s −1 evolves in a similar manner to that in the ideal MHD case (compare the top-left and bottom-right panels or the top-right and middle-left panels of Fig. 2).
The third row of Fig. 3 shows the evolution of the magnetic-field profile near the equatorial plane for a resistive model (Bm-Rm-Y). Figure 3 illustrates that the magnetic-field strength depends strongly on σ c : For σ c ≤ 10 8 s −1 , the magnetic field is dissipated in 0.1 s and the resulting field strength is by 1-2 orders of magnitude smaller than that in the ideal MHD case (compare the top and third-top panels of Fig. 3). For model Bm-Rm-Y, the toroidal field component varies from the positive to negative values in the neutron star for late times. As we mentioned our interpretation for this by an analysis at the end of Sec. II, this reflects the resistive evolution of the toroidal magnetic field in axial symmetry.
For the case that the resistive dissipation timescale is short, 0.1 s, the thermal energy is generated in the phase in which the neutrino-driven mass outflow is most efficient. However, for such cases, R τ is small, < 1, and hence, the electromagnetic energy caused by the winding is much smaller than the rotational kinetic energy and internal energy of the neutron star. Hence, the generated thermal energy plays only a minor role for the mass outflow and evolution of the neutron star (although this may contribute a bit to enhancing the mass ejection). For R τ 1, the electromagnetic energy reached is ∼ 10% of the rotational kinetic energy. For such cases, however, the dissipation timescale is quite long; in our present setting, it is longer than ∼ 1 s. Since the timescale for launching the neutrino-driven outflow is shorter than or as short as this timescale, the thermal energy generated by the resistive dissipation is not likely to play a major role for the evolution of the remnant neutron stars and mass ejection.
3. Outcome after the winding Figure 4 displays the profiles for the density, temperature, specific entropy, electron fraction (left for these four plots), poloidal and toroidal magnetic-field strength, the ratio of the rest-mass to magnetic energy density, and the ratio of the gas pressure to the magnetic pressure (right for these four plots) after the saturation of the electromagnetic energy is reached for models Bm-R0-Y (ideal MHD model: top panels) and Bm-Rm-Y (moderately high resistivity model: bottom panels). These plots show that the system is composed of a massive neutron star surrounded by a disk and an outflow driven from the polar region of the neutron star toward the z-direction irrespective of the presence or absence of the resistivity.
However, by comparing the profiles for the two models, two quantitative differences are found. First, the magnetic-field strength for Bm-R0-Y is stronger than for Bm-Rm-Y as we already mentioned in Sec. III A 2 (see Fig. 3). This is clearly reflected outside the neutron star in Fig. 4. Second, the property of the outflow driven from the polar surface of the neutron star depends on the strength of an MHD effect: In the presence of the strong MHD pressure (i.e., the ideal MHD case), the matter density in the polar region is smaller, and also, the velocity of the outflow is higher. This modifies the electron fraction of the matter in the polar region and ejecta: The electron fraction is lower for model Bm-R0-Y for which the outflow velocity is higher and the effect of the neutrino irradiation to the ejecta is weaker (see also Sec. III B).

Remark: Comparison with previous work
In the present work, we find that the massive neutron star does not contract significantly and its structure is not modified essentially, after the magnetic winding and braking processes. This conclusion is in contrast to previous work [39][40][41], which showed that differentially rotating neutron stars significantly contract, and if they are hypermassive, they collapse to a black hole, after the magnetic winding effect and associated angular momentum transport by the magnetic braking take place. The major reason for the absence of the significant contraction is that in the present work, a realistic angular velocity profile is determined by the merger simulation (i.e., dΩ/d > 0 in the central region of the neutron star) while in Refs. [39][40][41], they assumed that the angular ve- locity decreases steeply with the cylindrical radius. For dΩ/d > 0, the centrifugal force in the central region of the neutron star increases after the magnetic winding and subsequent magnetic braking proceed, whereas for dΩ/d < 0 (in the previous work), the centrifugal force in the central region decreases by the magnetic-field effects, and hence, the eventual collapse of the hypermassive neutron stars to a black hole can occur. Moreover, the outward angular momentum transport proceeds efficiently for dΩ/d < 0, while for dΩ/d > 0 (the present case) an inward angular momentum transport occurs. Our present work indicates that the choice of the angular velocity profile is an important factor for exploring a realistic evolution process of the hypermassive neutron stars in numerical simulations.

B. Ejecta
As we reported in our previous papers [11,17,25], the mass ejection proceeds from the remnant of binary neutron star mergers through the neutrino irradiation/heating even in the absence of any other effects like the MHD and viscous effects. The rest mass of this neutrino-irradiation component is not very large as 10 −3 M [17,25]. In the context of the present paper, the mass ejection may be enhanced by the MHD effect that results primarily from the amplified toroidal magnetic field. In this subsection, we pay attention to this enhancement. Because we finally find that the mass ejection by the magnetic pressure enhanced by the magnetic winding in the neutron star is a minor effect, we here pay attention to this topic focusing only on the ideal MHD results.
The ejecta component is determined using the same criterion as in Refs. [16,17]; we identify a matter component with |hu t | > h min as the ejecta. Here h min denotes the minimum value of the specific enthalpy in the adopted equation-of-state table, which is ≈ 0.9987c 2 . For the matter escaping from a sphere of r = r ext , we define the ejection rates of the rest mass and energy (kinetic energy plus internal energy) at a given radius and time byṀ whereê := hαu t − P/(ραu t ). The surface integral is performed at r = r ext with dS i = δ ir r 2 ext sin θdθdϕ for the ejecta component. r ext is chosen to be 1000 km in this work.
ρ √ −gu t (= ρ * ) obeys the continuity equation for the rest mass (see Eq. (2.1)), and thus, the time integration for it is a conserved quantity. Also in the absence of gravity and magnetic-field effects, ρê √ −gu t obeys the energy conservation equation, and far from the central region, the sum of its time integration and its gravitational potential energy are approximately conserved (assuming that the electromagnetic energy is much smaller than the kinetic energy). Thus, the total rest mass and energy (excluding the gravitational potential energy and electromagnetic energy) of the ejecta (which escape away from a sphere of r = r ext ) are calculated by where U and T kin are the values of the internal energy and 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 the matter at r = r ext with M being the total gravitational mass of the system [16]. Since the ratio of the internal energy to the kinetic energy of the ejecta decreases with its expansion, we may approximate U/T kin ≈ 0, and hence, E eje,esc by E eje,esc ≈ M eje,esc c 2 + T kin + GM M eje,esc /r ext for the region far from the central object. We then define the average velocity of the ejecta (for the component that escapes from a sphere of r = r ext ) by v eje := 2(E eje,esc − M eje,esc c 2 − GM M eje,esc /r ext ) M eje,esc .
(3.7) In the setting of the present paper, the mass outflow and resulting ejection occur only from the polar region of the neutron star (cf. Fig. 4), because we do not take into account the realistic evolution of the disk surrounding the neutron star, and mass outflow does not occur from it. Thus, in the present work, we perform the surface integral of Eqs. (3.2) and (3.3) only for θ ≤ 15 • where θ denotes the polar angle. Since the initial condition of our simulations is prepared using the result of a merger simulation for binary neutron stars, the dynamical ejecta component is present from the beginning in the computational domain [17]. However, the dynamical ejecta component is located primarily near the equatorial plane. Thus, with the restriction for the surface integral of θ ≤ 15 • , we can focus approximately on the postmerger mass ejection, which comes primarily from the polar region of the neutron star in the present context. Figure 5 shows the rest mass of the post-merger ejecta component as a function of time for four ideal MHD simulations. For comparison, we also plot the result in the absence of the magnetic field (but with the neutrino effects: model B0-Y) [17]. Figure 5 shows that the mass ejection is driven primarily by the neutrino irradiation/heating and pair-annihilation heating toward the polar region, because the rest mass of the ejecta in the absence of these neutrino effects is less than half of that in their presence and also less than that in the absence of the magnetic-field effect (model B0-Y). By comparing the results of models Bh-R0-Y and Rm-R0-Y with that of model B0-Y, it is still found that the magnetic pressure enhanced by the magnetic winding increases the ejecta mass. However, as obviously found from Fig. 5, the rest mass of this ejecta component is of the order of 10 −4 M and by two orders of magnitude smaller than that of the viscosity-driven ejecta, which comes from disks (tori) around the neutron star [11,17]. Therefore, the MHD effect that stems from the winding in the neutron star does not contribute appreciably to increasing the ejecta mass. This is likely to be due to the facts that (i) the MHD effect can strip the material only in a thin polar surface layer of the neutron star for which the total rest mass is tiny and (ii) the gravitational potential near the neutron star is so deep that the mass ejection cannot efficiently occur from its surface.
The present work indicates that the mass ejection by the MHD effect from the neutron star is much less significant than the viscosity-driven mass ejection from disks (tori) surrounding the neutron star. These results suggest that the major source of the mass ejection from the remnant of binary neutron star mergers may be the disk, not the remnant neutron star. However, this speculation should be explored in the future by the simulations in which other important effects such as dynamo effects are taken into account.
A noticeable MHD effect is found in the average velocity and electron fraction of the ejecta (see Fig. 6; cf. Fig. 4 as well). In the absence of the MHD effect (model B0-Y), the average velocity of the ejecta is at most 0.15c. However, in the presence of the strong MHD effect (models Bh-R0-Y and Bm-R0-Y) the average velocity of the ejecta is enhanced to be ∼ 0.3c at t ∼ 1 s. The comparison of the results for models Bh-R0-Y and Bm-R0-Y with Bh-R0-N and Bm-R0-N (models with no irradiation and pair-annihilation of neutrinos) shows that the neutrino effects also enhance the average velocity by the magnitude similar to that by the MHD effect.
Due to the acceleration by the MHD effect, the influence of the irradiation/heating by neutrinos to the electron fraction of the ejecta becomes weak: In the absence of the MHD effect, the electron fraction of the ejecta becomes always high, 0.45, due to the irradiation both by electron neutrinos and anti electron neutrinos, while in the presence of the MHD effect, the electron fraction becomes slightly lower, ∼ 0.4. For example, for the ideal MHD models Rh-B0-Y and Rm-B0-Y, the electron fraction in the polar region is as high as 0.4 before the magnetic field is amplified, but it decreases below 0.4 after the growth of the magnetic field and associated onset of the accelerated mass outflow. This is due to the fact that the neutron-rich nature of the matter outflowed originally from the neutron-star surface is preserved stronger. Thus, the MHD effect in the neutron star does not contribute much to increasing the ejecta mass, but still has the influence to modify the property of the outflow and ejecta. However, this effect is not as strong as the neu-trino irradiation/heating effect (compare models Bh-R0-Y and Bm-R0-Y with Bh-R0-N and Bm-R0-N) as already mentioned.

IV. SUMMARY
We performed ideal and resistive MHD simulations for a remnant neutron star of binary neutron star mergers in general relativity with neutrino effects. As a first step, we paid attention to the effect of the magnetic winding for the evolution of the remnant neutron star and resulting mass ejection. The initial matter profile for the simulations was obtained from a merger simulation. A seed poloidal magnetic field, for which the electromagnetic energy is much smaller than the rotational kinetic and internal energy of the system, was initially superimposed inside the remnant neutron star, and we focused only on the evolution of it; we did not pay attention to the MHD evolution of the disk (torus) surrounding the neutron star.
Because of the magnetic winding effect, the toroidal magnetic field is generated and amplified in this setting. Since the toroidal magnetic-field strength increases linearly with time, the electromagnetic-field energy increases as ∝ t 2 in the early growth stage. When the electromagnetic-field energy exceeds ∼ 3% of the rotational kinetic energy, the magnetic braking plays an important role for the redistribution of the angular momentum inside the neutron star. For the ideal MHD case, this always occurs and, as a result of the angular momentum redistribution, the angular velocity, which is initially an increase function of the cylindrical radius, is enforced to approach a rigidly rotating state for the region in which the magnetic braking works well. The maximum electromagnetic energy reached is ∼ 10% of the rotational kinetic energy of the neutron star, i.e., the maximum value of E B is ∼ 10 52 erg. By the angular-momentum redistribution, the rotational kinetic energy of the neutron star is not significantly changed, while at the late stage, the electromagnetic energy relaxes to ∼ 10 51 erg after the magnetic braking is exerted, because the magnetic flux escapes from the neutron star associated with the mass outflow.
For the resistive MHD simulation, the maximum electromagnetic energy reached by the winding effect depends strongly on the magnitude of the conductivity; specifically, the maximum value is determined by the ratio, R τ = τ dis /τ wind . For R τ 1, the magnetic winding proceeds until the electromagnetic energy reaches ∼ 10% of the rotational kinetic energy as in the ideal MHD case. On the other hand, for R τ < 1, the resistive dissipation plays a role for suppressing the growth of the electromagnetic energy. If the electromagnetic energy does not reach ∼ 3% of the rotational kinetic energy, the magnetic braking effect is weak, and hence, the differentially rotating state of the neutron star is preserved; the increase of the angular velocity near the central region is suppressed.
It is also found that the magnetic field amplified by the winding in the neutron star does not contribute much to enhancing the mass ejection; for the mass ejection from the neutron star, the neutrino effects such as neutrino irradiation/heating and pair-annihilation heating are more important than the MHD effect. Our interpretation for this result is that (i) the MHD effect in the neutron star can strip the material only in the thin polar surface region of the neutron star for which the total rest mass is tiny and (ii) the gravitational potential near the neutron star is so deep that the mass outflow cannot be efficient. The present result suggests that the major source of the mass ejection from the remnant of binary neutron star mergers may not be remnant neutron star but the disk surrounding it. However, this speculation should be examined in more detail by the simulations in which other important effects such as dynamo effects are taken into account.
In the present work, we performed axisymmetric simulations. This implies that we might overlook important MHD effects. One important missing effect is the turbulence induced by Parker [42] and Taylor instabilities [43] (e.g., Ref. [44]). As we showed in this paper, the toroidal magnetic field is inevitably amplified in the presence of a seed poloidal magnetic field. It is well-known that in the presence of the strong toroidal magnetic fields, (nonaxisymmetric) Parker and Taylor instabilities can induce a turbulence, which could further induce the magneticfield amplification through the dynamo effect. It is not clear at all what happens in such a situation. In addition, it is not very clear what the final configuration of the magnetic field profile is after the onset of the turbulence (e.g., Ref. [45]). To explore these questions, we obviously need three-dimensional MHD simulations in the subsequent work. An alternatively phenomenological approach is to incorporate a dynamo term to induce the turbulence state [35,36]. By this prescription, we do not need three-dimensional simulations: An axisymmetric simulation would be reasonable to obtain some insight on the effect of the turbulence. We plan to perform simulations with a dynamo term in the subsequent work.

ACKNOWLEDGMENTS
We thank Kenta Kiuchi and Federico Carrasco for helpful discussions. MS and SF thank Yukawa Institute for Theoretical Physics, Kyoto University for their hospitality during the first corona pandemic time in Germany, in which this project was started. This work was in part supported by Grant-in-Aid for Scientific Research (Grant Nos. JP16H02183 and JP20H00158) of Japanese MEXT/JSPS. Numerical computations were performed on Sakura and Cobra clusters at Max Planck Computing and Data Facility.

Appendix A: Test-bed simulations for resistive MHD implementation
We here present the results for a suit of test-bed simulations performed with our resistive MHD implementation. We employed several test-bed problems introduced in Ref. [36] (see also Refs. [46][47][48][49]): the problems of selfsimilar current sheet, resistive shock tube, and resistive rotor. We also performed test simulations for the propagation of the electromagnetic wave packet in the flat spacetime and for a dynamo closure.

Self-similar current sheet
This test-bed problem was first proposed in Ref. [46] and subsequently employed in many references [36,[47][48][49]. This is the test to check whether the magnetic field correctly diffuses out with the special-relativistic resistive MHD implementation in the limit that the gas pressure is much higher than the magnetic pressure. In the one dimensional problem in which all the quantities depend only on the x coordinate, an approximate (nearly exact) solution is written, e.g., in the following form: where erf denotes the error function and C is a constant with C 2 /8π much smaller than the gas pressure. Here we set C = 1 as the choice of the units for simplicity. Note that our notation for the basic MHD equations is different from the previous papers because of the presence of the factor 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 = (Γ − 1)ρε with Γ = 4/3 is employed. σ c is set to  Figure 7 plots B y at t = 1 and 10. For t = 10, we also plot the functional form of B y described in Eq. (A1). Figure 7 illustrates that our implementation can reproduce the solution well. In this setting the maximum relative error is ≈ 2.5 × 10 −4 .

Resistive shock tube
This problem was first proposed in Ref. [47]. The initial condition for this problem proposed in Ref. [36] is 4π (this renormalization is necessary to align the unit with the previous work) and v i = u i /u t . The initial electric field is given by E i = − ijk v j B k . The Γ-law equation of state with Γ = 5/3 is employed. σ c is chosen to be 10 −3 , 1, 10, 10 2 , and 10 10 . The computation is performed from t = 0 to 0.55 covering the computational domain of x = [−0.5 : 0.5] by 401 uniform grid points. The results are shown for ρ and B y in Fig. 8 as in Ref. [36]. A small unphysical oscillation could be seen for the curves of ρ and B y in the absence of the Kreiss-Oliger dissipation for evolving Maxwell's equations. However, this dissipation term with a small coefficient of 1/640 cures the unphysical oscillation. We find that the numerical solutions are quite similar to those found in Ref. [36].

Resistive rotor
This is a two-dimensional problem with the coordinates (x, z) and the initial condition proposed in Ref. [36] is given in the following manner. Inside a radius of r = √ x 2 + z 2 ≤ 0.1, a uniformly rotating medium with ρ = 10 and the uniform angular velocity of Ω = 8.5 is prepared (i.e., v x = −zΩ and v z = xΩ). Outside the circle of r = 0.1, on the other hand, we set ρ = 1 and Ω = 0. For the entire region, the pressure and the magnetic field are initially uniform as p = 1 and (B x ,B y ,B z ) = (1, 0, 0). The initial electric field is again determined by E i = − ijk v j B k . The Γ-law equation of state with Γ = 4/3 is employed, and the system is evolved from t = 0 to 0.3.
The computational domain is chosen to be x = [−0.5 : 0.5] and z = [−0.5 : 0.5] and is covered by a uniform grid of 401 × 401 points. This problem is numerically solved for σ c = 10 10 (nearly ideal MHD case) and 1. The results for the profiles of the pressure and z-component of the electric field are shown in Fig. 9. By comparing these results with Fig. 3 of Ref. [36], we find a good (qualitative/semi-quantitative) agreement.

Evolving electromagnetic fields
To check that our implementation for evolving the electromagnetic fields works well in a practical computational domain, we also performed several test-bed simulations in three spatial dimensions varying the value of σ c from 0 to high values. For σ c = 0 and σ c → ∞ in the flat spacetime of u i = 0, it is easy to derive analytic solutions for the electromagnetic equations, and thus, we compare the results of the simulations with these analytic solutions.
For σ c = 0, the electromagnetic fields obey wave equations in vacuum, and general solutions for the basic equations are easily derived. For example, an axisymmetric general solution for the dipole radiation is written as and B ϕ = 0. For this case, E ϕ is only non-zero component of E i and it is derived straightforwardly from the field equation. Here, f (u) is an arbitrary regular function, and for f (u) = −ue −u 2 /2 , the initial condition becomes and E i = 0. Note here that we choose the width of the wave packet (i.e., unity) as the unit of the length in this problem.   11. The maximum value of B y for the steady dynamo test for σc = 10 and α d = 0.2 with λx = 2, 1, 1/2, 1/3, 2/7, 1/4, 2/9, a = 0.1, and b = 0 (no perturbation in the z direction) and with λx = 1/4, λz = 1/2, a = 0.1, and b = 0.01. Note that the fastest growing mode has the wavelength of 1/2 while for the marginally stable mode it is 1/4. λ in the inset of the plot implies λx.
The relative error of the numerical solution to the analytic one, e.g., at t = 5, is of order 10 −3 with this setting (except for the outer boundaries). This illustrates that our implementation can well follow the propagation of the electromagnetic waves.
For large values of σ c , on the other hand, the field equations relax to parabolic ones. For the parabolic equations with the initial condition of Eq. (A3), a solution for the dipole magnetic field is written as We again evolve this wave packet in the domain of  Fig. 10 plots the profiles of B z along the x axis at t = 0, 50, and 100 (solid curves) together with the analytic solutions (dashed curves). With the chosen grid resolution, the maximum relative error size is ≈ 4% at t = 100 (except for the region near the outer boundaries). The largest error is always located near the region at which B z = 0, and for other regions, the accuracy is much better.

Steady dynamo
Following Ref. [36], we also performed this simple test to check that our implementation works well in the presence of a dynamo term. In this test-bed problem, again, we pay attention only to solving the electromagnetic field setting 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 are the angular frequency and the wave number. Then, the following dispersion relation is derived: For this equation, we find that the exponential growth mode is present for the case that k < 4πσ c α d , and the resultant expression of ω is and for the fastest growing mode with k = 2πσ c α d , this becomes Because the basic equation is linear in E i and B i , it is straightforward to extend this analysis in the multidimensional case as long as we focus only on the transverse component. That is, for the case that we initially 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 is determined simply by analyzing the dispersion relations in the x and z direction independently. Taking into account this fact, we perform two-dimensional simulations setting the region of x = [−1 : 1] and z = [−1 : 1] with the periodic boundary conditions in both coordinates.
We then prepare the following initial condition: and E i is determined by the ideal MHD condition. Here, a and b are constants. In this setting, for the case that 2λ i σ c α d > 1, the initial seed field should grow exponentially with time, where λ i := 2π/k i (i is x or z) denotes the wave length. Note that our setting of the computational domain can follow the waves only with λ i ≤ 2. Thus, the simulations are performed for σ c = 10 and α d = 0.2. For this setting, the fastest growing mode has λ i = 1/2 while the marginally stable mode has λ i = 1/4 (i.e., the mode with λ i > 1/4 is unstable). Figure 11 plots the evolution of the maximum value of B y for σ c = 10 and α d = 0.2 with λ x = 2, 1, 1/2, 1/3, 2/7, 1/4, 2/9, a = 0.1, and b = 0 (no perturbation in the z direction) and with λ x = 1/4, λ z = 1/2, a = 0.1, and b = 0.01. The growth rates of the unstable modes are well captured in each simulation. This figure illustrates that our implementation can derive the predicted results irrespective of the prepared initial conditions. Evolution of the magnetic-field energy for the remnant massive neutron star initially with a purely toroidal magnetic field for σc = 10 7 , 10 8 , 10 9 , and 10 11 s −1 (solid curves). The dashed lines denote ∝ exp(−t/τ diss ) with τ diss = 0.01, 0.1, and 1.0 s for σc = 10 7 , 10 8 , and 10 9 s −1 , which are consistent with Eq. (2.59). The result by the ideal MHD implementation is also shown by the dotted curve, which agrees approximately with the curve for σc = 10 11 s −1 .
Appendix B: Resistive evolution of the massive neutron star with purely toroidal magnetic field In this section, we present the results for the resistive MHD evolution of the merger remnant neutron star with purely toroidal magnetic fields and demonstrate that our implementation can follow the resistive dissipation of the toroidal magnetic field successfully. Figure 12 plots the evolution of the electromagnetic energy for the remnant massive neutron star, for which a purely toroidal magnetic field of Eq. (2.54) is superimposed at t = 0, in the axisymmetric resistive MHD simulations with σ c = 10 7 , 10 8 , 10 9 , and 10 11 s −1 (solid curves). The result of an ideal MHD simulation is shown together (dashed curve). Since the toroidal magnetic field is simply superimposed and thus the initial condition is not exactly in an equilibrium state, the electromagnetic energy initially decreases by 10-20% even in the absence of the resistivity during the early relaxation stage. The subsequent long-term gradual decrease of the electromagnetic energy for the ideal MHD case would be partly due to the numerical dissipation or diffusion. 2 Note again that this evolution process is valid only in the assumption of the axial symmetry. In general cases, non-axisymmetric instability such as Parker and Taylor instability [42,43] could occur and the evolution process could be significantly modified. For the resistive MHD case, the magnetic field decreases exponentially with time after the early relaxation stage. The short-dashed slopes denote the relation of ∝ exp(−t/τ diss ) where τ diss is the decay timescale, and are plotted for approximately fitting the curves for the resistive MHD results. It is found that the dissipation timescale of the magnetic field is ≈ 0.01, 0.1, and 1.0 s for σ c = 10 7 , 10 8 , and 10 9 s −1 , which are consistent with the timescales estimated by Eq. (2.59). For σ c = 10 11 s −1 , the dissipation timescale would be ∼ 10 2 s, and much longer than the simulation time. Therefore, the result for this case agrees approximately with that in the ideal MHD simulation: Our ideal and resistive MHD implementations can derive approximately the identical results.

Appendix C: Grid resolution dependence
As a convergence test, we performed numerical simulations 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 as Fig. 1 but for the comparison in the evolution of the electromagnetic energy for two different resolutions. The solid and dashed curves show the results for the cases that 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 magneticfield strength by the winding effect is followed with the lower grid resolution fairly well. However, the amplification of the magnetic-field strength is suppressed for the lower resolution, and thus, the maximum magneticfield energy becomes smaller. This effect is not appreciable for the case that the effects of the neutrino irradiation/heating and pair-annihilation are switched off.
However, in the presence of these neutrino effects, the suppression is quite large: The maximum magnetic-field strength for the lower resolution is by a factor of ∼ 2 smaller than for the higher one. It appears that the mass ejection associated with these neutrino effects enhance the magnetic-field dissipation or magnetic flux outflow from the neutron star. Thus the poor convergence is likely to be associated with a poor convergence in the neutrino transfer (cf. Ref. [25]).
For the resistive MHD simulations, we also find the tendency similar to that found for the ideal MHD case with neutrino effects. That is, the dissipation or magnetic flux escape from the neutron star is spuriously enhanced for the lower grid resolution. However, overall, the evolution process of the neutron star is not modified qualitatively by the grid resolution.