Electromagnetic Precursors to Gravitational-wave Events: Numerical Simulations of Flaring in Pre-merger Binary Neutron Star Magnetospheres

Possessing some of the highest densities and strongest magnetic fields in the universe neutron stars are an ideal tool to study strong gravity, nuclear physics, and high-energy plasma physics alike. The recent multimessenger observation of the neutron star merger GW170817 has demonstrated how the different observational channels (gravitational waves, kilonova afterglow, and short gamma-ray burst) can be used to constrain the properties of neutron stars, e.g., their masses, radii, and, in turn, nuclear physics beyond saturation (Abbott et al. 2018; De et al. 2018; Most et al. 2018b; Raithel 2019). While the observed electromagnetic counterparts have all been emitted following the merger, the presence of strong magnetic field configurations in radio pulsars indicates that the magnetospheres of the two stars might interact nontrivially prior to the merger. Even if the two neutron stars themselves would have inactive magnetospheres due to the spin-down over their long lifetime, the interaction of the two stars can reignite pair creation and establish a nearly force-free (Goldreich & Julian 1969) magnetosphere filled with pair plasma at the time of merger (Lyutikov 2019b). This scenario has been shown analytically to dissipate electromagnetic energy either in terms of a unipolar inductor scenario (Lai 2012; Piro 2012), where one of the stars has a higher magnetic field than its companion and also in the case of comparable magnetization (Hansen & Lyutikov 2001; Lyutikov 2019b). In addition, this scenario has also been invoked to drive powerful fireballs (Metzger & Zivancev 2016) and fast radio bursts (FRBs; Wang et al. 2016). Since the highly dynamical electromagnetic field configurations present in the inspiraling binary are too involved to be studied using purely analytical approaches a few numerical studies have been performed in order to study force-free magnetospheric interactions, either in binary black hole mergers (Palenzuela et al. 2010; Alic et al. 2012), in neutron star binaries (Palenzuela et al. 2013a, 2013b; Ponce et al. 2014), in mixed binaries (Paschalidis et al. 2013), and in collapsing neutron stars (Lehner et al. 2012; Palenzuela 2013; see also Nathanail et al. 2017; Most et al. 2018a for electrovacuum simulations). While these studies have been performed self-consistently in full general relativity they have not studied the main source of energy dissipation in current sheets, which are important sources of broadband electromagnetic emission. In this Letter we show that if all of the magnetic field dynamics is fully resolved the interaction of the magnetic fields in the binary can lead to the launching of powerful magnetic flares similar to magnetars (Kaspi & Beloborodov 2017) and coronal mass ejection in the Sun (Forbes 2000), where magnetic energy dissipation occurs in the main current sheet trailing the flare.

This work studies the emission of a electromagnetic precursor prior to the merger of a double neutron star system. We model the neutron stars as spherical conductors with a circumferential radius of 13 km and a spin axis aligned with the orbital angular momentum. The neutron stars are equipped with dipole magnetic fields having a magnetic field strength B₀ at the surface. We solve the covariant equations of general-relativistic force-free electrodynamics (Palenzuela 2013) using the newly developed GReX code (E. R. Most et. al. 2020, in preparation). We decompose the four-dimensional metric g_μν within a 3 + 1 split into ${{ds}}^{2}=\left(-{\alpha }^{2}+{\beta }_{k}{\beta }^{k}\right){{dt}}^{2}+2{\beta }_{k}{{dx}}^{k}{dt}+{\gamma }_{{ij}}{{dx}}^{i}{{dx}}^{j}$. Since, for simplicity, we only incorporate special-relativistic effects we adopt a corotating Minkowski frame α = 1, γ_ij = δ_ij, and ${\beta }_{i}=-{\varepsilon }_{{ijk}}{{\rm{\Omega }}}^{j}{x}^{k}$, where ${\boldsymbol{\Omega }}=(0,0,0,\omega )$ is the orbital angular momentum vector (see also Schiff 1939; Carrasco & Shibata 2020). The interior of the neutron star obeys the ideal-MHD condition

$$\begin{eqnarray}&&{E}^{i}=-{\varepsilon }^{{ijk}}{v}_{j}{B}_{k},\end{eqnarray} \tag{ 1 }$$

where ${v}_{i}={\varepsilon }_{{ijk}}{{\rm{\Omega }}}_{s}^{j}{x}_{s}^{k}$, Ω_s being the spin vector of the neutron star and x_s being the coordinate vector centered on the neutron star. We emphasize that adopting a corotating frame allows us to cleanly separate the orbital motion from the spin of the individual neutron stars so that only the spin velocity enters in (1). While simulations of pulsar magnetospheres typically do not adopt a fully general-relativistic framework to solve the Maxwell equations in the corotating frame (see, e.g., Bai & Spitkovsky 2010), it would be straightforward and interesting to also implement this approach within relativistic particle-in-cell codes for magnetospheric modeling (Parfrey et al. 2019; Crinquand et al. 2020). More details can be found in the Appendix. The exterior is evolved according to the force-free conditions ${E}_{i}{B}^{i}=0$ and E² < B² (Komissarov 2004). These are then imposed using a stiff constraint relaxation scheme (Alic et al. 2012) for the evolution of the electric and magnetic fields Eⁱ and Bⁱ (Baumgarte & Shapiro 2003). To handle the stiff current in the Ohm's law we use the third-order accurate implicit–explicit (IMEX) Runge–Kutta (RK) SSP3 (4, 3, 3) scheme (Pareschi & Russo 2005) and compute the explicit fluxes using a fourth-order accurate finite volume scheme (McCorquodale & Colella 2011), combining a fifth-order WENO-Z reconstruction (Borges et al. 2008) with a Rusanov Riemann solver (Rusanov 1961) similar to the approach in Most et al. (2019). We similarly enforce the ideal-MHD condition inside the star by computing a stiff instantaneous correction current that exactly enforces (1) at every substep of the IMEX-RK scheme.

The computational grid is provided by a set of nested boxes using the AMReX (Zhang et al. 2019) highly parallel adaptive mesh-refinement framework, on which GReX is built. The outermost box extends up to ≃1200 km, whereas the innermost box spans ≃40 km in all directions with a highest resolution of ≃400 m and a total number of six refinement levels. Since the angular component of the shift ${\beta }_{\phi }=-\omega R$ diverges at large distances R from the origin, we damp the shift β_ϕ ≃ R⁻² close to the boundary and impose simple third-order extrapolation boundary conditions on the electromagnetic fields.

We investigate the evolution of the common magnetosphere of a neutron star binary in close orbit to demonstrate under which conditions powerful electromagnetic flares can be expected shortly before merger. We will first study a fiducial system in which both stars are in a synchronized orbit, with one star having an additional spin Ω_s, and the energy for the flare is extracted from the rotational energy of the spinning neutron star. In order to establish that this flaring effect is also present in binaries endowed with different field configurations we in addition consider binaries where the magnetic moment of one of the stars is misaligned with the orbital angular momentum. Finally, we also establish the robustness of our results by studying the fiducial binary at various separations ranging up to 100 km. As we show below, for the case of the misaligned binary the flaring occurs even in the case of non-spinning neutron stars. Moreover, the flaring luminosity depends only on the separation between the two stars, not on the actual values of the stellar spin. This proves that flaring events are driven by the energy stored in the twisted magnetic field loop, and not by the Poynting flux driven by the orbital or rotational motion (Carrasco et al. 2019; Carrasco & Shibata 2020).

We consider two stars endowed with strong dipole magnetic fields having a field strength B₀ ≃ 10¹² G at the surface. In particular, we focus on anti-aligned field line configurations with their magnetic moments pointing in opposite directions, which result in a closed-loop magnetic field configuration. Since the field strengths of the interacting dipoles become strongest at close separation shortly before merger, we study a fiducial binary at 45 km separation, which corresponds to an orbital light cylinder ${R}_{\mathrm{LC}}^{\mathrm{orbit}}\simeq 170\,\mathrm{km}$. The light cylinder here denotes the distance from the origin of the binary beyond which it is no longer possible for the field to corotate as this would have to happen at speeds larger than the light speed (Goldreich & Julian 1969). We include differential motion in the fiducial binary by adding spin to one of the stars, corresponding to a stellar light cylinder ${R}_{\mathrm{LC}}^{\mathrm{star}}\simeq 470\,\mathrm{km}$. As outlined in the previous section we then proceed and perform a fully special-relativistic force-free simulation of this orbiting binary to illustrate how powerful electromagnetic flares can be launched in a binary system in close contact prior to merger.

After an initial transient necessary for the force-free constraint relaxation scheme to establish an initial force-free magnetic field configuration, the spinning star in the fiducial binary continuously twists the magnetic flux tube connecting the two stars. This is illustrated in Figure 1, which shows the out-of-plane magnetic field component B^ϕ in a corotating frame. The twist causes the magnetic field lines to inflate, transferring energy from the rotation of the neutron star into the magnetic field. It is important that for this process to work the stars need to be in sufficiently close contact since the twist is established by an Alfvén wave propagating between the two stars along the flux tube, which requires that the separation a of the binary is smaller than the stellar light cylinder ${R}_{\mathrm{LC}}^{\mathrm{star}}$. Assuming realistic dimensionless spins χ ≤ 0.05,³ i.e., ${R}_{\mathrm{LC}}^{\mathrm{star}}\simeq 350\,\mathrm{km}$, in the binary (Zhu et al. 2018), this will always happen during the last orbits before the merger.

Zoom In Zoom Out Reset image size

As can already be anticipated in Figure 1, at some point the built-up pressure of the toroidal magnetic field B^ϕ will be so strong that the twisted magnetic flux tube that connects the two stars, the blue and red regions in Figure 1, has to open up. This is shown in Figure 2, which presents the flare at the time when a magnetic bubble gets ejected together with a reconnecting current sheet trailing it, similar to magnetar flares (for example, in 2D force-free simulations of Parfrey et al. 2012). Looking at the out-of-plane magnetic field component it can easily be seen that magnetic islands are formed in the current sheet, suggesting that magnetic reconnection and the plasmoid instability are taking place (Loureiro et al. 2007; Bhattacharjee et al. 2009). In order to better quantify the energy contained in the outgoing bubble as well as understand the angular distribution of the emission on large scales we extract the Poynting flux $4\pi {{\boldsymbol{S}}}_{\mathrm{EM}}={\boldsymbol{E}}\times {\boldsymbol{B}}-\tfrac{1}{2}\left({E}^{2}+{B}^{2}\right){\boldsymbol{\beta }}$ on a spherical shell placed at a radius of 450 km from the origin. We then integrate the outgoing energy flux over one burst and show its angular distribution in Figure 3. We can see that the bubble begins to widen at large scales from the binary indicating its quasi-isotropic structure at infinity. Further, the energy in the bubble is about 10⁴¹ erg per surface angle. In order to quantify the amount of energy dissipated in the current sheet we consider Poynting's theorem of electromagnetic energy conservation in the simulation domain

$$\begin{eqnarray}&&{\partial }_{t}\left(\displaystyle \frac{1}{8\pi }\left({E}^{2}+{B}^{2}\right)\right)+{\rm{\nabla }}\cdot {{\boldsymbol{S}}}_{\mathrm{EM}}={{ \mathcal L }}_{\mathrm{dissipation}},\end{eqnarray} \tag{ 2 }$$

where ${{ \mathcal L }}_{\mathrm{dissipation}}$ is the dissipative luminosity, driven by the resistive terms in Ohm's law. Equation (2) states the the electromagnetic energy ${E}_{\mathrm{EM}}=\int {d}^{3}x\tfrac{1}{8\pi }\left({E}^{2}+{B}^{2}\right)$ can only change either by a transport of energy via a Poynting flux ${{\boldsymbol{S}}}_{\mathrm{EM}}$ or via dissipation ${{ \mathcal L }}_{\mathrm{dissipation}}$. In order to estimate the amount of dissipation we consider two spherical surfaces centered on the origin of the binary. The inner surface encloses just the two stars, while the outer shell is placed at a large distance from the origin. We can then estimate the amount of dissipation by computing the energy balance in the shell between the two surfaces according to Equation (2). This is shown in the left panel of Figure 4 which shows the evolution of the energy E_shell contained in the shell during one flaring event and the amount of energy dissipation ΔE_dissipation in the current sheet. We anticipate these results not to depend strongly on numerical resolution as it has been shown that the nonlinear development of the plasmoid instability causes the magnetic reconnection rate to become independent of the physical resistivity η at high Lundquist number S_L ≥ η⁻¹ ≃ 10⁴ (Loureiro et al. 2007; Bhattacharjee et al. 2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally we highlight how the peak luminosity ${{ \mathcal L }}_{\mathrm{EM}}^{\max }$ of the precursor flare for the fiducial binary scales with the orbital separation a. This is shown in Figure 4 for both the Poynting flux contained in the outgoing magnetic bubble as well as for the dissipative luminosity in the current sheet. We find a clean scaling of ${{ \mathcal L }}_{\mathrm{EM}}\propto {a}^{-7/2}$ in both cases. In order to better understand the scaling we estimate the free energy ${\rm{\Delta }}{E}_{\mathrm{twist}}$ in the twisted flux tube in the limit of small twist, ψ, following Parfrey et al. (2013),

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{twist}}\approx \displaystyle \frac{1}{8}{\psi }^{2}{u}^{3}{E}_{0}={\psi }^{2}\displaystyle \frac{{R}^{3}}{{a}^{3}}{E}_{0},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\approx \ 1.4\times {10}^{41}{B}_{12}^{2}{R}_{13}^{3}{\psi }^{2}\,\mathrm{erg},\end{eqnarray} \tag{ 4 }$$

where u characterizes the fraction of twisted field lines, u = 2R/a, and E₀ is the energy of the unperturbed dipole configuration. In addition we have introduced B₁₂ = B/10¹² G and R₁₃ = R/13 km. If we further assume that reconnection in the flaring sheet happens during the time ${\rm{\Delta }}t\simeq 2a/{v}_{\mathrm{rec}}$ (Parfrey et al. 2013), where v_rec ∼ 0.1c is the reconnection speed, we find that

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{EM}}^{\max }=\displaystyle \frac{{\rm{\Delta }}{E}_{\mathrm{twist}}}{{\rm{\Delta }}t}\approx \eta \displaystyle \frac{{\psi }^{2}}{2}{v}_{\mathrm{rec}}\displaystyle \frac{{R}^{3}}{{a}^{4}}{E}_{0}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\approx \ 4.6\times {10}^{44}\eta {B}_{12}^{2}{R}_{13}^{3}{\psi }_{\pi }^{2},\end{eqnarray} \tag{ 6 }$$

where the pre-factor η ≤ 1 quantifies the fraction of the free energy in the twisted flux tube that is available for reconnection, and flaring happens at twist value ∼π, i.e., ψ_π = ψ/π. Although not exact, this scaling is very similar to the one obtained from the simulations, and the differences might be associated with the assumption of small twist in Equation (6). We further find that the energy dissipated in the current sheet is always an order of magnitude below the energy carried away by the flare. Although the scaling found in our simulations seems very clean we caution that we have not included any form of orbital decay caused by the inspiral. Since this will be subdominant at larger separations we still believe our results to be applicable to most binaries.

While our analysis has so far been focused on studying a single system with equal magnetization and a nonzero stellar spin of one of the neutron stars we now show that the flaring effect is quite general and will occur for a variety of orbital configurations. In Figure 5 we show the flaring process for an equally magnetized binary, where both stars are non-spinning, but the left star has a magnetic moment that is misaligned by 45° with the direction of the orbital angular momentum. Whereas in the fiducial case the differential motion was caused by a relative difference in spin between the two stars. A non-spinning binary with a misaligned magnetic field will twist its common magnetosphere because of the orbital motion itself. In other words, while the energy that drives the flare in the fiducial case was provided by the rotational motion of the star, it is in this case provided purely by the orbital motion. Since we expect pulsar magnetic fields to be generally misaligned with their spin axis, this proves that magnetic flaring on the timescale of the orbital period is the most generic result of the interaction of magnetospheres with the comparable field strength of both stars.

Zoom In Zoom Out Reset image size

We have presented the first force-free electrodynamics simulations demonstrating how powerful electromagnetic flares can be launched as precursors to the neutron star merger events. We have found that these are produced by a buildup of twist in the common force-free magnetosphere of the binary system, caused by differential motion. We have shown that this can be either caused by a relative spin difference of the two stars or misalignment of the magnetosphere, which are common for pulsars.

While predicting both high-energy and coherent radio emission signatures of the magnetic flare requires first-principle kinetic plasma simulations, lessons learned in the pulsar magnetosphere research allow us to describe potential outcomes. For example, we expect most of the dissipated power in the current sheet to go into accelerated particles and, at typical magnetic field strengths ≥10⁶ G at the sheet location, to be quickly radiated away as high-energy synchrotron radiation. Kinetic simulations of relativistic radiative reconnection (similar to ones done for pulsar current sheets in Hakobyan et al. 2019) that take into account synchrotron cooling of emitting particles and pair production due to collisions of high-energy photons are needed to calculate the expected high-energy signature. While this study adopts a simplified resistive force-free prescription for the dissipation in the current sheets, the ≈10% efficiency of converting the outgoing Poynting flux into the dissipated power, ${{ \mathcal L }}_{\mathrm{dissipation}}$, in magnetic reconnection is not uncommon in full kinetic plasma simulations of magnetospheres (Philippov & Spitkovsky 2014; Cerutti et al. 2015; Brambilla et al. 2018). Future kinetic studies of flaring events in binary magnetospheres will be able to verify this conclusion (Crinquand et al. 2019).

Given relatively low expected luminosities in the high-energy band, coherent radio emission is a best bet for a potential precursor signal (Lyutikov 2019a). We anticipate two potential mechanisms: the first channel is a radio afterglow of the magnetic reconnection in the current sheet trailing the outgoing bubble, and the second channel is the synchrotron maser emission model of the outgoing magnetized bubble shocking the ambient plasma. Lyubarsky (2019) and Philippov et al. (2019) have shown that merging plasmoids in the pulsar current sheet beyond the light cylinder can produce coherent radio emission if the magnetic field strength in the upstream of the sheet is ≃10⁶–10⁸ G.⁴ Rescaling the results shown in Figures 2 and 5 we find that in our simulation this corresponds to fields strengths of B₀ ≃ 10⁸–10¹⁰ G at the surface of the star. For such a field strength the outgoing Poynting luminosity would, hence, be ${{ \mathcal L }}_{\mathrm{EM}}^{\max }\simeq {10}^{39}\mbox{--}{10}^{41}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$. Given typically very low efficiency, ≤10⁻⁴, of converting reconnecting magnetic flux into escaping coherent electromagnetic waves it makes it unlikely to expect an FRB-strength radio signal from the current sheet trailing the flare. As noted by Lyubarsky (2020) in a similar context of magnetar flares, reconnection in the collision of the escaping magnetic bubble with the large-scale magnetospheric current sheet is more likely to produce a powerful radio signal. Another coherent radio emission channel is linked to the escaping magnetic bubble. Since the ambient of the binary is not empty but filled with a highly conducting electron–positron plasma with insignificant baryon pollution, the bubble can produce a magnetized shock at some distance from the binary, which may drive a synchrotron maser instability and associated electromagnetic emission (Gallant et al. 1992). The quasi-isotropic structure of the escaping magnetized bubble that we find in this work should help with the detectability of these events. This scenario has been studied in the case of flares of young magnetars, which is one possibility to explain FRBs (Lyubarsky 2014; Beloborodov 2017; Metzger et al. 2019). We plan to explore these possibilities in the future.

One caveat of our study is that for now we have neglected effects of the inspiral motion and general-relativistic corrections, such as redshifts. These will be particularly important to understand potential transients from the last orbits before the merger, and we plan to incorporate those effects in an upcoming work. In order to better illustrate the observational prospect it will be important to cover the vast parameter space of magnetic field configurations, e.g., to explore unequal magnetization and inclination effects. This investigation will be reported in a follow-up paper.

During the final preparation of this work we became aware of Carrasco & Shibata (2020), who study the case of a single orbiting neutron star with a force-free magnetosphere using a setup very similar to the one described in this work. Their results focus on a continuous electromagnetic emission from a large-scale magnetospheric current sheet resulting from orbital motion, similar to the one that occurs in the magnetosphere of a usual rotating pulsar. While they are well applicable to systems with large orbital separations, i.e., greater than the stellar and orbital light cylinders, this work discusses the nontrivial magnetospheric interaction of two neutron stars prior to merger. It leads to the emission of powerful electromagnetic flares, which are more likely to be detectable as electromagnetic precursors of neutron star mergers (Callister et al. 2019).

E.R.M. and A.P. would like to thank Federico Carrasco, William East, Hayk Hakobyan, Luis Lehner, Yuri Levin, Maxim Lyutikov, Brian Metzger, Kohta Murase, Eliot Quataert, Luciano Rezzolla, Masaru Shibata, Lorenzo Sironi, Anatoly Spitkovsky, and James Stone for valuable and insightful discussions. E.R.M. gratefully acknowledges support and hospitality from the Simons Foundation through the pre-doctoral program at the Center for Computational Astrophysics, Flatiron Institute. This research was supported by the National Science Foundation under grant No. AST-1909458. Research at the Flatiron Institute is supported by the Simons Foundation.

In this appendix we briefly state the equations of general-relativistic electrodynamics augmented with divergence cleaning in order to maintain the D_iBⁱ = 0 and ${D}_{i}{E}^{i}=4\pi q$ constraints. Written in covariant form using the field strength tensor F^μν the augmented Maxwell equations read (Palenzuela 2013)

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left({F}^{\mu \nu }+\psi {g}^{\mu \nu }\right)=-4\pi {{ \mathcal J }}^{\nu }+{\kappa }_{\psi }\psi {n}^{\nu }\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left({}^{* }{F}^{\mu \nu }+\phi {g}^{\mu \nu }\right)={\kappa }_{\phi }{\phi n}^{\nu }\end{eqnarray} \tag{ A2 }$$

where ϕ, ψ are generalized Lagrange multipliers and κ_ϕ, κ_ψ are their damping constants and the 4-current ${{ \mathcal J }}^{\mu }={{qn}}^{\mu }+{J}^{\mu }$, where J^μ is the spatial part of the current, q is the charge density, and ${\boldsymbol{n}}={\alpha }^{-1}\left(1,-{\beta }^{i}\right)$ is the normal vector of the three-dimensional hypersurface of the spacetime foliation.

Within the 3 + 1 split of the metric introduced earlier in the text,

$$\begin{eqnarray}&&{{ds}}^{2}=\left(-{\alpha }^{2}+{\beta }_{k}{\beta }^{k}\right){{dt}}^{2}+2{\beta }_{k}{{dx}}^{k}{dt}+{\gamma }_{{ij}}{{dx}}^{i}{{dx}}^{j},\end{eqnarray} \tag{ A3 }$$

the Maxwell Equations (A1) and (A2) become (Palenzuela 2013)

$$\begin{eqnarray}&&{\partial }_{t}\left(\sqrt{\gamma }q\right)+{\partial }_{i}\left(\alpha {J}^{i}-{\beta }^{i}q\right)=0,\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left(\sqrt{\gamma }{B}^{i}\right)+{\partial }_{k}\left(-\sqrt{\gamma }{\beta }^{k}{B}^{i}+\alpha \sqrt{\gamma }{\varepsilon }^{{ikj}}{E}_{j}+\alpha \sqrt{\gamma }\phi \right)\\ \quad =\,-\sqrt{\gamma }{B}^{k}{\partial }_{k}{\beta }^{i}+\sqrt{\gamma }\phi \left({\gamma }^{{ij}}{\partial }_{j}\alpha -\alpha {\gamma }^{{jk}}{{\rm{\Gamma }}}_{{jk}}^{i}\right)\end{array}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left(\sqrt{\gamma }\phi \right)+{\partial }_{k}\left(-{\beta }^{k}\sqrt{\gamma }\phi +\alpha \sqrt{\gamma }{B}^{k}\right)\\ \quad =\,-\alpha \sqrt{\gamma }\phi K+\sqrt{\gamma }{B}^{k}{\partial }_{k}\alpha -\alpha {\kappa }_{\phi }\sqrt{\gamma }\phi ,\end{array}\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left(\sqrt{\gamma }{E}^{i}\right)+{\partial }_{k}\left(-\sqrt{\gamma }{\beta }^{k}{E}^{i}-\alpha \sqrt{\gamma }{\varepsilon }^{{ikj}}{B}_{j}+\alpha \sqrt{\gamma }\psi \right)\\ \quad =\,-\sqrt{\gamma }{E}^{k}{\partial }_{k}{\beta }^{i}+\sqrt{\gamma }\psi \left({\gamma }^{{ij}}{\partial }_{j}\alpha -\alpha {\gamma }^{{jk}}{{\rm{\Gamma }}}_{{jk}}^{i}\right)-4\pi \alpha \sqrt{\gamma }{J}^{i}\end{array}\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left(\sqrt{\gamma }\psi \right)+{\partial }_{k}\left(-{\beta }^{k}\sqrt{\gamma }\psi +\alpha \sqrt{\gamma }{E}^{k}\right)\\ \quad =\,-\alpha \sqrt{\gamma }\psi K+4\pi \alpha \sqrt{\gamma }q+\sqrt{\gamma }{E}^{k}{\partial }_{k}\alpha -\alpha {\kappa }_{\psi }\sqrt{\gamma }\psi \end{array}\end{eqnarray} \tag{ A8 }$$

Here ${{\rm{\Gamma }}}_{{ij}}^{k}$ is the metric compatible Christoffel symbol of the 3-metric γ_ij, $\gamma =\det {\gamma }_{{ij}}$ and K is the trace of the extrinsic curvature. These equations are covariant and are valid in any frame, in particular also in the corotating frame (Schiff 1939). In order to clarify the implications of this approach and similar (but not identically) to Carrasco & Shibata (2020), we adopt a corotating Minkowski frame (Schiff 1939) α = 1, ${\gamma }_{{ij}}={\delta }_{{ij}}$, and ${\beta }_{i}=-{\varepsilon }_{{ijk}}{{\rm{\Omega }}}^{j}{x}^{k}$, where ${\boldsymbol{\Omega }}=(0,0,0,\omega )$ is the orbital angular momentum vector. In order to understand the meaning of this construction we consider the advection velocity u^μ of a fluid element inside the stars.

There,

$$\begin{eqnarray}&&\displaystyle \frac{{u}^{i}}{{u}^{0}}=\alpha {v}^{i}-{\beta }^{i}={\varepsilon }^{{ijk}}{{\rm{\Omega }}}_{sj}{x}_{sk}+{\varepsilon }^{{ijk}}{{\rm{\Omega }}}_{j}{x}_{k},\end{eqnarray} \tag{ A9 }$$

where vⁱ is the 3-velocity seen by the corotating observer. In simple words, since the observer corotates with the star, the local velocity he sees does not contain the ${\boldsymbol{\Omega }}\times x$ part of the orbital motion and, hence, does not enter the electric field

$$\begin{eqnarray}&&{E}^{i}=-{\varepsilon }^{{ijk}}{v}_{j}{B}_{k}\end{eqnarray} \tag{ A10 }$$

seen by the corotating observer. Nonetheless, because we are solving the covariant form of the Maxwell equations, we can see that the actual field that enters the induction equation is ${\boldsymbol{E}}-{\boldsymbol{\beta }}\times {\boldsymbol{B}}$, and hence the rotational contribution does enter into the magnetic field evolution consistently. Since ${\boldsymbol{\Omega }}$ and ${{\boldsymbol{\Omega }}}_{{\bf{s}}}$ are aligned we can express this as

$$\begin{eqnarray}&&\displaystyle \frac{{u}^{i}}{{u}^{0}}={\varepsilon }^{{ijk}}\left({{\rm{\Omega }}}_{s}+{\rm{\Omega }}\right){z}_{j}{x}_{sk}+{\varepsilon }^{{ijk}}{{\rm{\Omega }}}_{j}{x}_{\mathrm{star}{\rm{k}}},\end{eqnarray} \tag{ A11 }$$

where ${\boldsymbol{z}}$ is the z-coordinate vector and ${{\boldsymbol{x}}}_{\mathrm{star}}={\boldsymbol{x}}-{{\boldsymbol{x}}}_{{\boldsymbol{s}}}$ is the coordinate vector of the stellar center. Hence, a fluid element inside the star would spin with $\left({{\rm{\Omega }}}_{s}+{\rm{\Omega }}\right)$ around the spin axis of the star and at the same time comove with the orbital motion of the star. If we would in addition subtract ${\varepsilon }_{{ijk}}{{\rm{\Omega }}}^{j}{x}_{\mathrm{star}}^{k}$ from the 3-velocity vⁱ, which enters the electric field inside the star via Equation (A10), we would indeed adopt a truly corotating frame in which the star no longer moves since then

$$\begin{eqnarray}&&\displaystyle \frac{{u}^{i}}{{u}^{0}}=\alpha {v}^{i}-{\beta }^{i}={\varepsilon }^{{ijk}}\left({{\rm{\Omega }}}_{s}+{\rm{\Omega }}\right){z}_{j}{x}_{sk},\end{eqnarray} \tag{ A12 }$$

which corresponds to the purely rotational motion of the star. We can also see from this expression that Ω_s = 0 implies a residual rotation in the true comoving frame, corresponding to the case of synchronized orbital motion of the stars.