General-Relativistic Gauge-Invariant Magnetic Helicity Transport:
Basic Formulation and Application to Neutron Star Mergers

Jiaxi Wu TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA Elias R. Most TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA

Abstract

Dynamo processes are ubiquitous in astrophysical systems. In relativistic astrophysical systems, such as accretion disks around black holes or neutron stars, they may critically affect the launching of winds and jets that can power electromagnetic emission. Dynamo processes are governed by several microscopic parameters, one of them being magnetic helicity. As a conserved quantity in nonresistive plasmas, magnetic helicity is transported across the system. One important implication of helicity conservation is, that in the absence of helicity fluxes some mean-field dynamos can be quenched, potentially affecting the large-scale evolution of the magnetic field. One of the major challenges in computing magnetic helicity is the need to fix a meaningful electromagnetic gauge. We here present a fully covariant formulation of magnetic helicity transport in general-relativistic plasmas based on the concept of relative helicity by Berger & Field and Finn & Antonsen. This formulation is separately invariant under gauge-transformation of the Maxwell and Einstein equations. As an application of this new formalism we present the first analysis of magnetic helicity transport in the merger of two neutron stars. We demonstrate the presence of global helicity fluxes into the outer layers of the stellar merger remnant, which may impact subsequent large-scale dynamo amplification in these regions.

I Introduction

Astrophysical dynamos are important for a variety of systems, including stars [1], accretion disks [2], planets [3] and galaxies [4]. Dynamos can operate on different scales, causing turbulent magnetic field amplification on small (microscopic scales) to mean-field processes in large scale shear flows [5, 6]. In relativistic systems, turbulent dynamo amplification may aid the launching of jets and outflows from neutron stars [7, 8, 9, 10, 11, 12] and black hole accretion disks [13, 14, 15].

While the details of these individual processes are complicated and will depend on the specific dynamo mechanism operating, several of them including the $\alpha\Omega$ -dynamo [16] potentially active in neutron star mergers [8], are affected by a magnetic field invariant, magnetic helicity.

Magnetic helicity is a topological invariant of the plasma [17]. In simple terms, it quantifies the amount of magnetic field line linkage in a given volume. A high degree of helicity implies very intricately intertwined field geometries. While helicity can be dissipated by resistivity, in nonresistive plasmas (such as in many relativistic systems) it is a globally conserved quantity. In other words, any magnetic field amplification or rearrangement will be subject to the constraints imposed by the initial helicity content in a given volume. Prominent examples of this concern the relaxation of a given field into its lowest energy configuration [18], which can affect large scale structures of the field. In addition, Refs. [16, 19] (see also Ref. [20]) have argued that some dynamo actions may be quenched in the presence of finite helicity, such that the field may not be able to rearrange itself into large scale structures. This effect of $\alpha-$ quenching can particularly affect $\alpha\Omega$ -dynamos, though other dynamo models [21] potentially active in accretion disks [15] are intrinsically non-helical.

A potential way to circumvent the limitations of a fixed helicity inside a given volume are helicity fluxes, which allow helicity to be transported in and out of this volume. Since many of the (relativistic) systems, we are interested in, feature complicated large scale flow structures, it is conceivable that such a situation may not be uncommon. With a wealth of numerical simulations available nowadays, this question can in principle be answered.

In practice, a major obstacle in computing helicity fluxes is the gauge dependence of any local helicity measurement [22]. Put differently, while the global value of the helicity over a volume enclosing all magnetic field lines of the system is well defined and gauge invariant, the value of its contributions throughout the volume is gauge dependent.

Addressing this problem is not straightforward. One possibility is to relax the helicity definition to mean the linkage with respect to a single field line, which can be more easily computed [23]. Previous numerical studies have commonly considered a fixed gauge, which may be appropriate given special geometries of the problem [24]. For a general flow, as we plan to consider, such an approach may not always be feasible. Apart from uncertainties in quantifying transport, helicity conservation may also be affected by artificial dissipation of the numerical scheme [25].

Leveraging advances in solar dynamo physics [26], in this work we formulate a general-relativistic version of gauge-invariant magnetic helicity transport in arbitrary spacetimes. This formulation relies on a clean separation of field lines enclosed in a subvolume from those leaving this volume in a formulation first developed by Berger & Field [22] and Finn & Antonsen [27].

Our paper is structured as follows: In Sec. II, we present a summary of the general idea of helicity and its transport, which we then mathematically formulate in Sec. III. In Sec. IV, we demonstrate how this formulation can be applied to a neutron star merger simulation. Throughout this work, we use geometric units, $G=c=1$ , and Heaviside-Lorentz units for the electromagnetic fields. We further use Greek indices $(\mu,\nu,\alpha,\dots)$ to indicate indices running from $0$ to $3$ , and Latin indices $(i,j,k,\dots)$ to represent spatial indices from $1$ to $3$ .

II Basic picture

Magnetic helicity is a measurement of magnetic flux tube linkage, including terms of “self-helicity” and “mutual helicity”. The self-helicity quantifies how magnetic field lines are twisted and writhed with respect to itself, while the mutual helicity accounts for crossings of different field lines [28]. Magnetic helicity can be defined in terms of a helicity density,

\displaystyle\mathcal{H}=\boldsymbol{A}\cdot\boldsymbol{B}\,,

(1)

where $\boldsymbol{A}$ is the magnetic vector potential, and $\boldsymbol{B}=\boldsymbol{\nabla}\times\boldsymbol{A}$ is the magnetic field. While $\mathcal{H}$ can be computed everywhere in the domain, due to the gauge degree of freedom of the magnetic vector potential, $\boldsymbol{A}\rightarrow\boldsymbol{A}+\boldsymbol{\nabla}\Lambda$ , the local helicity density $\mathcal{H}$ is fundamentally not gauge invariant.

At the same time, the global helicity,

H=\int_{V}{\text{d}^{3}x\,\mathcal{H}}=\int_{V}{\text{d}^{3}x\,\boldsymbol{A}\cdot\boldsymbol{B}}\,,

(2)

is gauge invariant, provided that the volume $V$ encloses all magnetic field lines, i.e.,

H^{\prime}=H+\int_{V}{\text{d}^{3}x\,\boldsymbol{\nabla}\Lambda\cdot\boldsymbol{B}}=H-\oint_{\partial V}{\text{d}S\,{\boldsymbol{s}}\cdot\Lambda\boldsymbol{B}},

(3)

where the last term vanishes if $B$ is fully enclosed in the volume $V$ . Here ${\boldsymbol{s}}$ is the unit normal pointing inward to the boundary $\partial V$ .

However, in realistic situations one might be interested in a more local measure, i.e., the amount of helicity in the vicinity of a black hole rather than of the entire accretion disk. The helicity concept above can be generalized to the case, where one separates magnetic field lines entirely contained inside a volume $V$ from those leaving it through the boundary. Since from a measurement in $V$ alone one cannot say anything about these “open” field lines, asking about linkage numbers for magnetic field lines entirely contained in a volume is meaningful within the helicity framework introduced above.

This concept of relative helicity goes back to Berger & Field [22] and Finn & Antonsen [27]. Schematically, the situation we are after is shown in Fig. 1. Here we can see the background field $\boldsymbol{B}_{R}$ , which consists solely of field lines leaving the volume $V$ in question, separated from those entirely contained in it $\boldsymbol{B}-\boldsymbol{B}_{R}$ . Following Ref. [22], the background field can be constructed as a potential flow, satisfying

	$\displaystyle\boldsymbol{B}_{R}$	$\displaystyle=-\boldsymbol{\nabla}\psi_{R}\,,$		(4)
	$\displaystyle\boldsymbol{\nabla}^{2}\psi_{R}$	$\displaystyle=0\,,$		(5)

such that,

\left.{\boldsymbol{s}}\cdot\boldsymbol{B}\right|_{\partial V}=\left.{\boldsymbol{s}}\cdot\boldsymbol{B}_{R}\right|_{\partial V}=-\left.\frac{\partial\psi_{R}}{\partial s}\right|_{\partial V}\,,

(6)

which ensures that the reference fields (red lines in Fig. 1) match the outgoing magnetic field lines on the boundary $\partial V$ of the volume $V$ . Given the background magnetic field, one can uniquely define an associated magnetic vector potential $\boldsymbol{A}_{R}$ , such that $\boldsymbol{B}_{R}=\boldsymbol{\nabla}\times\boldsymbol{A}_{R}$ and $\boldsymbol{\nabla}\cdot\boldsymbol{A}_{R}=0$ . Equipped with this split into background and foreground fields, the relative helicity density is defined as [22, 27],

h=(\boldsymbol{A}+\boldsymbol{A}_{R})\cdot(\boldsymbol{B}-\boldsymbol{B}_{R})\,.

(7)

Refer to caption — Figure 1: A sketch showing the concept of reference fields. Magnetic fields shown as orange field lines are completely confined within the volume boundary, and thus can be used to calculate gauge-independent relative helicity. In contrast, red field lines are crossing the boundary and can be subtracted out by introducing the reference field $\boldsymbol{B}_{R}$ corresponding to these field lines.

Following Refs. [22, 27], we can express the relative helicity as,

H=\int_{V}{\text{d}^{3}x\,h}=\int_{V}{\text{d}^{3}x\,(\boldsymbol{A}+\boldsymbol{A}_{R})\cdot(\boldsymbol{B}-\boldsymbol{B}_{R})}.

(8)

It can be shown that Eq. (8) defined in this way is gauge-invariant and retains its usual meaning [22, 27].

Computing this term in practice is rather involved. First one needs to compute the background field, which will be the solution to a three-dimensional elliptic equation with the boundary conditions given by Eq. (6). This is rather unwieldy in practice, especially if one wanted to do it for many subdomains of a numerical simulation. Relative magnetic helicity is also not additive, meaning that it has to be recomputed separately for every volume $V$ under considerations [29].

Computing the full relative helicity integral (8) is, however, not the only way to track helicity. It can be shown that the original helicity expression gives rise to a conservation law,

\displaystyle\frac{\partial\mathcal{H}}{\partial t}+\boldsymbol{\nabla}\cdot\left(\boldsymbol{E}\times\boldsymbol{A}-\Phi\boldsymbol{B}\right)=-2\boldsymbol{E}\cdot\boldsymbol{B}\,,

(9)

where $\Phi$ is the scalar potential, and the RHS vanishes for a nonresistive plasma, for which $\boldsymbol{E}=-\boldsymbol{v}\times\boldsymbol{B}$ . While the above transport equation is not locally gauge invariant, we can take a similar approach and re-formulate the relative helicity as a transport problem (relative to a given volume $V$ ) [26],

	$\displaystyle\frac{\partial H}{\partial t}$	$\displaystyle=$	$\displaystyle-2\int_{{V}}{\text{d}^{3}x\,\boldsymbol{E}\cdot\boldsymbol{B}}$		(10)
			$\displaystyle-\oint_{\partial V}{\text{d}S\,{\boldsymbol{s}}\cdot[(\boldsymbol{A}+\boldsymbol{A}_{R})\times(\boldsymbol{E}-\boldsymbol{E}_{R})]}.$		(10)

Due to the inclusion of the background field, this expression fundamentally maintains gauge invariance (see Appendix B).

While the above expression does at first glance seemingly not simplify our problem, it turns out that after substantial algebraic manipulation in the case of an ideal (nonresistive) plasma one can recast this expression as [26],

\frac{\partial H_{S}}{\partial t}=2\oint_{\partial V}{\text{d}S\,\zeta B_{n}},

(11)

where $\zeta$ is the lamellar part (i.e. the part with zero normal component of the curl of the surface gradient) of the electric field, which can be calculated by solving a surface elliptic partial differential equation on the boundary $\partial V$ of the enclosing volume,

\nabla_{S}^{2}\zeta=-\boldsymbol{\nabla}_{S}\cdot\boldsymbol{E}_{S}.

(12)

We denote the induced two-dimensional operators and fields on the boundary with a subscript $S$ . In particular, for a spherical volume $\zeta$ is uniquely determined by the solution of Eq. (12). Remarkably, this equation only depends on the ideal electric field, not on the vector potential, $\boldsymbol{A}$ , making it suitable also for most numerical codes that only evolve $\boldsymbol{B}$ [30]. It should not come as a surprise that this expression must involve a surface elliptic equation, since the reference magnetic field approach in essence counts the number of open magnetic field lines, which is a global problem on the boundary.

In the remainder of this work, we will generalize this approach [26] to general-relativity in arbitrary spacetimes. We will show that when carefully defining the relevant electric and magnetic fields within a 3+1 decomposition of spacetime, we can retain the the simplicity character of Eqs. (11) and (12).

III General-relativistic helicity transport

In this section, we want to derive a general-relativistic expression for gauge-invariant magnetic helicity transport. Before we do so, we will briefly review the covariant form of the Maxwell equations, as well as decomposition of spacetime useful for our purposes.

III.1 General-relativistic electrodynamics

Covariant electrodynamics can be formulated in terms of a covariant vector potential $\mathcal{A}_{\mu}$ , which gives rise to a field strength tensor,

\displaystyle F^{\mu\nu}

\displaystyle=\nabla^{\mu}\mathcal{A}^{\nu}-\nabla^{\nu}\mathcal{A}^{\mu}\,,

(13)

where $\nabla_{\mu}$ denotes the covariant derivative relative to the four-dimensional spacetime metric $g_{\mu\nu}$ . It will turn out to be convenient to further introduce the dual of the field strength tensor,

\displaystyle{~\!{}^{\ast}\!F}^{\alpha\beta}=\frac{1}{2}\varepsilon^{\alpha\beta\mu\nu}F_{\mu\nu},

(14)

where $\varepsilon^{\alpha\beta\mu\nu}=-\frac{1}{\sqrt{-g}}\epsilon^{\alpha\beta\mu\nu}$ is the four-dimensional Levi-Civita tensor, and $\epsilon^{\alpha\beta\mu\nu}$ is the permutation symbol. The corresponding covariant Levi-Civita tensor is $\varepsilon_{\alpha\beta\mu\nu}=\sqrt{-g}\epsilon_{\alpha\beta\mu\nu}$ .

The evolution equations for the field strength tensor are given by the covariant form of the Maxwell equations,

	$\displaystyle\nabla_{\beta}F^{\alpha\beta}$	$\displaystyle=\mathcal{J}^{\alpha}\,,$		(16)
	$\displaystyle\nabla_{\beta}\prescript{\ast}{}{F}^{\alpha\beta}$	$\displaystyle=0\,,$		(17)

where $\mathcal{J}^{\alpha}$ is the electric current.

We can now define a helicity current [31],

\displaystyle\mathcal{H}^{\mu}=\mathcal{A}_{\nu}{~\!{}^{\ast}\!F}^{\mu\nu}\,.

(18)

Using the Maxwell equations (16), it follows that,

\displaystyle\nabla_{\mu}\mathcal{H}^{\mu}={~\!{}^{\ast}\!F}^{\alpha\beta}F_{\alpha\beta}\,,

(19)

which is the covariant generalization of the helicity transport equation (9).

The energy evolution of the electromagnetic field can be associated with an energy momentum tensor,

\displaystyle T_{\rm EM}^{\mu\nu}=F^{\mu\alpha}F^{\nu}_{\alpha}-g^{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}\,,

(20)

which satisfies [32],

\displaystyle\nabla_{\mu}T^{\nu\mu}_{\rm EM}=-F^{\mu\nu}\mathcal{J}_{\mu}\,.

(21)

III.1.1 3+1 decomposition of spacetime

In order to make a connection to the formulation in terms of electric and magnetic fields discussed in Sec. II, we need to introduce a normal observer, $n_{\mu}=(-\alpha,0,0,0)$ , whose worldline trajectory will provide a time direction. Here $\alpha$ is the lapse function. We do so using the conventional 3+1 decomposition of spacetime [33], where the metric is expressed as,

g_{\mu\nu}\text{d}x^{\mu}\text{d}x^{\nu}=-\left(\alpha^{2}-\beta_{i}\beta^{i}\right)\text{d}t^{2}+2\beta_{i}\text{d}x^{i}\text{d}t+\gamma_{ij}\text{d}x^{i}\text{d}x^{j}\,,

(22)

with spatial metric $\gamma_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$ , $n^{\nu}=(1/\alpha,-\beta^{i}/\alpha)$ , and $\beta^{i}$ being the coordinate shift.

Within this coordinate choice, we can split the covariant vector potential $\mathcal{A}_{\mu}$ and electric current $\mathcal{J}_{\mu}$ into,

	$\displaystyle\mathcal{A}_{\mu}$	$\displaystyle=\Phi n_{\mu}+A_{\mu}\,,$		(23)
	$\displaystyle\mathcal{J}_{\mu}$	$\displaystyle=qn_{\mu}+J_{\mu}\,,$		(24)

where $\Phi$ is the scalar potential, $A_{\mu}$ the vector potential, $q$ the electric charge density, and $J_{\mu}$ the electric current on the three dimensional hypersurface induced by $n_{\mu}$ .

Using this split, we can introduce electric, $E^{\mu}$ , and magnetic, $B^{\mu}$ , fields via [32],

	$\displaystyle E^{\mu}$	$\displaystyle=n_{\nu}F^{\mu\nu}=\alpha F^{0\mu}\,,$		(25)
	$\displaystyle B^{\mu}$	$\displaystyle=n_{\nu}{~\!{}^{\ast}\!F}^{\mu\nu}=\alpha{~\!{}^{\ast}\!F}^{0\mu}\,.$		(26)

One can show that,

\displaystyle B^{i}=\varepsilon^{ijk}\partial_{j}A_{k}\,,

(27)

where $\varepsilon^{\mu\nu\kappa}=n_{\alpha}\varepsilon^{\mu\nu\kappa\alpha}$ is the three-dimensional Levi-Civita tensor. This establishes consistency between the 3+1 formulation of electrodynamics and those in flat spacetimes.

We can similarly decompose the field strength tensor into [32],

	$\displaystyle F^{\alpha\beta}$	$\displaystyle=n^{\alpha}E^{\beta}-n^{\beta}E^{\alpha}-\varepsilon^{\alpha\beta\mu\nu}B_{\mu}n_{\nu}\,,$		(28)
	$\displaystyle\prescript{\ast}{}{F}^{\alpha\beta}$	$\displaystyle=-n^{\alpha}B^{\beta}+B^{\alpha}n^{\beta}-\varepsilon^{\alpha\beta\mu\nu}E_{\mu}n_{\nu}\,.$		(29)

Within the 3+1 split, the Maxwell equations take the familiar form,

	$\displaystyle\mathcal{D}_{i}E^{i}=q,$		(30)
	$\displaystyle\mathcal{D}_{i}B^{i}=0,$		(31)
	$\displaystyle\partial_{t}(\sqrt{\gamma}E^{i})-\epsilon^{ijk}\partial_{j}\tilde{B}_{k}=-\sqrt{\gamma}J^{i},$		(32)
	$\displaystyle\partial_{t}(\sqrt{\gamma}B^{i})+\epsilon^{ijk}\partial_{j}D_{k}=0.$		(33)

Here $\mathcal{D}_{i}x^{i}=\frac{1}{\sqrt{\gamma}}\partial_{i}(\sqrt{\gamma}x^{i})$ is the three-dimensional covariant derivative with respect to $\gamma_{ij}$ . Finally, we have introduced effective expressions for the electric and magnetic fields that enter the flux terms in the above expressions,

	$\displaystyle D^{i}=\alpha E^{i}+\varepsilon^{ijk}\beta_{j}B_{k},$		(34)
	$\displaystyle\tilde{B}^{i}=\alpha B^{i}-\varepsilon^{ijk}\beta_{j}E_{k}.$		(35)

When using such definitions, it can be shown that the Maxwell equations take a form closest to their flat spacetime counterparts [34].

III.1.2 Surface Helicity Transport

With the formalism for covariant electrodynamics in place, we are now in a position to formulate relative helicity transport within the 3+1 language. Following our discussion in Sec. II, we need to define a reference magnetic field, $B^{i}_{R}=-\mathcal{D}^{i}\psi_{R}$ , for a suitable scalar potential $\psi_{R}$ , where

	$\displaystyle-\mathcal{D}_{i}\mathcal{D}^{i}\psi_{R}=\mathcal{D}_{i}B_{R}^{i}=0,$		$\displaystyle s_{i}B^{i}\|_{S}=s_{i}B_{R}^{i}\|_{\partial V},$		(36)
	$\displaystyle\mathcal{D}_{i}A_{R}^{i}=0,$		$\displaystyle s_{i}A_{R}^{i}\|_{\partial V}=0\,,$		(37)

with a corresponding vector potential, $\mathcal{A}^{\mu}_{R}=A^{\mu}_{R}$ , such that $B^{i}_{R}=\varepsilon^{ijk}\partial_{j}\left(A_{R}\right)_{k}$ , and vanishing scalar potential $\Phi_{R}=0$ . The presence of relativity does almost not affect these equations, except through the necessity of using three-dimensional covariant derivatives, $\mathcal{D}_{i}$ in order to preserve the constraint condition (31).

In relativity, we cannot express the helicity density in terms of the magnetic field, but need the full field strength tensor. We therefore define the corresponding field strength tensor of the background field as

\displaystyle F^{\mu\nu}_{R}=\nabla^{\mu}\mathcal{A}^{\nu}_{R}-\nabla^{\nu}\mathcal{A}^{\mu}_{R}\,.

(38)

Using these definitions, we propose to generalize the helicity density (8) to a relativistic relative helicity current, $h^{\mu}$ in the following way

h^{\mu}=(\mathcal{A}_{\nu}+\mathcal{A}_{R\nu})(\prescript{\ast}{}{F}^{\mu\nu}-\prescript{\ast}{}{F_{R}}^{\mu\nu})\,.

(39)

For easier comparison with their non-relativistic expressions, we now expand this expression into the normal electric and magnetic fields,

	$\displaystyle h^{0}$	$\displaystyle=$	$\displaystyle\frac{1}{\alpha}(A_{j}+A_{Rj})(B^{j}-B_{R}^{j}),$		(40)
	$\displaystyle h^{i}$	$\displaystyle=$	$\displaystyle\Phi(B^{i}-B_{R}^{i})-\frac{1}{\alpha}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-D_{Rk}).$		(41)

We point out that although $h^{\mu}$ in Eq. (39) may appear to be a regular vector field, it is only well defined relative to a given volume $V$ , on the boundary $\partial V$ of which we have calibrated the background field.

Next, we want to compute the corresponding transport equation for the relative helicity current $h^{\mu}$ . We now compute the four-covariant divergence of $h^{\mu}$ ,

\partial_{\mu}(\alpha\sqrt{\gamma}h^{\mu})=-2\sqrt{\gamma}(D_{i}B^{i}-D_{Ri}B_{R}^{i})\,.

(42)

which conceptually takes the same form as (19).

Based on this equation, we can now generalize the definition of relative helicity to the relativistic context,

	$\displaystyle\mathcal{H}$	$\displaystyle=\alpha\sqrt{\gamma}h^{0}\,,$		(43)
	$\displaystyle H$	$\displaystyle=\int_{V}\!{\rm d}^{3}x\,\mathcal{H}=\int_{V}\!{\rm d}^{3}x\,\sqrt{\gamma}(A_{j}+A_{Rj})(B^{j}-B_{R}^{j})\,.$		(44)

Combining this definition with the transport equation (42), we find

	$\displaystyle\partial_{t}H$	$\displaystyle=$	$\displaystyle-2\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}(D_{i}B^{i}-D_{Ri}B_{R}^{i})}$		(45)
			$\displaystyle-\oint_{\partial V}{\text{d}S_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-D_{Rk})}\,,$		(45)

which is a direct generalization of the main transport equation (10) of Ref. [26].

In the last stage of obtaining the final form of the transport equation, we want to convert Eq. (45) into an equation resembling the surface Laplace problem (11). The main obstacle is to remove the direct appearance of the vector potential.

Following Ref. [26], we perform two decompositions with respect to the potential field $D_{Ri}$ . Below we only outline the decompositions and corresponding boundary conditions, with details given in Appendix A. We first split the electric field into a solenoidal part, $\Sigma_{R}^{i}$ , and an irrotational part, $\mathcal{D}^{i}\Lambda_{R}$ , using a Helmholtz decomposition [26], i.e.,

D_{R}^{i}=\Sigma_{R}^{i}+\mathcal{D}^{i}\Lambda_{R},

(46)

with boundary constraints listed below

	$\displaystyle\mathcal{D}_{i}\Sigma_{R}^{i}=0,$		$\displaystyle{s}_{i}\Sigma_{R}^{i}\|_{\partial V}=0,$		(47)
	$\displaystyle\mathcal{D}_{i}\mathcal{D}^{i}\Lambda_{R}=\mathcal{D}_{i}D_{R}^{i},$		$\displaystyle\Lambda_{R}\|_{\partial V}=\text{constant},$		(48)

so that $\boldsymbol{\Sigma}_{R}$ is uniquely determined and the irrotational part has no contribution to the helicity transport. In this case, Eq. (45) can be written as

	$\displaystyle\partial_{t}H$	$\displaystyle=$	$\displaystyle-2\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}D_{i}B^{i}}$		(49)
		$\displaystyle-$	$\displaystyle\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-\Sigma_{Rk}}).$		(49)

This expression contains two parts, $\partial_{t}H=\partial_{t}H_{V}+\partial_{t}H_{S}$ , namely the volume change in helicity due to dissipation

\displaystyle\partial_{t}H_{V}=-2\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}D_{i}B^{i}}=-2\int_{V}{\text{d}^{3}x\,\alpha\sqrt{\gamma}E_{i}B^{i}}\,,

(50)

which vanishes for an ideal plasma, as well as the surface transport term,

\displaystyle\partial_{t}H_{S}=-\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-\Sigma_{Rk}})\,.

(51)

In turbulent systems, helicity will be produced on resistive scales due to the volume term, whereas on larger scales transport will take over. The volume term is gauge invariant and can be linked to magnetic current helicity [35]. The transport part is gauge dependent, and it is the latter we want to analyze here.

Following Ref. [26], we then decompose $\boldsymbol{D}$ and $\boldsymbol{\Sigma}_{R}$ using the Helmholtz-Hodge theorem into a normal part $\tau\hat{\boldsymbol{s}}$ , a solenoidal part $\hat{\boldsymbol{s}}\times\partial_{S}\frac{\partial\chi}{\partial t}$ , and a lamellar component (i.e. have vanishing normal component of surface curl) $\partial_{S}\zeta$ with surface covariant derivative operator defined as $\partial_{Si}=\partial_{i}-s_{i}s^{j}\partial_{j}$ [26],

	$\displaystyle D_{i}=\tau s_{i}-\varepsilon_{ijk}s^{j}\partial_{S}^{k}\frac{\partial\chi}{\partial t}-\partial_{Si}\zeta,$		(52)
	$\displaystyle\Sigma_{Ri}=\tau_{R}s_{i}-\varepsilon_{ijk}s^{j}\partial_{S}^{k}\frac{\partial\chi_{R}}{\partial t}-\partial_{Si}\zeta_{R}.$		(53)

The three components for $\boldsymbol{D}$ (and similarly for $\boldsymbol{\Sigma}_{R}$ ) can be calculated from [26],

	$\displaystyle\tau=s^{i}D_{i},$		(54)
	$\displaystyle\partial_{S}^{2}\partial_{t}\chi=-s_{i}\varepsilon^{ijk}\partial_{Sj}(D_{k}-\tau s_{k}),$		(55)
	$\displaystyle\partial_{S}^{2}\zeta=-\partial_{Si}(D^{i}-\tau s^{i}).$		(56)

Substituting these expressions into Eq. (49) and after some algebra (see Appendix C), we find

	$\displaystyle\partial_{t}H_{S}$	$\displaystyle=$	$\displaystyle\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})\partial_{k}(\zeta-\zeta_{R})}$		(57)
		$\displaystyle=$	$\displaystyle 2\oint_{\partial V}{\text{d}S\,\sqrt{\gamma}(\zeta-\zeta_{R})B_{n}}\,,$		(57)

where $B_{n}=s_{i}B^{i}$ is the magnetic field component normal to the surface $\partial V$ . Refs. [26] and [36] showed that for simple geometries (including spherical surfaces we consider in this work), $\zeta_{R}$ is a constant on the boundary. Thus, the term $\oint_{\partial V}{\text{d}S\,\sqrt{\gamma}B_{n}\zeta_{R}}$ vanishes due to solenoidal constraint on the magnetic field. The surface helicity transport is simply

\partial_{t}H_{S}=2\oint_{\partial V}{\text{d}S\,\sqrt{\gamma}\zeta B_{n}},

(58)

where $\zeta$ can be solved by the surface Laplacian Eq. (56). This is the main equation for transport of relativistic gauge-invariant helicity through a given closed surface.

III.2 Ideal magnetohydrodynamics limit

We now consider the case of a nonresistive, ideal plasma with a fluid four-velocity $u^{\mu}$ . Within the 3+1 split, we can define an induced three-velocity, $v^{i}$ , on the hypersurface, which satisfies

\displaystyle\tilde{u}^{i}\equiv\frac{u^{i}}{u^{0}}=\alpha v^{i}-\beta^{i}\,,

(59)

where we call $\tilde{u}^{i}$ the advection velocity. In ideal magnetohydrodynamics, the electric field is given as

	$\displaystyle E^{i}=$	$\displaystyle-\varepsilon^{ijk}v_{j}B_{k}\,,$		(60)
	$\displaystyle D^{i}=$	$\displaystyle-\varepsilon^{ijk}\tilde{u}_{j}B_{k}\,.$		(61)

We can further split the advection velocity into a part perpendicular to the magnetic field $\boldsymbol{\tilde{u}}_{\perp}$ , and a part parallel to the field $\boldsymbol{\tilde{u}}_{\parallel}$ . Then we can rewrite Eq. (61) into

D^{i}=-\varepsilon^{ijk}\tilde{u}_{\perp j}B_{k}\,,

(62)

where the perpendicular component can be calculated from

\tilde{u}_{\perp j}=\tilde{u}_{j}-\frac{\tilde{u}_{i}B^{i}}{B_{k}B^{k}}B_{j}\,,

(63)

We can further decompose the electric field into an emerging (em) term $D_{\rm em}^{i}$ (representing the transport of linked flux across the surface) and a shearing (sh) term $D_{\rm sh}^{i}$ (representing the twisting and tangling of footpoints by motions in the surface). These two terms depends on the normal and tangent components of the fluid advection velocity $\tilde{u}_{\perp j}$ , respectively. Therefore, we have

D_{\text{em}}^{i}=-\varepsilon^{ijk}\tilde{u}_{\perp_{n}j}(B_{nk}+B_{Sk})=-\varepsilon^{ijk}\tilde{u}_{\perp_{n}j}B_{Sk}\,,

(64)

D_{\text{sh}}^{i}=-\varepsilon^{ijk}\tilde{u}_{\perp_{S}j}(B_{nk}+B_{Sk})\,,

(65)

where $\tilde{u}_{\perp j}$ and $B_{k}$ are separated into

\tilde{u}_{\perp_{n}j}=\tilde{u}_{\perp i}s^{i}s_{j}\qquad\tilde{u}_{\perp_{S}j}=\tilde{u}_{\perp j}-\tilde{u}_{\perp_{n}j}\,,

(66)

B_{nk}=B_{i}s^{i}s_{k}\qquad B_{Sk}=B_{k}-B_{nk}\,,

(67)

Combining Eqs. (64) and (65) with Eq. (56), we obtain the equations for the corresponding lamellar part $\zeta_{\rm em}$ and $\zeta_{\rm sh}$

	$\displaystyle\partial_{S}^{2}\zeta_{\text{em}}$	$\displaystyle=\partial_{Si}\varepsilon^{ijk}\tilde{u}_{\perp_{n}j}B_{Sk}\,,$		(68)
	$\displaystyle\partial_{S}^{2}\zeta_{\text{sh}}$	$\displaystyle=\partial_{Si}\varepsilon^{ijk}\tilde{u}_{\perp_{S}j}B_{nk}\,,$		(69)

Note that $\varepsilon^{ijk}\tilde{u}_{\perp_{S}j}B_{Sk}$ in Eq. (65) is discarded since that the surface divergence has no components along the unit normal direction. Combining Eqs. (68) and (69) into the surface transport Eq. (58), we obtain the final surface transport equation for magnetic helicity

\partial_{t}H_{S}=2\oint_{\partial V}{\text{d}S\,\sqrt{\gamma}\left(\zeta_{\rm em}+\zeta_{\rm sh}\right)B_{n}}\,,

(70)

and each term can be calculated separately in a gauge invariant way.

IV General-Relativistic Helicity transport in numerical simulations

Having derived a gauge-invariant formulation of relativistic relative helicity transport, we now want to demonstrate how these expressions can be used as a diagnostic tool in numerical simulations. While the formulation we present is generic and can be applied to any simulation of turbulence, black hole accretion or other relativistic system, as a test problem we focus on neutron star mergers [37, 38]. These feature complex (non-steady state) flow structures [39], small [40, 12] and large-scale [8, 41] magnetic dynamo amplification, and present an ideal test bed to investigate global flows of magnetic helicity.

We begin by formulating the discrete problem of helicity transport on a coordinate sphere, before we apply it to a neutron star merger simulation.

IV.1 Helicity transport through spherical surfaces

In simulations, it is convenient to extract data on coordinate spheres, either because the simulations are already using spherical coordinate, or because infrastructure readily exists to monitor winds and outflows on spheres. Furthermore, the transport equations we aim to solve require the existence of a Laplace operator, which is trivially fulfilled on such a smooth surface.

In spherical polar coordinates the surface divergence operator can be written as

\partial_{Si}f^{i}=\frac{1}{r\sin{\theta}}\left[\frac{\partial}{\partial\theta}(\sin{\theta}f_{\theta})+\frac{\partial f_{\phi}}{\partial\phi}\right]\,,

(71)

whereas the Laplace operator takes the following form

\partial_{S}^{2}\zeta=\frac{1}{r^{2}}\left[\frac{\partial}{\sin{\theta}\partial\theta}\left(\sin{\theta}\frac{\partial\zeta}{\partial\theta}\right)+\frac{1}{\sin^{2}{\theta}}\frac{\partial^{2}\zeta}{\partial\phi^{2}}\right]\,.

(72)

Here $\theta$ is the longitudinal and $\phi$ is azimuthal coordinate.

We adopt a second-order finite difference method to discretize these surface operator in spherical polar coordinates. Introducing grid coordinates $(i,j)$ on the numerical lattice, we can write

	$\displaystyle\partial_{Si}f^{i}$	$\displaystyle=$	$\displaystyle\frac{1}{r}\left(\frac{f_{\theta}^{i+1,j}-f_{\theta}^{i-1,j}}{2\Delta\theta}+\cot{\theta^{i,j}}f_{\theta}^{i,j}\right)$		(73)
			$\displaystyle+\frac{1}{r\sin{\theta^{i,j}}}\frac{f_{\phi}^{i,j+1}-f_{\phi}^{i,j-1}}{2\Delta\phi}\,,$		(73)

	$\displaystyle\partial_{S}^{2}\zeta$	$\displaystyle=$	$\displaystyle\frac{1}{r^{2}}\left(\frac{\zeta^{i+1,j}-2\zeta^{i,j}+\zeta^{i-1,j}}{\Delta\theta^{2}}+\cot{\theta^{i,j}}\frac{\zeta^{i+1,j}-\zeta^{i-1,j}}{2\Delta\theta}\right)$		(74)
			$\displaystyle+\frac{1}{r^{2}\sin^{2}{\theta^{i,j}}}\frac{\zeta^{i,j+1}-2\zeta^{i,j}+\zeta^{i,j-1}}{\Delta\phi^{2}}\,,$		(74)

where $\Delta\phi$ and $\Delta\theta$ are the discrete spacings of the numerical grid.

We then numerically solve Eqs. (68) and (69). For this first test, we use $236\times 116$ points, which makes a direct inversion of the resulting linear problem easily possible. In higher-resolution simulations with finer scale helicity structures, one would have to use higher angular resolutions. The helicity fluxes then straightforwardly follow from (70).

IV.2 Helicity transport in neutron star mergers

We here present a first demonstration of the helicity transport framework in a binary neutron star merger. To this end, we have numerically evolved an equal mass binary neutron star system through merger using full numerical relativity to capture the spacetime and magnetohydrodynamics of the system. A detailed discussion of the numerical setup can be found in Appendix D.

The magnetic field evolution in a neutron star merger has been studied extensively (see,e.g., [37] for a review). Rather than giving a full account of the merger and post-merger dynamics, we only give a brief summary and defer to dedicated studies in the literature [42, 7, 43, 44], since our main focus is on helicity transport.

Starting out with a mixed poloidal and toroidal field in each initial neutron star, we evolve the system through merger. Due to the disruption dynamics, the field is stretched leading to large toroidal fields in the remnant [43], which likely dominate over the pre-merger topology [45]. This is because strong turbulent dynamo amplification during merger will likely source the bulk of the magnetic field strength [40, 12, 46]. After merger, this field produced largely at the interface of the two merging stars will redistribute itself following large scale flow patterns in the remnant [39]. During all these processes, helicity will be conserved except on smallest scales where numerical resistivity can act. This implies that the merger remnant should contain large scale helicity currents. In the following, we aim to provide a first demonstration of this.

Our simulations have insufficient numerical resolution to fully capture the magnetic field amplification dynamics described above. At the same time, large scale stretching, winding and braking dynamics [47] are fully captured, as are flow structures since these happen on macroscopic scales. This means that while our value of helicity will be off (since it is produced at the resistive scale), the bulk transport flows should be meaningful.

In order to track how helicity evolves throughout the merger, we introduce a set of nested spheres centered on the origin of the domain. Since we consider and equal mass merger (see Appendix D for details) this is where the center of mass of the remnant will be. We then compute helicity fluxes through these concentric shells (see left column of Fig. 4).

We present the results of these simulations in Fig. 2. We can see that the magnetic field after merger shows turbulent structures in both the star and the disk. At the same time we also show the gauge-dependent magnetic helicity density $\mathcal{H}$ normalized to the background field strength $B$ . We can see that in the mixed poloidal/toroidal field geometry, we start with a net helicity density inside the two stars. This helicity is then transported through merger. Finally we can see that both the stellar remnant and the disk have substantial turbulent helicity patches. However, due to the gauge dependence of this quantity, it is difficult to interpret the precise meaning of these patches. Using our gauge-invariant transport formulation, we will provide a more detailed discussion later in this Section, see also Fig. 4. As an additional diagnostic quantity, we also track the magnetic cross helicity [31]. Cross-helicity is defined via

	$\displaystyle\mathcal{H}_{c}^{\mu}$	$\displaystyle=h_{T}b^{\mu}\,,$		(75)
	$\displaystyle\nabla_{\mu}\mathcal{H}_{c}^{\mu}$	$\displaystyle=Tb^{\mu}\nabla_{\mu}s\,,$		(76)

where $b^{\mu}=\left(B^{\mu}+B^{\nu}u_{\nu}u^{\mu}\right)/(-n_{\kappa}u^{\kappa})$ is the comoving magnetic field, $h_{T}$ is the specific enthalpy, $s$ being the entropy ber baryon and $T$ being the temperature. Cross-helicity measures the amount of linkage of fluid streamlines with magnetic field lines [20].

We can see that the merger produces substantial cross-helicity throughout the remnant and disk with similar turbulent patches as we found for the magnetic helicity density (see also Fig. 3). Since the cross-helicity density does not suffer from gauge ambiguities, we will loosely use it as a reference quantity to contrast it with the magnetic helicity evolution.

From Eq. (76) we can see that cross-helicity is conserved for isentropic flows, and along fluid stream lines with constant entropy. In the merger, entropy is mainly produced during the collision with the flow patterns in the remnant being largely isentropic [39, 48]. It is therefore meaningful to additionally consider transport of volume cross-helicity, $H_{c}$ , through the remnant

\displaystyle\partial_{t}H_{c}=\partial_{t}\int_{V}{\rm d}x^{3}\sqrt{\gamma}\alpha\mathcal{H}_{c}^{0}\approx-\oint_{\partial V}{\rm d}S\,s_{i}\sqrt{\gamma}\alpha h_{T}b^{i}\,.

(77)

This expression parallels that of surface magnetic helicity transport (70), but without the need for an elliptic equation.

We are now in a position to demonstrate the helicity transport framework presented in this paper. To this end, we calculate the magnetic helicity and cross helicity fluxes on surfaces of constant radius $r/M_{\odot}=[1,3,5,8,10,15]$ from the origin (see Fig. 4). We can see that during merger the helicity begins to move throughout the remnant. We can best see this when comparing the fluxes through consecutive shells (right column, Fig. 4). We first see that there is a net outflow of helicity from the center. This is consistent with regions of high temperature (and entropy) rearranging into a torus shape after being initial produced at the collision interface [39]. Consequently, both cross helicity and magnetic helicity feature a strong inflow into these regions (green curves) and outflow from the inner regions (orange curves). The inflow/outflow of helicity also correlates with an outflow/inflow of electromagnetic energy, $E_{\rm EM}$ into these regions. Most importantly, we see a clear outflow of helicity into the outer regions of the remnant neutron star, which may affect the presence of large scale $\alpha\Omega-$ dynamos in that region [8, 41].

Overall, this leads to regions of oppositely signed helicity inside the merger remnant. This behavior does not seem apparent in the gauge dependent magnetic helicity density $\mathcal{H}$ , underlining the importance of gauge-invariant transport analysis.

V Discussion

Magnetic helicity is an essential quantity of magnetohydrodynamic flows and dynamos. For nonresistive plasmas present in many relativistic astrophysical systems it is conserved. As a topological invariant, it can affect the formation of large scale magnetic fields [20], as well as the feasibility of some mean-field dynamo models [16].

Magnetic helicity follows a continuity-type transport equation, but its interpretation is complicated by an apparent lack of gauge invariance. In essence, helicity is only meaningfully defined as a global integral over linked field lines solely confined in a given volume. The need to measure helicity in arbitrary patches of a given system therefore requires separating field lines contained in that volume, from field lines leaving it.

By generalizing the concepts of relative helicity by Berger & Field [22] and Finn & Antonsen [27] to the relativistic context, we have provided such a formulation that is gauge-invariant under electromagnetic and general-relativistic gauge transformations. In doing so, we have adopted the approach of Ref. [26] proposed for solar plasmas to formulate the transport fluxes as a two-dimensional elliptic problem. These depend only on the fluid velocity and local magnetic field, making the suitable for a broad variety of general-relativistic magnetohydrodynamics (GRMHD) codes that do not evolve the magnetic vector potential [30].

We have then applied this formulation to a neutron star merger, showing that after merger radial zones of different helicity form, which are associated with global currents redistributing small scale turbulent fields throughout the remnant.

It has recently been suggested that the outer layers of the hypermassive neutron star remnant are subject to strong magnetic field amplification from large-scale $\alpha\Omega-$ dynamos [8]. These could affect break out of the magnetic field from the star and subsequent jet and wind launching [10, 9, 8, 41]. Since the $\alpha\Omega-$ dynamo is subject to quenching depending on the local helicity value, it is imperative to understand the helicity background on which it is operating on. Our results indicate that global helicity currents are present, leading to an inflow of net helicity into these regions.

While our application to neutron star mergers has so far been mainly a proof-of-concept, we expect our formulation of helicity transport to be especially relevant for high resolution global simulations of dynamo action in neutron star mergers and black hole accretion disks. We plan to carry out such studies and analysis in future work.

Acknowledgements.

The authors are grateful to N. Vu for helpful advice concerning numerical solutions of elliptic problems. ERM is grateful for insightful discussions with J. Beattie, A. Bhattacharjee, and A. Philippov. ERM acknowledges support by the National Science Foundation under grant No. PHY-2309210. This work mainly used Delta at the National Center for Supercomputing Applications (NCSA) through allocation PHY210074 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. Additional simulations were performed on the NSF Frontera supercomputer under grant AST21006. ERM gratefully acknowledges discussions and participation at a workshop at the Kavli Institute for Theoretical Physics. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP).

Appendix A Decomposition of the reference field

The Helmholtz decomposition theorem (see Ref. [26] for an extensive discussion) guarantees that any finite, integrable, and continuously differentiable vector function defined in a simply connected volume can be uniquely decomposed into a solenoidal component (i.e. a curl of vector function) plus an irrotational component (i.e. a gradient of scalar). Performing such a decomposition with respect to the vector field $\boldsymbol{D}_{R}$ results in Eq. (46) with boundary conditions Eqs. (47) and (48). Thus, Eq. (45) turns into

	$\displaystyle\int_{V}{\text{d}^{3}x\,\frac{\partial}{\partial t}(\alpha\sqrt{\gamma}h^{0})}$	(78)
$\displaystyle=$	$\displaystyle-2\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}D_{i}B^{i}}+2\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}(\Sigma_{Ri}+\mathcal{D}_{i}\Lambda_{R})B_{R}^{i}}$
$\displaystyle-$	$\displaystyle\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-\Sigma_{Rk}-\mathcal{D}_{k}\Lambda_{R})}$

The second term on the RHS of Eq. (78) is

	$\displaystyle\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}(\Sigma_{Ri}+\mathcal{D}_{i}\Lambda_{R})B_{R}^{i}}$	(79)
$\displaystyle=$	$\displaystyle-\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}\Sigma_{Ri}\mathcal{D}^{i}\psi_{R}}+\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}\mathcal{D}_{i}\Lambda_{R}B_{R}^{i}}$
$\displaystyle=$	$\displaystyle\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}\left[\psi_{R}\mathcal{D}^{i}\Sigma_{Ri}-\mathcal{D}^{i}(\psi_{R}\Sigma_{Ri})\right]}$
	$\displaystyle+\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}\left[\mathcal{D}_{i}(\Lambda_{R}B_{R}^{i})-\Lambda_{R}\mathcal{D}_{i}B_{R}^{i}\right]}$
$\displaystyle=$	$\displaystyle\oint_{\partial V}{\text{d}S\sqrt{\gamma}s_{i}(\psi_{R}\Sigma_{R}^{i}-\Lambda_{R}B_{R}^{i})}$
$\displaystyle=$	$\displaystyle\Lambda_{R}\|_{S}\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}\mathcal{D}_{i}B_{R}^{i}}=0.$

where we use the boundary conditions in Eqs. (36), (47), and (48). We can also show that the third term on the RHS of Eq. (78) can be simplified as $\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-\Sigma_{Rk})}$ , since $\mathcal{D}_{k}\Lambda_{R}\propto s_{k}$ due to its constancy on the boundary and the antisymmetric feature of $\varepsilon^{ijk}s_{i}s_{k}$ . After this decomposition, we reach Eq. (49)

Appendix B Gauge invariance of relative helicity transport

Next we show the gauge invariance of Eq. (49). We add a guage transformation to the vector potential

A_{j}\rightarrow A_{j}+\partial_{j}\Lambda

(80)

The first term on the RHS of Eq. (49), which represent the magnetic helicity inside the probing volume, are automatically gauge-invariant. The second term, which is the surface transport ( $\frac{\partial\mathcal{H}_{S}}{\partial t}$ ), can also be shown to be gauge-invariant. For simplicity, we still use $D_{k}-D_{Rk}$ rather than $D_{k}-\Sigma_{Rk}$ to verify the gauge invariance.

$\displaystyle\frac{\partial H_{S}^{\prime}}{\partial t}$	$\displaystyle=$	$\displaystyle\frac{\partial H_{S}}{\partial t}-\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}\partial_{j}\Lambda(D_{k}-D_{Rk})}$	(81)
	$\displaystyle=$	$\displaystyle\frac{\partial H_{S}}{\partial t}-\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}\partial_{j}[\Lambda(D_{k}-D_{Rk})]}$
		$\displaystyle+\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}\Lambda\partial_{j}(D_{k}-D_{Rk})}$
	$\displaystyle=$	$\displaystyle\frac{\partial H_{S}}{\partial t}+\int_{V}{\text{d}^{3}x\,\epsilon^{ijk}\partial_{i}\partial_{j}[\Lambda(D_{k}-D_{Rk})]}$
		$\displaystyle-\oint_{\partial V}{\text{d}S\,s_{i}\Lambda\partial_{t}\left[\sqrt{\gamma}(B^{i}-B_{R}^{i})\right]}$
	$\displaystyle=$	$\displaystyle\frac{\partial H_{S}}{\partial t}$

where we use Eq. (36) and antisymmetry of three-dimensional Levi-Civita tensor.

Appendix C Helmholtz-Hodge surface decomposition

The Helmholtz-Hodge decomposition [49] asserts that any finite, square integrable vector function on a $\mathcal{C}^{2}$ surface, which is satisfied in our problem as spherical surface, can be separated into a normal part plus a surface part, which is further splitted into a solenoidal component and a lamallar (i.e. whose surface curl resides in the tangent space of the surface) component. The decomposition is presented in Eqs. (52) and (53) with the three components determined by Eqs. (54)-(56). Below we show a proof of Eq. (57)

First, we calculate a relation between the solenoidal parts $\chi$ and $\chi_{R}$ . Using the definition of $\chi$ in Eq. (55), we find

$\displaystyle\partial_{S}^{2}\partial_{t}\chi$	$\displaystyle=$	$\displaystyle-s_{i}\varepsilon^{ijk}\partial_{Sj}(D_{k}-\tau s_{k})$	(82)
	$\displaystyle=$	$\displaystyle-s_{i}\varepsilon^{ijk}(\partial_{j}-s_{j}s^{l}\partial_{l})(D_{k}-\tau s_{k})$
	$\displaystyle=$	$\displaystyle-s_{i}\varepsilon^{ijk}\partial_{j}D_{k}+s_{i}\varepsilon^{ijk}\partial_{j}(\tau s_{k})$
	$\displaystyle=$	$\displaystyle s_{i}\frac{1}{\sqrt{\gamma}}\partial_{t}(\sqrt{\gamma}B^{i})+\varepsilon^{ijk}s_{i}s_{k}\partial_{j}\tau+\tau s_{i}\varepsilon^{ijk}\partial_{j}s_{k}$
	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{\gamma}}\partial_{t}(\sqrt{\gamma}s_{i}B^{i}),$

where we use the fact that the normal vector, $s^{i}$ , is irrotational.

Similarly, for $\chi_{R}$ we have

\partial_{S}^{2}\partial_{t}\chi_{R}=\frac{1}{\sqrt{\gamma}}\partial_{t}(\sqrt{\gamma}s_{i}B_{R}^{i}).

(83)

Combining the two equation together and using the boundary condition Eq. (36), we conclude that

\partial_{S}^{2}\left(\frac{\partial\chi_{P}}{\partial t}-\frac{\partial\chi}{\partial t}\right)=0

(84)

which means

\frac{\partial\chi}{\partial t}\equiv\frac{\partial\chi_{P}}{\partial t}+C

(85)

on boundary $S$ , where $C$ is a constant.

We now substitute the decomposition of Eq. (52) and (53) into the surface helicity transport term in Eq. (57) and find

	$\displaystyle\frac{\partial H_{S}}{\partial t}=-\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(D_{k}-\Sigma_{Rk})}$	(86)
$\displaystyle=$	$\displaystyle-\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(\tau-\tau_{R})s_{k}}$
	$\displaystyle+\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})}\times$
	$\displaystyle\qquad\qquad\qquad\varepsilon_{klm}s^{l}(\partial^{m}-s^{m}s_{n}\partial^{n})\left(\frac{\partial\chi}{\partial t}-\frac{\partial\chi_{R}}{\partial t}\right)$
	$\displaystyle+\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})(\partial_{k}-s_{k}s_{l}\partial^{l})(\zeta-\zeta_{R})}$
$\displaystyle=$	$\displaystyle\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(A_{j}+A_{Rj})\partial_{k}(\zeta-\zeta_{R})}$

The first term on the right-hand-side of the second equation vanishes due to antisymmetric $\varepsilon^{ijk}s_{i}s_{k}$ , and the second term also vanishes bacause Eq. (85) ensures that the surface divergence operator $\partial^{m}-s^{m}s_{n}\partial^{n}$ produces zero when acting on a surface constant object. Further, we can show

$\displaystyle\frac{\partial H_{S}}{\partial t}$	$\displaystyle=$	$\displaystyle-\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}(\zeta-\zeta_{R})\partial_{k}(A_{j}+A_{Rj})}$	(87)
		$\displaystyle+\oint_{\partial V}{\text{d}S\,s_{i}\sqrt{\gamma}\varepsilon^{ijk}\partial_{k}[(A_{j}+A_{Rj})(\zeta-\zeta_{R})]}$
	$\displaystyle=$	$\displaystyle\oint_{\partial V}{\text{d}S\,\sqrt{\gamma}(\zeta-\zeta_{R})s_{i}(B^{i}+B_{R}^{i})}$
		$\displaystyle-\int_{V}{\text{d}^{3}x\,\sqrt{\gamma}\varepsilon^{ijk}\partial_{i}\partial_{k}[(A_{j}+A_{Rj})(\zeta-\zeta_{R})]}$
	$\displaystyle=$	$\displaystyle 2\oint_{\partial V}{\text{d}S\,\sqrt{\gamma}s_{i}B^{i}(\zeta-\zeta_{R})}$

where we use the boundary condition Eq. (36) in the last equation.

Appendix D Numerical relativity simulations

We numerically solve the coupled Einstein-GRMHD system to evolve the dynamical merger phase of a binary neutron star system. We do so by evolving the spacetime dynamics using the Z4c formulation of the Einstein equations [50, 51] in moving puncture gauge [52]. Furthermore, we solve the GRMHD equations in dynamical spacetimes [53] using a vector potential formulation in Lorenz gauge [54, 55]. This allows us to have direct access to the vector potential and to compute local quantities like the gauge-dependent helicity density, $\mathcal{H}$ .
We solve the above system using the Frankfurt/IllinoisGRMHD (FIL) code [56, 57]. FIL utilizes the EinsteinToolkit infrastructure [58], and implements fourth order unlimited finite-difference for the spacetime [59], as well as conservative finite difference for the GRMHD sector following the ECHO scheme [60]. The initial data for the binary system is chosen to be an equal mass binary with a total mass of $M=2.7\,M_{\odot}$ using the DD2 equation of state [61]. The initial data is computed using the spectral Kadath[62]/FUKA[63] framework. Our simulations use a nested domain with 8 refinement levels (see Ref. [64] for details), and a finest resolution of $\Delta x=260\,\rm m$ . The magnetic field geometry is initialized as a mixed poloidal-toroidal field, with $A_{\phi}=A^{0}_{\phi}\max\left(p-0.04\,p_{\rm max}\right)^{2}$ , $A_{z}=A^{0}_{z}\max\left(p-0.04\,p_{\rm max}\right)^{2}$ , where $p$ is the fluid pressure and $p_{\rm max}$ its maximum value inside each star. We choose the overall normalization such that we start out with a maximum of $B_{\rm max}=10^{15}\,\rm G$ inside the star, as well as a ratio of $A^{0}_{z}/A^{0}_{\phi}=0.05$ of initial toroidal to poloidal field.

References

Charbonneau [2014] P. Charbonneau, Solar dynamo theory, Annual Review of Astronomy and Astrophysics 52, 251 (2014).
Ruediger et al. [1995] G. Ruediger, D. Elstner, and T. F. Stepinski, The standard-accretion disk dynamo., Astron. & Astrophys. 298, 934 (1995).
Roberts and Glatzmaier [2000] P. H. Roberts and G. A. Glatzmaier, Geodynamo theory and simulations, Rev. Mod. Phys. 72, 1081 (2000).
Brandenburg and Ntormousi [2023] A. Brandenburg and E. Ntormousi, Galactic dynamos, Annual Review of Astronomy and Astrophysics 61, 561 (2023).
Brandenburg et al. [2012] A. Brandenburg, D. Sokoloff, and K. Subramanian, Current status of turbulent dynamo theory: from large-scale to small-scale dynamos, Space Science Reviews 169, 123 (2012).
Schekochihin [2022] A. A. Schekochihin, Mhd turbulence: a biased review, Journal of Plasma Physics 88, 155880501 (2022).
Kiuchi et al. [2018] K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Global simulations of strongly magnetized remnant massive neutron stars formed in binary neutron star mergers, Phys. Rev. D 97, 124039 (2018), arXiv:1710.01311 [astro-ph.HE] .
Kiuchi et al. [2024] K. Kiuchi, A. Reboul-Salze, M. Shibata, and Y. Sekiguchi, A large-scale magnetic field produced by a solar-like dynamo in binary neutron star mergers, Nature Astron. 8, 298 (2024), arXiv:2306.15721 [astro-ph.HE] .
Combi and Siegel [2023] L. Combi and D. M. Siegel, Jets from Neutron-Star Merger Remnants and Massive Blue Kilonovae, Phys. Rev. Lett. 131, 231402 (2023), arXiv:2303.12284 [astro-ph.HE] .
Most and Quataert [2023] E. R. Most and E. Quataert, Flares, Jets, and Quasiperiodic Outbursts from Neutron Star Merger Remnants, Astrophys. J. Lett. 947, L15 (2023), arXiv:2303.08062 [astro-ph.HE] .
Aguilera-Miret et al. [2023] R. Aguilera-Miret, C. Palenzuela, F. Carrasco, and D. Viganò, Role of turbulence and winding in the development of large-scale, strong magnetic fields in long-lived remnants of binary neutron star mergers, Phys. Rev. D 108, 103001 (2023), arXiv:2307.04837 [astro-ph.HE] .
Palenzuela et al. [2022] C. Palenzuela, R. Aguilera-Miret, F. Carrasco, R. Ciolfi, J. V. Kalinani, W. Kastaun, B. Miñano, and D. Viganò, Turbulent magnetic field amplification in binary neutron star mergers, Phys. Rev. D 106, 023013 (2022), arXiv:2112.08413 [gr-qc] .
Liska et al. [2020] M. T. P. Liska, A. Tchekhovskoy, and E. Quataert, Large-Scale Poloidal Magnetic Field Dynamo Leads to Powerful Jets in GRMHD Simulations of Black Hole Accretion with Toroidal Field, Mon. Not. Roy. Astron. Soc. 494, 3656 (2020), arXiv:1809.04608 [astro-ph.HE] .
Hogg and Reynolds [2018] J. D. Hogg and C. Reynolds, The influence of accretion disk thickness on the large-scale magnetic dynamo, Astrophys. J. 861, 24 (2018), arXiv:1805.05372 [astro-ph.HE] .
Jacquemin-Ide et al. [2023] J. Jacquemin-Ide, F. Rincon, A. Tchekhovskoy, and M. Liska, Magnetorotational dynamo can generate large-scale vertical magnetic fields in 3d grmhd simulations of accreting black holes (2023), arXiv:2311.00034 [astro-ph.HE] .
Gruzinov and Diamond [1994] A. V. Gruzinov and P. H. Diamond, Self-consistent theory of mean-field electrodynamics, Phys. Rev. Lett. 72, 1651 (1994).
Berger [1999] M. A. Berger, Introduction to magnetic helicity., Plasma Physics and Controlled Fusion 41, B167 (1999).
Taylor [1974] J. B. Taylor, Relaxation of Toroidal Plasma and Generation of Reverse Magnetic Fields, Phys. Rev. Lett. 33, 1139 (1974).
Bhattacharjee and Yuan [1995] A. Bhattacharjee and Y. Yuan, Self-Consistency Constraints on the Dynamo Mechanism, Astrophys. J. 449, 739 (1995).
Blackman [2015] E. G. Blackman, Magnetic Helicity and Large Scale Magnetic Fields: A Primer, Space Sci. Rev. 188, 59 (2015), arXiv:1402.0933 [astro-ph.SR] .
Squire and Bhattacharjee [2015] J. Squire and A. Bhattacharjee, Generation of large-scale magnetic fields by small-scale dynamo in shear flows, Phys. Rev. Lett. 115, 175003 (2015), arXiv:1506.04109 [astro-ph.SR] .
Berger and Field [1984] M. A. Berger and G. B. Field, The topological properties of magnetic helicity, Journal of Fluid Mechanics 147, 133–148 (1984).
Moraitis et al. [2019] K. Moraitis, E. Pariat, G. Valori, and K. Dalmasse, Relative magnetic field line helicity, Astron. & Astrophys. 624, A51 (2019), arXiv:1902.10410 [astro-ph.SR] .
Brandenburg et al. [2009] A. Brandenburg, S. Candelaresi, and P. Chatterjee, Small-scale magnetic helicity losses from a mean-field dynamo, Mon. Not. Roy. Astron. Soc. 398, 1414 (2009), arXiv:0905.0242 [astro-ph.SR] .
Zenati and Vishniac [2023] Y. Zenati and E. T. Vishniac, Conserving Local Magnetic Helicity in Numerical Simulations, Astrophys. J. 948, 11 (2023), arXiv:2106.06078 [physics.plasm-ph] .
Schuck and Antiochos [2019] P. W. Schuck and S. K. Antiochos, Determining the transport of magnetic helicity and free energy in the sun’s atmosphere, The Astrophysical Journal 882, 151 (2019).
Finn and Antonsen [1985] J. Finn and T. Antonsen, CoPPC (1985).
Moffatt [2014] H. K. Moffatt, Helicity and singular structures in fluid dynamics, Proceedings of the National Academy of Sciences 111, 3663 (2014), https://www.pnas.org/doi/pdf/10.1073/pnas.1400277111 .
Valori et al. [2020] G. Valori, P. Démoulin, E. Pariat, A. Yeates, K. Moraitis, and L. Linan, Additivity of relative magnetic helicity in finite volumes, Astron. & Astrophys. 643, A26 (2020), arXiv:2008.00968 [astro-ph.SR] .
Font [2008] J. A. Font, Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity, Living Rev. Rel. 11, 7 (2008).
Bekenstein [1987] J. D. Bekenstein, Helicity Conservation Laws for Fluids and Plasmas, Astrophys. J. 319, 207 (1987).
Baumgarte and Shapiro [2003] T. W. Baumgarte and S. L. Shapiro, General - relativistic MHD for the numerical construction of dynamical space - times, Astrophys. J. 585, 921 (2003), arXiv:astro-ph/0211340 .
Arnowitt et al. [2008] R. L. Arnowitt, S. Deser, and C. W. Misner, The Dynamics of general relativity, Gen. Rel. Grav. 40, 1997 (2008), arXiv:gr-qc/0405109 .
Komissarov [2004] S. S. Komissarov, Electrodynamics of black hole magnetospheres, Mon. Not. Roy. Astron. Soc. 350, 407 (2004), arXiv:astro-ph/0402403 .
Russell et al. [2019] A. J. B. Russell, P. Demoulin, G. Hornig, D. I. Pontin, and S. Candelaresi, Do Current and Magnetic Helicities Have the Same Sign?, Astrophys. J. 884, 55 (2019).
Backus [1986] G. Backus, Poloidal and toroidal fields in geomagnetic field modeling, Reviews of Geophysics 24, 75 (1986).
Baiotti and Rezzolla [2017] L. Baiotti and L. Rezzolla, Binary neutron star mergers: a review of Einstein’s richest laboratory, Rept. Prog. Phys. 80, 096901 (2017), arXiv:1607.03540 [gr-qc] .
Radice et al. [2020] D. Radice, S. Bernuzzi, and A. Perego, The Dynamics of Binary Neutron Star Mergers and GW170817, Ann. Rev. Nucl. Part. Sci. 70, 95 (2020), arXiv:2002.03863 [astro-ph.HE] .
Kastaun et al. [2016] W. Kastaun, R. Ciolfi, and B. Giacomazzo, Structure of Stable Binary Neutron Star Merger Remnants: a Case Study, Phys. Rev. D 94, 044060 (2016), arXiv:1607.02186 [astro-ph.HE] .
Kiuchi et al. [2015] K. Kiuchi, P. Cerdá-Durán, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Efficient magnetic-field amplification due to the Kelvin-Helmholtz instability in binary neutron star mergers, Phys. Rev. D 92, 124034 (2015), arXiv:1509.09205 [astro-ph.HE] .
Most [2023] E. R. Most, Impact of a mean field dynamo on neutron star mergers leading to magnetar remnants, Phys. Rev. D 108, 123012 (2023), arXiv:2311.03333 [astro-ph.HE] .
Palenzuela et al. [2015] C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. L. Caballero, E. O’Connor, and M. Anderson, Effects of the microphysical Equation of State in the mergers of magnetized Neutron Stars With Neutrino Cooling, Phys. Rev. D 92, 044045 (2015), arXiv:1505.01607 [gr-qc] .
Ciolfi et al. [2019] R. Ciolfi, W. Kastaun, J. V. Kalinani, and B. Giacomazzo, First 100 ms of a long-lived magnetized neutron star formed in a binary neutron star merger, Phys. Rev. D 100, 023005 (2019), arXiv:1904.10222 [astro-ph.HE] .
Ruiz et al. [2016] M. Ruiz, R. N. Lang, V. Paschalidis, and S. L. Shapiro, Binary Neutron Star Mergers: a jet Engine for Short Gamma-ray Bursts, Astrophys. J. Lett. 824, L6 (2016), arXiv:1604.02455 [astro-ph.HE] .
Aguilera-Miret et al. [2022] R. Aguilera-Miret, D. Viganò, and C. Palenzuela, Universality of the Turbulent Magnetic Field in Hypermassive Neutron Stars Produced by Binary Mergers, Astrophys. J. Lett. 926, L31 (2022), arXiv:2112.08406 [gr-qc] .
Chabanov et al. [2023] M. Chabanov, S. D. Tootle, E. R. Most, and L. Rezzolla, Crustal Magnetic Fields Do Not Lead to Large Magnetic-field Amplifications in Binary Neutron Star Mergers, Astrophys. J. Lett. 945, L14 (2023), arXiv:2211.13661 [astro-ph.HE] .
Shapiro [2000] S. L. Shapiro, Differential rotation in neutron stars: Magnetic braking and viscous damping, Astrophys. J. 544, 397 (2000), arXiv:astro-ph/0010493 .
Most et al. [2023] E. R. Most, A. Motornenko, J. Steinheimer, V. Dexheimer, M. Hanauske, L. Rezzolla, and H. Stoecker, Probing neutron-star matter in the lab: Similarities and differences between binary mergers and heavy-ion collisions, Phys. Rev. D 107, 043034 (2023), arXiv:2201.13150 [nucl-th] .
Hodge and Hodge [1989] W. V. D. Hodge and W. V. D. Hodge, The theory and applications of harmonic integrals (CUP Archive, 1989).
Bernuzzi and Hilditch [2010] S. Bernuzzi and D. Hilditch, Constraint violation in free evolution schemes: Comparing BSSNOK with a conformal decomposition of Z4, Phys. Rev. D 81, 084003 (2010), arXiv:0912.2920 [gr-qc] .
Hilditch et al. [2013] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy, and B. Bruegmann, Compact binary evolutions with the Z4c formulation, Phys. Rev. D 88, 084057 (2013), arXiv:1212.2901 [gr-qc] .
Alcubierre et al. [2003] M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Gauge conditions for long term numerical black hole evolutions without excision, Phys. Rev. D 67, 084023 (2003), arXiv:gr-qc/0206072 .
Duez et al. [2005] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Relativistic magnetohydrodynamics in dynamical spacetimes: Numerical methods and tests, Phys. Rev. D 72, 024028 (2005), arXiv:astro-ph/0503420 .
Etienne et al. [2010] Z. B. Etienne, Y. T. Liu, and S. L. Shapiro, Relativistic magnetohydrodynamics in dynamical spacetimes: A new AMR implementation, Phys. Rev. D 82, 084031 (2010), arXiv:1007.2848 [astro-ph.HE] .
Etienne et al. [2012] Z. B. Etienne, V. Paschalidis, Y. T. Liu, and S. L. Shapiro, Relativistic MHD in dynamical spacetimes: Improved EM gauge condition for AMR grids, Phys. Rev. D 85, 024013 (2012), arXiv:1110.4633 [astro-ph.HE] .
Most et al. [2019] E. R. Most, L. J. Papenfort, and L. Rezzolla, Beyond second-order convergence in simulations of magnetized binary neutron stars with realistic microphysics, Mon. Not. Roy. Astron. Soc. 490, 3588 (2019), arXiv:1907.10328 [astro-ph.HE] .
Etienne et al. [2015] Z. B. Etienne, V. Paschalidis, R. Haas, P. Mösta, and S. L. Shapiro, IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes, Class. Quant. Grav. 32, 175009 (2015), arXiv:1501.07276 [astro-ph.HE] .
Loffler et al. [2012] F. Loffler et al., The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics, Class. Quant. Grav. 29, 115001 (2012), arXiv:1111.3344 [gr-qc] .
Zlochower et al. [2005] Y. Zlochower, J. G. Baker, M. Campanelli, and C. O. Lousto, Accurate black hole evolutions by fourth-order numerical relativity, Phys. Rev. D 72, 024021 (2005), arXiv:gr-qc/0505055 .
Del Zanna et al. [2007] L. Del Zanna, O. Zanotti, N. Bucciantini, and P. Londrillo, ECHO: an Eulerian Conservative High Order scheme for general relativistic magnetohydrodynamics and magnetodynamics, Astron. Astrophys. 473, 11 (2007), arXiv:0704.3206 [astro-ph] .
Hempel and Schaffner-Bielich [2010] M. Hempel and J. Schaffner-Bielich, Statistical Model for a Complete Supernova Equation of State, Nucl. Phys. A 837, 210 (2010), arXiv:0911.4073 [nucl-th] .
Grandclement [2010] P. Grandclement, Kadath: A Spectral solver for theoretical physics, J. Comput. Phys. 229, 3334 (2010), arXiv:0909.1228 [gr-qc] .
Papenfort et al. [2021] L. J. Papenfort, S. D. Tootle, P. Grandclément, E. R. Most, and L. Rezzolla, New public code for initial data of unequal-mass, spinning compact-object binaries, Phys. Rev. D 104, 024057 (2021), arXiv:2103.09911 [gr-qc] .
Most and Raithel [2021] E. R. Most and C. A. Raithel, Impact of the nuclear symmetry energy on the post-merger phase of a binary neutron star coalescence, Phys. Rev. D 104, 124012 (2021), arXiv:2107.06804 [astro-ph.HE] .

General-Relativistic Gauge-Invariant Magnetic Helicity Transport: Basic Formulation and Application to Neutron Star Mergers