New code for equilibriums and quasiequilibrium initial data of compact objects. IV. Rotating relativistic stars with mixed poloidal and toroidal magnetic fields

TL;DR

This paper introduces a new computational code for modeling strongly magnetized, rapidly rotating compact stars with mixed magnetic field components in full general relativity, revealing potential internal magnetic structures.

Contribution

The paper develops a novel code that solves Einstein-Maxwell-MHD equations for magnetized rotating stars with mixed poloidal and toroidal fields in non-circular spacetimes, advancing modeling capabilities.

Findings

Successfully computed highly magnetized star solutions

Found matter expulsion from regions of strongest toroidal fields

Conjecture of internal toroidal vacuum regions in extreme cases

Abstract

A new code for computing fully general relativistic solutions of strongly magnetized rapidly rotating compact stars is developed as a part of the COCAL (Compact Object CALculator) code. The full set of Einstein's equations, Maxwell's equations and magnetohydrodynamic equations are consistently solved assuming perfect conductivity, stationarity, and axisymmetry, and strongly magnetized solutions associated with mixed poloidal and toroidal components of magnetic fields are successfully obtained in generic (non-circular) spacetimes. We introduce the formulation of the problem and the numerical method in detail, then present examples of extremely magnetized compact star solutions and their convergence tests. It is found that, in extremely magnetized stars, the stellar matter can be expelled from the region of strongest toroidal fields. Hence we conjecture that a toroidal electro-vacuum region may appear inside of the extremely magnetized compact stars, which may seem like the neutron star becoming the strongest toroidal solenoid coil in the universe.

Tables

Table 1. Table 1: Summary of grid parameters.

$r_a$	:	Radial coordinate where the radial grids start.
$r_b$	:	Radial coordinate where the radial grids end.
$r_c$	:	Radial coordinate between $r_a$ and $r_b$ where
		the radial grid spacing changes.
$N_r$	:	Number of intervals $\Delta r_i$ in $r\in [r_a,r_b]$.
$N_r^f$	:	Number of intervals $\Delta r_i$ in $r\in [r_a,1]$.
$N_r^m$	:	Number of intervals $\Delta r_i$ in $r\in [r_a,r_c]$.
$N_{\theta }$	:	Number of intervals $\Delta {\theta }_j$ in $\theta \in [0,\pi ]$.
$N_{\varphi }$	:	Number of intervals $\Delta {\varphi }_k$ in $\varphi \in [0,2\pi ]$.
$L$	:	Order of included multipoles.

Table 2. Table 2: Grid parameters used for computing magnetized rotating compact stars. Resolution types SD12-SD3 are used for computing model P2, and SE12-SE3 are for P1 and P3. Normalized radial coordinates r a subscript 𝑟 𝑎 r_{a} , r b subscript 𝑟 𝑏 r_{b} , and r c subscript 𝑟 𝑐 r_{c} are in the unit of equatorial radius R 0 subscript 𝑅 0 R_{0} (in coordinate length).

Type	$r_a$	$r_b$	$r_c$	$N_r^f$	$N_r^m$	$N_r$	$N_{\theta }$	$N_{\varphi }$	$L$
SD12	0.0	${10}^6$	1.1	60	66	144	72	48	30
SD2	0.0	${10}^6$	1.1	80	88	192	96	48	30
SD23	0.0	${10}^6$	1.1	120	132	288	144	48	30
SD3	0.0	${10}^6$	1.1	160	176	384	192	48	30
SE12	0.0	${10}^6$	1.1	60	66	144	96	48	40
SE2	0.0	${10}^6$	1.1	80	88	192	128	48	40
SE23	0.0	${10}^6$	1.1	120	132	288	192	48	40
SE3	0.0	${10}^6$	1.1	160	176	384	256	48	20,30,40
SE3p	0.0	${10}^6$	1.1	160	176	384	256	60	50
SE3t	0.0	${10}^6$	1.1	160	176	384	384	60	50
SE3tp	0.0	${10}^6$	1.1	160	176	384	384	72	60

Table 3. Table 3: Quantities of a TOV solution in G = c = M ⊙ = 1 𝐺 𝑐 subscript 𝑀 direct-product 1 G=c=M_{\odot}=1 units for the polytropic EOS p = K ​ ρ Γ 𝑝 𝐾 superscript 𝜌 Γ p=K\rho^{\Gamma} with Γ = 2 Γ 2 \Gamma=2 and 3 3 3 . The values of maximum mass models of the corresponding EOS parameters are listed, where p c subscript 𝑝 𝑐 p_{c} and ρ c subscript 𝜌 𝑐 \rho_{c} are the pressure and the rest mass density at the center, M 0 subscript 𝑀 0 M_{0} is the rest mass, M 𝑀 M the gravitational mass, and M / R 𝑀 𝑅 M/R the compactness (a ratio of the gravitational mass to the circumferential radius). The adiabatic constant K 𝐾 K is chosen so that the value of M 0 subscript 𝑀 0 M_{0} becomes M 0 = 1.5 subscript 𝑀 0 1.5 M_{0}=1.5 at the compactness M / R = 0.2 𝑀 𝑅 0.2 M/R=0.2 . To convert a unit of ρ c subscript 𝜌 c \rho_{\rm c} to cgs, multiply the values by M ⊙ ​ ( G ​ M ⊙ / c 2 ) − 3 ≈ 6.176393 × 10 17 ​ g ​ cm − 3 subscript 𝑀 direct-product superscript 𝐺 subscript 𝑀 direct-product superscript 𝑐 2 3 6.176393 superscript 10 17 g superscript cm 3 M_{\odot}(GM_{\odot}/c^{2})^{-3}\approx 6.176393\times 10^{17}{\rm g}\ {\rm cm}^{-3} .

$\Gamma$	${\left(p/\rho \right)}_c$	${\rho }_c$	$M_0$	$M$	$M/R$	Models
$2$	$0.318244$	$0.00448412$	$1.51524$	$1.37931$	$0.214440$	P1, P2
$3$	$0.827497$	$0.00415972$	$2.24295$	$1.84989$	$0.316115$	P3

Table 4. Table 4: Parameters of functions in the integrability conditions ( 157 ) and ( 159 ), and of EOS ( 191 ), used in computing solutions presented in Fig.  1 , and Tables 5 and 6 .

Models	${\Lambda }_0$	${\Lambda }_1$	${\Lambda }_{\varphi 0}$	$b$	$c$	$\Gamma$
P1	$-3.0$	$0.3$	$2.3$	$0.2$	$0.5$	$2$
P2	$-1.7$	$0.1$	$1.7$	$0.2$	$0.5$	$2$
P3	$-0.2$	$0.3$	$1.0$	$0.2$	$0.5$	$3$

Table 5. Table 5: Selected solutions for extremely magnetized rotating compact stars are presented. Models P1 – P3 are calculated using the corresponding parameters in Table 4 . Listed quantities include the equatorial and polar radii in proper length R ¯ 0 subscript ¯ 𝑅 0 \bar{R}_{0} and R ¯ z subscript ¯ 𝑅 𝑧 \bar{R}_{z} , the ratio of the maximum values of the pressure to the rest mass density ( p / ρ ) c subscript 𝑝 𝜌 c (p/\rho)_{\rm c} , the angular velocity near the rotation axis Ω c subscript Ω c \Omega_{\rm c} , the ADM mass M ADM subscript 𝑀 ADM M_{\rm ADM} , the rest mass M 0 subscript 𝑀 0 M_{0} , the angular momentum J 𝐽 J , the virial constant I vir subscript 𝐼 vir I_{\rm vir} , and a residual of the equality of the Komar mass M K subscript 𝑀 K M_{\rm K} and the ADM mass M ADM subscript 𝑀 ADM M_{\rm ADM} . The definitions of these quantities are found in Appendix.  F (see also, cocal ). To convert a unit of length from G = c = M ⊙ = 1 𝐺 𝑐 subscript 𝑀 direct-product 1 G=c=M_{\odot}=1 to [km], multiply G ​ M ⊙ / c 2 = 1.477 𝐺 subscript 𝑀 direct-product superscript 𝑐 2 1.477 GM_{\odot}/c^{2}=1.477 [km].

Model	${\overline{R}}_0$	${\overline{R}}_z/{\overline{R}}_0$	${\left(p/\rho \right)}_c$	${\rho }_c\left[{\text{g/cm}}^3\right]$	${\Omega }_c$	$M_{\mathrm{ADM}}$	$M_0$	$J/M_{\mathrm{ADM}}^2$	$\left|1-M_K/M_{\mathrm{ADM}}\right|$
P1   ($\Gamma =2$, normal mass)	11.0609	0.71996	0.12322	$1.0717\times {10}^{15}$	0.026384	1.35908	1.46223	0.52809	$1.96\times {10}^{-5}$
P2   ($\Gamma =2$, supramassive)	11.0787	0.64818	0.25582	$2.2250\times {10}^{15}$	0.043648	1.58645	1.74178	0.57113	$4.06\times {10}^{-5}$
P3   ($\Gamma =3$, normal mass)	8.8439	0.71839	0.18830	$1.2248\times {10}^{15}$	0.033124	1.59179	1.80318	0.51282	$7.26\times {10}^{-5}$

Table 6. Table 6: Continuing from Table 5 , listed for the same solutions are the maximum values of poloidal and toroidal magnetic fields, B pol max subscript superscript 𝐵 max pol {B^{\rm max}_{\rm pol}} and B tor max subscript superscript 𝐵 max tor {B^{\rm max}_{\rm tor}} , the ratios of poloidal and toroidal magnetic field energies, ℳ pol subscript ℳ pol {\cal M}_{\rm pol} and ℳ tor subscript ℳ tor {\cal M}_{\rm tor} , and electric field energy ℳ ele subscript ℳ ele {\cal M}_{\rm ele} to the total electromagnetic field energy ℳ ℳ {\cal M} , the ratios of the kinetic, internal, and electromagnetic field energies to the gravitational energy, 𝒯 / | 𝒲 | 𝒯 𝒲 {\cal T}/|{\cal W}| , Π / | 𝒲 | Π 𝒲 \Pi/|{\cal W}| , and ℳ / | 𝒲 | ℳ 𝒲 {\cal M}/|{\cal W}| , respectively, and the virial constant I vir subscript 𝐼 vir I_{\rm vir} , and the electric charge contribution from the volume integral of the star Q M subscript 𝑄 𝑀 Q_{M} . Details of the definitions are found in Appendix.  F .

Model	$B_{\mathrm{pol}}^{\max }$[G]	$B_{\mathrm{tor}}^{\max }$[G]	${ℳ}_{\mathrm{pol}}/ℳ$	${ℳ}_{\mathrm{tor}}/ℳ$	${ℳ}_{\mathrm{ele}}/ℳ$	$𝒯/\left|𝒲\right|$	$\Pi /\left|𝒲\right|$	$ℳ/\left|𝒲\right|$	$I_{\mathrm{vir}}/\left|𝒲\right|$	$Q_M$
P1	$6.5382\times {10}^{17}$	$6.5133\times {10}^{17}$	0.93381	0.043905	0.022284	0.063996	0.29637	0.019832	$3.2202\times {10}^{-4}$	0.037041
P2	$6.2207\times {10}^{17}$	$2.2065\times {10}^{17}$	0.92609	0.033001	0.040905	0.083577	0.29918	0.001624	$2.4514\times {10}^{-4}$	0.024655
P3	$1.7797\times {10}^{18}$	$1.4487\times {10}^{18}$	0.93967	0.040486	0.019840	0.068800	0.29176	0.043991	$1.3324\times {10}^{-4}$	0.068080

Equations

$$v\cdot \omega =v^{\gamma }{\omega }_{\gamma \alpha \dots \beta },\omega \cdot v={\omega }_{\alpha \dots \beta \gamma }{v}^{\gamma }.$$

$$\text{mbox}{\pounds }_v\omega =v\cdot d\omega +d(v\cdot \omega ).$$

$$t^{\alpha }=\alpha n^{\alpha }+{\beta }^{\alpha },$$

$${\gamma }_{\alpha \beta }=g_{\alpha \beta }+n_{\alpha }{n}_{\beta },$$

$$d{s}^2$$

$${\gamma }_{ab}={\psi }^4{\tilde{\gamma }}_{ab},$$

$${\tilde{\gamma }}_{ab}=f_{ab}+h_{ab},\text{mbox}and{\tilde{\gamma }}^{ab}=f^{ab}+h^{ab},$$

$$\sqrt{-g}=\alpha {\psi }^6\sqrt{\tilde{\gamma }}=\alpha {\psi }^6\sqrt{f}.$$

$$K_{ab}=-\frac{1}{2}{\gamma }^{\alpha }{}_a{\gamma }^{\beta }{}_b\text{mbox}{\pounds }_n{\gamma }_{\alpha \beta }=-\frac{1}{2\alpha }\left({\partial }_t{\gamma }_{ab}-\text{mbox}{\pounds }_{\beta }{\gamma }_{ab}\right),$$

$$A_{ab}=K_{ab}-\frac{1}{3}{\gamma }_{ab}K.$$

$$A_{ab}$$

$$K=\frac{6}{\alpha \psi }\left(\text{mbox}{\pounds }_{\beta }\psi -{\partial }_t\psi \right)+\frac{1}{2\alpha \tilde{\gamma }}\left(\text{mbox}{\pounds }_{\beta }\tilde{\gamma }-{\partial }_t\tilde{\gamma }\right).$$

$$T^{\alpha \beta }=T_M^{\alpha \beta }+T_F^{\alpha \beta },$$

$$T_M^{\alpha \beta }=ϵu^{\alpha }{u}^{\beta }+p{q}^{\alpha \beta },$$

$$T_F^{\alpha \beta }=\frac{1}{4\pi }\left(F^{\alpha \gamma }{F}^{\beta }{}_{\gamma }-\frac{1}{4}{g}^{\alpha \beta }{F}_{\gamma \delta }{F}^{\gamma \delta }\right).$$

$$q^{\alpha \beta }=g^{\alpha \beta }+u^{\alpha }{u}^{\beta }$$

$$p=p(\rho ,s),ϵ=ϵ(\rho ,s),$$

$$F_{\alpha \beta }=(dA{)}_{\alpha \beta }:={\nabla }_{\alpha }{A}_{\beta }-{\nabla }_{\beta }{A}_{\alpha }.$$

$$(dF{)}_{\alpha \beta \gamma }:=3{\nabla }_{[\alpha }{F}_{\beta \gamma ]}={\nabla }_{\alpha }{F}_{\beta \gamma }+{\nabla }_{\beta }{F}_{\gamma \alpha }+{\nabla }_{\gamma }{F}_{\alpha \beta }=0.$$

$$\text{mbox}{\pounds }_{\eta }{g}_{\alpha \beta }=0,\text{mbox}{\pounds }_{\eta }{A}_{\alpha }=0,$$

$$\text{mbox}{\pounds }_{\eta }{u}_{\alpha }=0,\text{mbox}{\pounds }_{\eta }h=0,\text{mbox}{\pounds }_{\eta }s=0,$$

$$\left(G_{\alpha \beta }-8\pi T_{\alpha \beta }\right)n^{\alpha }{n}^{\beta }=0$$

$$\left(G_{\alpha \beta }-8\pi T_{\alpha \beta }\right){\gamma }^{\alpha }{}_a{n}_{\beta }=0$$

$$\left(G_{\alpha \beta }-8\pi T_{\alpha \beta }\right)\left({\gamma }^{\alpha \beta }+\frac{1}{2}{n}^{\alpha }{n}^{\beta }\right)=0$$

$$\left(G_{\alpha \beta }-8\pi T_{\alpha \beta }\right)\left({\gamma }^{\alpha }{}_a{\gamma }^{\beta }{}_b-\frac{1}{3}{\gamma }_{ab}{\gamma }^{\alpha \beta }\right)=0$$

$$\text{raise}4.30554pt{}^{\circ }\Delta \Phi =\mathcal{S}.$$

$$K=0\text{mbox}and\text{raise}4.30554pt{}^{\circ }{D}_b{\tilde{\gamma }}^{ab}=0,$$

$${\partial }_t{\tilde{\gamma }}^{ab}=O(r^{-3}).$$

$${\partial }_t{\tilde{\gamma }}^{ab}=0.$$

$$K=K_{\mathcal{G}}\text{mbox}and\text{raise}4.30554pt{}^{\circ }{D}_b{\tilde{\gamma }}^{ab}={\mathcal{G}}^a,$$

$$\delta {\beta }^a$$

$$\delta {\gamma }^{ab}$$

$${\beta }^a{}^{\prime }$$

Full text

I Introduction

A magnetar, a neutron star (NS) associated with very strong surface magnetic fields around $10^{14}-10^{15}$G, has become a widely accepted model for Soft Gamma Repeaters (SGR) and Anomalous X-ray pulsars (AXP) Magnetar . Although electromagnetic fields of observed magnetars are very strong, their electromagnetic energy may not be expected to dominate over internal or gravitational energies. Therefore in most theoretical models of magnetars, the electromagnetic fields are treated separately from the hydrostatic equilibrium of the compact stars as, for example, in Oron:2002gs , or they are treated as perturbations. With the perturbative techniques, general relativistic stars having purely poloidal magnetic fields and both toroidal and poloidal magnetic fields were calculated in Konno:1999 and PolTorMS , respectively. Effects of stable stratification on structures of stars with mixed poloidal-toroidal magnetic fields were included in SSMS .

However, the electromagnetic fields of newly born magnetars could be strong enough to have a comparable amount of energy, or could be highly concentrated and distributed anisotropically, so that the fields may largely alter the hydrostatic equilibrium of stars globally, or locally, respectively. From a theoretical viewpoint, it is also interesting to compute extreme solutions such as compact stars associated with the electromagnetic fields in their strongest limit, and to investigate their impact onto the hydrostatic as well as the spacetime structure RNS .

Several numerical methods have been developed in the last three decades for computing such stationary and axisymmetric equilibriums of relativistic compact stars, which are largely deformed due to strong electromagnetic fields and rapid rotation. The first success was achieved by the Meudon group (LUTH) for computing compact stars associated with poloidal magnetic fields Bocquet:1995je . Those associated with purely toroidal magnetic fields were solved by Kiuchi and Yoshida Kiuchi:2008ch , and later by Frieben and Rezzolla Frieben:2012dz . More recently, the Florence group published a series of articles for computing magnetized compact stars with purely poloidal, toroidal, as well as mixed magnetic fields Firenze . In their computations, simplified formulations for the gravitational fields have been used, which enabled systematic computations of solutions in a wide region of parameter space.

In our previous paper Uryu:2014tda , we presented preliminary results for stationary and axisymmetric equilibriums of relativistic rotating stars associated with strong electromagnetic fields, in particular, with mixed toroidal and poloidal magnetic fields. Following Uryu:2014tda , we detail below the formulations and a numerical methods for computing such equilibriums, including improvements on our earlier work Uryu:2014tda . We then present a few examples of solutions associated with extremely strong electromagnetic fields and results of convergence tests. In the newly calculated solutions, it is found that the toroidal magnetic fields concentrate near, but well below, the equatorial surface, and that the fields expel the matter when their strength becomes of order $10^{17}$G or higher for typical neutron stars. From this finding we can conjecture that a neutron star associated with such extremely strong toroidal magnetic fields may have a toroidal magneto-vacuum tunnel in it, that is, such a neutron star may become a toroidal solenoid itself.

This paper is organized as follows. The formulation for stationary and axisymmetric equilibriums of relativistic stars associated with electromagnetic fields is described in Sec. II with emphasis on the 3+1 decomposition of Maxwell’s equations and the derivation of a system of first integrals and integrability conditions for ideal magnetohydrodynamic (MHD) flows. In Sec. III, the derived formulation is further modified into the form implemented in the present numerical code, the cocal code, and then the numerical method used in the code is briefly described. In Sec. IV, three new numerical solutions calculated from the latest version of the cocal code for magnetized rotating equilibriums are presented, and their convergence test with respect to resolution and number of multipoles included in the Poisson solver are presented.

II Formulation

II.1 Summary for formulation

In the following, relativistic rotating stars associated with electromagnetic fields are modeled in the framework of a stationary and axisymmetric Einstein-Maxwell charged and magnetized perfect fluid spacetime. We assume that the equilibriums of magnetized stars satisfy the ideal MHD condition. Because of the nature of mixed poloidal and toroidal components of magnetic fields as well as a possible existence of meridional flows of matter, the spacetime is no longer circular: it is not invariant under a simultaneous inversion of $t\rightarrow-t$ and $\phi\rightarrow-\phi$ circular . To incorporate all metric components that describe such non-circular spacetimes, we apply the waveless formulation which is developed for solving initial data sets for numerical relativity simulations Shibata:2004qz ; MeudonWL ; WLBNS . The waveless formulation is based on a 3+1 decomposition and conformal decomposition of the spatial metric, which are commonly used in numerical relativity. Under appropriate gauge conditions, and time and rotational symmetries, the metric components are obtained by solving a system of elliptic partial differential equations (PDEs) on an asymptotically flat spacelike slice $\Sigma$.

An analogous formulation is also applied to recast Maxwell’s equations into 3+1 form, with the electromagnetic 1-form obtained by solving elliptic PDEs. The formulation for the electromagnetic fields is detailed below, which differs from the standard formulation from which the well-known Grad-Shafranov equation is derived.

A formulation for a system of ideal MHD equations has been discussed in our previous paper Gourgoulhon:2011gz . In Gourgoulhon:2011gz , integrability conditions to guarantee consistency of the stationary and axisymmetric system and associated set of first integrals have been derived. The basic idea of the formulation used in Uryu:2014tda as well as in the present paper is essentially the same as that of Gourgoulhon:2011gz , but an alternative choice of variables results in somewhat different set of equations to be solved. In the formulation of Gourgoulhon:2011gz , the electromagnetic 2-form $F=dA$ and its Hodge dual $\star F$ are decomposed covariantly using the 1-form basis dual to symmetry vectors $t^{\alpha}$ and $\phi^{\alpha}$, and three scalar fields which are the same in both decompositions of $F$ and $\star F$ (see e.g., Eqs. (2.35) and (2.36) in Gourgoulhon:2011gz ). An analogous decomposition is applied to the vorticity 2-form $d(h{\underline{u}})$, then, after careful algebraic manipulations, the relativistic transfield equation (a generalized form of the Grad-Shafranov equation with meridional flows) is derived.

In the following formulation, unlike Gourgoulhon:2011gz , we use the contravariant tensor $F^{\alpha\beta}$ instead of $\star F$, and an orthogonal basis of a reference flat metric defined in Sec. II.6.1 to decompose the set of equations. This choice is probably more common in formulations of numerical relativity, and hence results in a more familiar form of the equations, although redundant components remain in the equations. Another difference is that we do not reduce the number of variables by imposing axial symmetry in our formalism. This allows enough generality in the new part of the code that will enable easy extension for computing, for example, non-axisymmetric configurations of electromagnetic fields, electromagnetic standing waves, or a magnetic dipole field misaligned with the rotation axis, in the future. This also minimizes the effort to develop and debug a new code for such a rather complex problem, as computing tools having already implemented in cocal, such as its multipole moment elliptic PDE solver, can be utilized.

The 3+1 decomposition and the waveless formulation for Einstein’s equations are briefly summarized below, whose details are found, for example, in RNS ; NR ; NRBS ; Gourg2012 , and in our previous papers Shibata:2004qz ; WLBNS ; Uryu:2016dqr , respectively. The derivations of the formulations for Maxwell’s equations and the first integrals of the ideal MHD equations are presented in full detail in the following subsections. In this paper, we use abstract index notation for tensors; the Greek letters $\alpha,\beta,\gamma,...$ are for abstract 4D indices, the Latin lowercase letters $a,b,c,...$ for 3D indices, and the Latin uppercase letters $A,B,C,...$ for 2D indices.

In the above, we expressed the 2-forms, $F,\star F,dA$, and $d(h{\underline{u}})$ omitting indices. Such index-free notation may also be used with caution, in particular, when calculations involve forms and vectors. A dot denotes an inner product, that is, a contraction between adjacent indices. For example, a vector $v$ and a $p$-form $\omega$ have inner product

[TABLE]

In particular, the Cartan identity for a $p$-form $\omega$ in index-free notation is written,

[TABLE]

Certain relations in index-free notation are summarized in Appendix A.2.

II.2 Framework and notations

II.2.1 3+1 decomposition of spacetime

We consider globally hyperbolic spacetimes $({\cal M},g_{\alpha\beta})$, ${\cal M}={\mathbb{R}}\times\Sigma$, admitting two symmetries: stationarity associated with a timelike Killing vector $t^{\alpha}$ and axisymmetry associated with a spacelike rotational Killing vector $\phi^{\alpha}$. The spacetime is foliated by spacelike hypersurfaces $\Sigma_{t}=\chi_{t}(\Sigma_{0})$ parametrized by a time coordinate $t$, where $\chi_{t}$ is a diffeomorphism generated by $t^{\alpha}$ and $\Sigma_{0}$ is an initial slice. Because of the time-translation symmetry, $\Sigma_{t}$ are identical for any $t$. The spacelike vector $\phi^{\alpha}$ generates a congruence of circles in $\Sigma_{t}$ parametrized by $\phi$ whose length is $2\pi$. Those parameters $t$ and $\phi$ are chosen as coordinates.

The future pointing normal 1-form $n_{\alpha}$ is defined by $n_{\alpha}=-\alpha\nabla_{\alpha}t$, and it is related to $t^{\alpha}$ as

[TABLE]

where $\alpha$ and $\beta^{\alpha}$ are the lapse function and the shift vector, respectively, and the shift is spacelike, $\beta^{\alpha}n_{\alpha}=0$. The projection tensor to a slice $\Sigma_{t}$ is defined,

[TABLE]

and its pullback to $\Sigma_{t}$ is written $\gamma_{ab}$. Then, the metric $g_{\alpha\beta}$ in a chart $\{x^{\alpha}\}$ is split into 3+1 form in a chart $\{t,x^{a}\}$,

[TABLE]

For the spatial metric $\gamma_{ab}$, a conformal decomposition is introduced as

[TABLE]

where $\psi$ is the conformal factor, and $\tilde{\gamma}_{ab}$ the spatial conformal metric. This decomposition is specified through a condition $\tilde{\gamma}=f$, where $\tilde{\gamma}$ is the determinant of $\tilde{\gamma}_{ab}$ and $f$ the determinant of the flat metric $f_{ab}$ in a chart $\{x^{a}\}$. The differences between the spatial conformal metric and the flat metric, $h_{ab}$ and $h^{ab}$, are defined by

[TABLE]

where $\tilde{\gamma}^{ab}$ and $f^{ab}$ are the inverses of the corresponding metrics. Because of the conformal decomposition, the weight of the Levi-Civita tensor becomes

[TABLE]

We denote spatial derivative operators $D_{a}$, $\tilde{D}_{a}$, and {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a} which are compatible with spatial metrics $\gamma_{ab}$, $\tilde{\gamma}_{ab}$, and $f_{ab}$, respectively, and a spacetime derivative operator compatible with the metric $g_{\alpha\beta}$ by $\nabla_{\alpha}$.

Since we write down all field equations, equations of motion, and other associated relations, including coordinate conditions, using the flat derivative operator {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}, we have freedom to choose $\{x^{a}\}$, a coordinate system of the reference frame associated with $f_{ab}$, without changing the spacetime geometry. In this paper, as in elementary vector analysis, we only choose one of Cartesian, cylindrical, or spherical coordinates associated with a set of orthogonal bases for $\{x^{a}\}$ (see, Appendix A.1). Under a choice of orthogonal basis, a difference in the weight (8) arises from $f$ which may or may not be included in $\psi$ depending on whether one chooses a coordinate or non-coordinate basis.

The extrinsic curvature of $\Sigma_{t}$ is defined by

[TABLE]

where the Lie derivatives £ are defined on either $\cal M$ or $\Sigma_{t}$ depending on the vector to derive along, and $\partial_{t}$ is the pull back of $\mbox{\pounds}_{t}$ defined on $\cal M$ to $\Sigma_{t}$. The trace of $K_{ab}$ is written $K$, and the trace free part of $K_{ab}$ is defined by

[TABLE]

Substitution of Eq. (9) to the trace free part (10) results in conformal Killing operators with respect to $t^{\alpha}$ and $\beta^{\alpha}$,

[TABLE]

while the trace of Eq. (9) becomes

[TABLE]

As a result, under the condition $\tilde{\gamma}=f$ with an assumption $\partial_{t}\tilde{\gamma}=\partial_{t}f=0$, the time derivatives of $\tilde{\gamma}_{ab}$ and $\psi$ are separately related to $A_{ab}$ and $K$, respectively.

In actual computations, we introduce conformally rescaled $K_{ab}$ and $A_{ab}$ defined by $\tilde{K}_{ab}=\psi^{-4}K_{ab}$ and $\tilde{A}_{ab}=\psi^{-4}A_{ab}$, for which their indices are raised (lowered) with $\tilde{\gamma}^{ab}$ ($\tilde{\gamma}_{ab}$).

II.2.2 Stress energy tensors for the perfect fluid and electromagnetic fields

A strongly magnetized (and possibly charged) compact star is described by Einstein-Maxwell, charged and magnetized, perfect fluid spacetime. The stress-energy tensor $T^{\alpha\beta}$ is the sum of the perfect-fluid stress-energy tensor $T_{M}^{\alpha\beta}$ and the electromagnetic stress-energy tensor $T_{F}^{\alpha\beta}$,

[TABLE]

where $T_{M}^{\alpha\beta}$ is defined by

[TABLE]

and $T_{F}^{\alpha\beta}$ by

[TABLE]

In the definition of $T_{M}^{\alpha\beta}$, $u^{\alpha}$ is the fluid 4-velocity, $p$ the pressure, $\epsilon$ the energy density, and

[TABLE]

is the projection tensor onto a surface orthogonal to $u^{\alpha}$. Here, we assume that the fluids satisfy equations of state (EOS) of the form

[TABLE]

where $\rho$ is the baryon-mass density,111That is, $\rho:=m_{B}n$, with $n$ the number density of baryons and $m_{B}$ the mean baryon mass. and $s$ the entropy per unit baryon mass, although later we assume a simpler one-parameter EOS.

In the definition of $T_{F}^{\alpha\beta}$, the electromagnetic field 2-form $F_{\alpha\beta}$ is related to the electromagnetic potential 1-form $A_{\alpha}$ by

[TABLE]

In this Eq. (18), $(dA)_{\alpha\beta}$ is the exterior derivative of the 1-form $A_{\alpha}$. Since $F_{\alpha\beta}$ is a closed 2-form,

[TABLE]

II.2.3 Stationarity and axisymmetry

For stationary and axisymmetric systems, the Lie derivatives of field and matter variables, $\{g_{\alpha\beta},A_{\alpha},u^{\alpha},h,s\}$, along the time and axial symmetry vectors, asymptotically timelike vector $t^{\alpha}$, and a spacelike rotation vector $\phi^{\alpha}$, vanish:

[TABLE]

where $\eta^{\alpha}=t^{\alpha}$ or $\phi^{\alpha}$, and $h$ is the relativistic enthalpy defined by $h=(\epsilon+p)/\rho$.

As mentioned earlier, we use the same set of field equations for the gravity as the waveless formulation derived and used in Shibata:2004qz ; MeudonWL ; Uryu:2016dqr ; WLBNS . In this formulation, we do not impose $\phi$-symmetry explicitly onto the field equations. The waveless condition becomes a part of the time symmetry conditions imposed on the time derivatives of field variables in the inertial frame. Consequently, our formalism for solving the field equations may also be applicable for computing quasi-equilibrium solutions without axial symmetry.

II.3 Formulation for gravitational fields

II.3.1 Summary of the waveless formulation

As in the common formulations of numerical relativity, we decompose Einstein’s equations, $G_{\alpha\beta}=8\pi T_{\alpha\beta}$, into normal and transverse components with respect to the hypersurface $\Sigma_{t}$ NR ; NRBS ; Gourg2012 . In our equilibrium (or quasiequilibrium) initial data formalism for numerical relativity, we choose the following combinations of components,

[TABLE]

and formally recast them into a system of elliptic PDEs for the 3+1 variables $\{\psi,\tilde{\beta}_{a},\alpha\psi,h_{ab}\}$, respectively. In the present computations, we separate the flat Laplacians, {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\Delta:={\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}{\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D^{a} operating on these variables, represented by $\Phi$, and moving other terms to the source terms, ${\cal S}$,222 The manner of determining $K_{ab}$ in our formulation is analogous to the one used in an initial data formulation often referred to as the conformal thin-sandwich formalism York:1998hy ; Cook:2000vr ; NR ; NRBS ; Gourg2012 .

[TABLE]

The source terms ${\cal S}$ of Eqs. (23) and (24) contain a time derivative of the trace and trace-free part of extrinsic curvature $\partial_{t}K$ and $\partial_{t}A_{ab}$, respectively, and, as mentioned earlier, $K$ and $A_{ab}$ contain $\partial_{t}\psi$ and $\partial_{t}\tilde{\gamma}_{ab}$, respectively, as Eqs. (12) and (11).

In most of formulations for numerical relativity simulations, gauge conditions are dynamically imposed through the so called $\alpha$-driver and $\Gamma$-driver that determine the lapse $\alpha$ and shift $\beta^{a}$ NR ; NRBS ; Gourg2012 . In our initial data formulation, $\alpha$ and $\beta^{a}$ are part of the metric potentials to be determined, and the gauge conditions are introduced by prescribing the values of the trace $K$ and the divergence {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{b}\tilde{\gamma}^{ab}. We normally choose maximal slicing and Dirac gauge conditions,

[TABLE]

for the four coordinate conditions. A method to impose these conditions has been described in WLBNS ; Uryu:2016dqr , and is repeated in the next subsection.

In paper Shibata:2004qz , a waveless condition is derived for the gravitational fields, which results in all metric potentials, including $h_{ab}$, having Coulomb type fall off in the asymptotics under the gauge (26). Such waveless condition is to impose an asymptotic behavior on the time derivative of spatial conformal metric,

[TABLE]

In WLBNS ; Uryu:2016dqr as well as the present calculations, we impose a stronger condition

[TABLE]

As mentioned above, the time derivative terms in the gravitational field equations are $\{\partial_{t}\psi,\partial_{t}\tilde{\gamma}_{ab},\partial_{t}K,\partial_{t}A_{ab}\}$. Since the value of $K$ is fixed by the gauge condition (26), the time derivative of the conformal factor, $\partial_{t}\psi$, as seen in Eq. (12) does not appear in the field equations.333It may appear in a gauge condition (see below). The maximal slicing condition, $K=0$, is assumed to be satisfied not only instantaneously on the initial hypersurface, but also on the neighboring slices, hence $\partial_{t}K=0$. The waveless condition (28) fixes the value of $\partial_{t}\tilde{\gamma}_{ab}$. The remaining $\partial_{t}A_{ab}$, and other time derivative terms appearing in the equations of motion for the matter, may be prescribed by stationarity as in Eq. (20).444 The discussion above is to elucidate that the formulation can be applicable for quasi-stationary non-axisymmetric data. One can assume global time symmetry and discard all time derivative terms from the beginning, which results in the same set of equations used in the following computations.

The waveless formulation has been successfully applied for computing equilibriums of single rotating stars, as well as quasi-equilibrium initial data of non-axisymmetric rotating stars and binary neutron stars Uryu:2016dqr ; WLBNS by replacing time symmetry vector with that in the rotating frame (the helical Killing vector $k^{\alpha}=t^{\alpha}+\Omega\phi^{\alpha}$ BD92D94 ; Friedman:2001pf ) except for $\partial_{t}\tilde{\gamma}^{ab}$ on which the waveless condition (27) is imposed.555Under helical symmetry, non-axisymmetric initial data associated with standing gravitational waves may be calculated by imposing the symmetry also to the time derivative $\partial_{t}\tilde{\gamma}^{ab}$, then rearranging the field equations to separate the Helmholtz operator for $\tilde{\gamma}_{ab}$ in the left hand side BD92D94 ; Friedman:2001pf ; Klein:2004be ; PSWconsortium ; YBRUF06 .

The concrete form of Eq. (25) for each metric potential $\{\psi,\tilde{\beta}_{a},\alpha\psi,h_{ab}\}$ is presented in Shibata:2004qz ; WLBNS ; Uryu:2016dqr .

II.3.2 Imposition of the gauge conditions

Recently we have developed a novel formulation for imposing arbitrary gauge conditions on the waveless initial data, and successfully computed a black hole toroid system in Kerr-Schild coordinates Tsokaros:2018zlf . The maximal slice and Dirac conditions (26) are replaced by generalized gauge conditions

[TABLE]

where $K_{\cal G}$ and ${\cal G}^{a}$ are, respectively, a function and a vector given arbitrarily. Our computation for the strongly magnetized rotating star is therefore not limited to the choice (26), that is, $K_{\cal G}=0$ and ${\cal G}^{a}=0$. Since the asymptotic flatness may be imposed in most of applications for computing astrophysical compact objects, it will be convenient to choose gauge conditions which become the maximal and Dirac gauges except for the vicinity of the sources.

Taking into account the gauge invariance of the linearized metric under transformations

[TABLE]

we introduce a gauge potential $\xi$ and gauge vector potentials $\xi^{a}$ to adjust $\beta^{a}$ and $h^{ab}$ as

[TABLE]

To impose gauge conditions (29), we solve $K^{\prime}=K_{\cal G}$ and {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{b}h^{ab}{}^{\prime}={\cal G}^{a} for these gauge potentials $\xi$ and $\xi^{a}$, respectively, where $K^{\prime}$ is Eq. (12) in which $\beta^{a}{}^{\prime}$ is substituted in place of $\beta^{a}$. Then, the (primed) new variables are reconstructed accordingly with Eqs. (32) and (33), and the original variables are replaced by the new ones, $\beta^{a}{}^{\prime}\rightarrow\beta^{a}$ and $h^{ab}{}^{\prime}\rightarrow h^{ab}$, during iterations for solving the field equations. Substituting Eqs. (32) and (33) to these conditions, the gauge vector potentials $\xi$ and $\xi^{a}$ are solved from elliptic equations,

[TABLE]

and

[TABLE]

A time derivative term $\partial_{t}\psi$ in the source (35) is prescribed in computing initial data on $\Sigma_{t}$, or may be absorbed in the gauge condition $K_{\cal G}$. In the following, it is set $\partial_{t}\psi=0$. The above system of elliptic equations (34)–(37) are solved simultaneously and iteratively together with the field equations WLBNS ; Uryu:2016dqr ; Tsokaros:2018zlf .

II.4 Maxwell’s and relativistic ideal MHD equations

Hereafter in this section , we describe the formulations for solving electromagnetic fields and equilibriums of magnetized matter in detail. Maxwell’s equations are written

[TABLE]

where $j^{\alpha}$ is the electric current density. The converse of Poincaré lemma implies the existence of a potential 1-form $A_{\alpha}$, such that $F_{\alpha\beta}=(dA)_{\alpha\beta}$. By construction, the current density is conserved,

[TABLE]

From the Bianchi identity

[TABLE]

and the rest mass conservation law,

[TABLE]

the relativistic MHD-Euler equations are derived;

[TABLE]

where $(d(h{\underline{u}}))_{{\alpha\beta}}:=\nabla_{\alpha}(hu_{\beta})-\nabla_{\beta}(hu_{\alpha})$ is the canonical vorticity 2-form. The system of equations for the matter is closed by adding an EOS for the thermodynamic variables and a relation for energy transport. Since we introduce a one-parameter EOS for computing equilibriums, we do not need to consider the latter relation.

Finally, we assume that the ideal MHD condition holds:

[TABLE]

This condition implies the conservation of the flux $F_{\alpha\beta}$ as recalled briefly in Appendix A.2.3.

In our previous paper Gourgoulhon:2011gz , we have shown that, under stationarity and axisymmetry, the above system of Maxwell’s equations and ideal MHD equations (38)–(44) can be recast in a system of a single elliptic PDE for a master potential, the relativistic master transfield equation, and first integrals, where the master potential may be related to a flux function, for example, the $\phi$-component of the potential $A_{\alpha}$. In the absence of a meridional flow field, the PDE becomes the well known Grad-Shafranov equation. The formulation in Gourgoulhon:2011gz is superior to the other formulations, since the single governing equation for the master potential is derived in a fully covariant form thanks to the use of exterior calculus and the orthogonal decomposition of a tangent space into subspaces spanned by the Killing vectors $(t^{\alpha},\phi^{\alpha})$ and a remaining “meridional” spacelike 2-surface. It is also shown that the system of the transfield equation and associated first integrals is the most general form which contains all types of limiting cases including purely poloidal/toroidal magnetic fields, no-magnetic fields with meridional flows, and purely circular flows. As mentioned earlier, however, we do not follow this style of formulation presented in Gourgoulhon:2011gz , in particular we do not solve the master transfield equation, but solve Maxwell’s equations to determine all (3+1 decomposed) components of $A_{\alpha}$.

II.5 Formulation for the electromagnetic field

In this subsection, we derive a 3+1 form of Maxwell’s equations, which are recast in a set of elliptic PDEs for the components of the (3+1 decomposed) electromagnetic potential 1-form $A_{\alpha}$. Although these calculations are straightforward (see e.g., NRBS ; Knapp:2002fm ), we present them in detail since the resulting equations are implemented in cocal code for actual computations.

II.5.1 3+1 decomposition of electromagnetic fields

In this subsection, we introduce the 3+1 decomposed variables for the electromagnetic potential 1-form $A_{\alpha}$ and the Faraday tensor $F^{\alpha\beta}$, as well as a decomposition of its divergence $\nabla_{\beta}F^{\alpha\beta}$. To avoid confusion, for a given 4D object, its projection to $\Sigma_{t}$ defined on $\cal M$ is denoted with a barred symbol, and the one defined on $\Sigma_{t}$ itself is denoted using the same symbol with 4D indices being replaced by 3D ones. For example, the 4D object $A_{\alpha}$ is related to $\bar{A}_{\alpha}$ and the 3D object $A_{a}$ as follows;

[TABLE]

As usual, we omit the bar on a projected 4D object when the 4D object is spatial.

We define the 3+1 variables of the electromagnetic potential 1-form $A_{\alpha}$ and Faraday tensor $F^{\alpha\beta}$ by

[TABLE]

and $F^{\alpha\beta}n_{\alpha}n_{\beta}=0$ by anti-symmetry. Then, $A_{\alpha}$ and $F^{\alpha\beta}$ are related to their spatial components by

[TABLE]

As shown in Appendix B, the projected Faraday tensor, $F_{a}$ and $F_{ab}$, and its divergence defined on $\Sigma$ become

[TABLE]

II.5.2 3+1 decomposition of Maxwell’s equations

Using Eq. (53), the projection of Maxwell’s equations along the hypersurface normal $n^{\alpha}$ is written

[TABLE]

where $\rho_{\Sigma}$ is the projection of the current $j^{\alpha}$ along the normal $n^{\alpha}$ defined by

[TABLE]

and $F^{a}$ is derived from Eq. (51) by raising the index,

[TABLE]

Note that the charge $\rho_{\Sigma}$ is related to the time component of the 4 current as

[TABLE]

Using Eq. (54), the projection of Maxwell’s equations onto the hypersurface $\Sigma_{t}$ is written

[TABLE]

where $j_{\Sigma}^{a}$ is defined by

[TABLE]

Note that the projected current $j_{\Sigma}^{a}$ is related to the components of the 4 current $j^{\alpha}$ as 666For later convenience, we use $j^{a}$ as a spatial part of 4D current $j^{\alpha}$, hence $j^{a}\neq\gamma^{a}{}_{\alpha}j^{\alpha}$.

[TABLE]

For later use, the dual of Eq. (59) is derived,

[TABLE]

where the relation $\gamma_{ab}\mbox{\pounds}_{n}F^{b}=\mbox{\pounds}_{n}F_{a}+2K_{ab}\,F^{b}$ is used.

II.5.3 Imposition of stationarity

We assume that the electromagnetic potential 1-form $A_{\alpha}$ respects the time and rotational symmetry as discussed in Sec.II.2.3. In our formulation, we impose the stationary condition explicitly on the equations, $\mbox{\pounds}_{t}A_{\alpha}=0$. Since $t^{\alpha}=\alpha n^{\alpha}+\beta^{\alpha}$, we have

[TABLE]

and similarly for the duals $A^{a}$ and $F^{a}$. With this symmetry, Maxwell’s equations (55) and the relation (51) are rewritten

[TABLE]

and Eq. (62) becomes

[TABLE]

where $A_{ab}$ is the trace free part of the extrinsic curvature $K_{ab}$ as defined in Eq. (10).

II.5.4 Conformal decomposition and equations for electromagnetic potentials

To write down the final form of Maxwell’s equations implemented in our actual numerical code, we introduce a conformal decomposition of the spatial metric Eq. (6) with the condition $\tilde{\gamma}=f$ as explained in Sec.II.2.1. We introduce conformally rescaled quantities of the spatial electromagnetic potential 1-form and vector,

[TABLE]

and for the spatial Faraday tensor

[TABLE]

where tilded objects are rescaled quantities whose indices are raised or lowered by $\tilde{\gamma}^{ab}$ or $\tilde{\gamma}_{ab}$, respectively.

Details on the conformal decomposition of 3+1 form of Maxwell’s equations (65) and (67) are provided in Appendix C. The projection of Eq. (65) along $n^{\alpha}$ results in an elliptic PDE for $\alpha\Phi_{\Sigma}$,

[TABLE]

where the source $S$ is written

[TABLE]

The projection of Eq. (67) to $\Sigma_{t}$ results in elliptic PDEs for $A_{a}$,

[TABLE]

where the source $S_{a}$ is written

[TABLE]

To shorten the expression of the source (73), we keep $\tilde{D}_{a}$ instead of replacing it with {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a} and a connection $C^{c}_{ab}$ in some terms, where C^{c}_{ab}:=\frac{1}{2}\tilde{\gamma}^{cd}({\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}\tilde{\gamma}_{db}+{\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{b}\tilde{\gamma}_{ad}-{\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{d}\tilde{\gamma}_{ab}). The fifth term in the right-hand side (RHS) of the source (73), $\tilde{D}_{a}\tilde{D}_{b}\tilde{A}^{b}$, may be expanded as follows,

[TABLE]

and then the Coulomb gauge condition {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D^{a}A_{a}=0 may be imposed explicitly, or a simpler expression of this term may be written applying the condition $\tilde{\gamma}=f$ explicitly,

[TABLE]

II.5.5 Imposition of Coulomb gauge

As discussed in Appendix G, we have freedom to choose the spatial gauge for the spatial part of the electromagnetic potentials. We impose Coulomb gauge analogously to that for coordinate conditions discussed in Sec. II.3.2:

[TABLE]

Although we have not tested gauge conditions other than (76), we formulate the method to impose more general gauge conditions analogously to the imposition of coordinate conditions discussed in Sec. II.3.2:

[TABLE]

where $A_{\cal G}$ is a given arbitrary function.

Considering the following gauge transformation,

[TABLE]

we let $A^{\prime}_{a}$ defined by

[TABLE]

satisfy the gauge condition {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D^{a}A^{\prime}_{a}=A_{\cal G} which leads to

[TABLE]

Same as the coordinate imposition Eqs. (34)–(37), the above equation (80) is solved simultaneously and iteratively together with the field equations. Then, substituting the solution $f$ to Eq. (79), $A^{\prime}_{a}$ is calculated and $A_{a}$ is replaced by $A^{\prime}_{a}\rightarrow A_{a}$.

II.6 First integrals of relativistic ideal MHD equations

A key to a successful formulation for computing equilibriums of compact stars is to derive a set of first integrals of a system of hydrostationary equations. For ideal MHD case, including ideal MHD condition (44), all equations for the magnetized matter have to be analytically and consistently integrated. In the formulation in Gourgoulhon:2011gz , those first integrals are thoroughly used to reduce the number of variables in deriving a single master transfield equation, in particular to eliminate the current $j^{\alpha}$ in Maxwell’s equations. In the present formulation we simply solve the system of the first integrals simultaneously with the field equations.

In what follows, we analyze the rest mass conservation law (42), each of the $t$, $\phi$ and meridional $x^{A}$ components of the relativistic MHD-Euler equations (43), the ideal MHD condition (44), as well as Maxwell’s equations in its original form (39), applying $t$ and $\phi$ symmetries.

II.6.1 Coordinates and basis

To begin with, we introduce an orthogonal basis, $\{t^{\alpha},\phi^{\alpha},e_{A}^{\alpha}\}$, associated with the coordinates $t$, $\phi$ and two other spatial coordinates $x^{A}$, where the $t$ and $\phi$ coordinates are adapted to the spacetime symmetries generated by the two Killing vectors $t^{\alpha}$ and $\phi^{\alpha}$. The remaining two spatial ‘meridional’ coordinates $x^{A}$, with indices denoted by uppercase Latin letters $A,B,\cdots=1,2$, may, for example, be $(\varpi,z)$ for cylindrical, or $(r,\theta)$ for spherical coordinates, for example. These bases and natural 1-form bases generated from these coordinates are normalized as

[TABLE]

where $\delta_{A}^{B}$ is the Kronecker delta, and otherwise orthogonal. In the following sections, we will use the 4D flat metric $\eta_{\alpha\beta}$ and 3D flat metric $f_{ab}$ associated with these bases,

[TABLE]

Objects with contravariant indices are expanded using these bases, and are denoted, for example, as

[TABLE]

for the 4-velocity $u^{\alpha}$. It is understood that the last term is summed over $A=1,2$. Similarly, objects with covariant indices (such as p-forms) are expanded with respect to the basis $\{\nabla_{\alpha}t,\nabla_{\alpha}\phi,\nabla_{\alpha}x^{A}\}$.

II.6.2 Rest mass conservation equation

The densitized rest mass conservation equation is written,

[TABLE]

where we have applied the symmetries,

[TABLE]

and similarly to a term associated with $\phi^{\alpha}$. This suggests to introduce a stream function $\sqrt{-g}\Psi$ for the meridional flow fields $u^{A}$,

[TABLE]

where, $\epsilon^{AB}$ is an antisymmetric matrix with a signature $\epsilon^{12}=-1$

[TABLE]

The scalar function $\Psi$ is introduced to make explicit that the stream function $\sqrt{-g}\Psi$ is a densitized scalar.

II.6.3 Components of electromagnetic 2-form

and vorticity 2-form

We write the components of $F=dA$, and impose symmetries. For economical reason, we omit spacetime indices $\alpha,\beta,\cdots$ in the following of this section and in the next section.

$(t,\phi)$-component:

[TABLE]

$(t,x^{A})$-components:

[TABLE]

$(\phi,x^{A})$ components:

[TABLE]

$(x^{A},x^{B})$ components are obviously,

[TABLE]

or explicitly,

[TABLE]

We introduce expressions for the spatial components of the 2-form $F$ and its dual as follows,

[TABLE]

where $\epsilon_{AB}$ and $\epsilon_{A}{}^{B}$ are also anti-symmetric matrices with signatures $\epsilon_{12}=-1$ and $\epsilon_{1}{}^{2}=-1$.

Analogously, the components of the vorticity two-form $d(h{\underline{u}})$ are written as follows:

$(t,\phi)$ component:

[TABLE]

$(t,x^{A})$ components:

[TABLE]

$(\phi,x^{A})$ components:

[TABLE]

$(x^{A},x^{B})$ components:

[TABLE]

We introduce expressions for the spatial components of the 2-form $d(h{\underline{u}})$ as follows,777 The magnetic flux density $B$ and the vorticity $\omega$ are related to the Hodge dual of $F$ and $d(h{\underline{u}})$ , respectively, as $B=u\cdot\star F$, and $\vec{\underline{\omega}}=\star d(h{\underline{u}})\cdot u$.

[TABLE]

II.6.4 Ideal MHD condition

Substituting the 4-velocity in terms of a basis (86) to the ideal MHD condition $F\cdot u=0$, and applying the symmetries to the 2-form $F=dA$ as discussed in the previous section, each component is written as follows:

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

$x^{A}$-components:

[TABLE]

Substituting the stream function $\sqrt{-g}\Psi$ defined by Eq. (89) to each of the $t$ and $\phi$ components (105) and (106), and then multiplying $\rho\sqrt{-g}$, we have,

[TABLE]

These relations imply that the constant surfaces of the scalars $A_{t}$, $A_{\phi}$, and the scalar density $\sqrt{-g}\Psi$ coincide.888 This means that the stream function $\sqrt{-g}\Psi$ depends on the choice of coordinate conditions. Therefore, introducing the master potential $\Upsilon$ as an independent variable, they are written

[TABLE]

These are part of the integrability conditions.

To obtain the first integral of the $x^{A}$-components (107), we again substitute the definition of the stream function (89), and $(dA)_{AB}=\epsilon_{AB}B_{\phi}$, to rewrite

[TABLE]

Because of the conditions (110), it is rewritten

[TABLE]

where the primes $A^{\prime}_{t}$, $A^{\prime}_{\phi}$ and $[\sqrt{-g}\Psi]^{\prime}$ stands for a derivative with respect to the master potential $\Upsilon$,

[TABLE]

Therefore, we have one of the first integrals, a consistency relation for components to be satisfied,

[TABLE]

In the absence of meridional flows, $[\sqrt{-g}\Psi]^{\prime}(\Upsilon)=0$, Eq. (114) implies a relativistic version of Ferraro’s isorotation law Ferraro(1937) ,

[TABLE]

II.6.5 Maxwell’s equations

For any coordinate basis of the 1-form $\nabla_{\alpha}x$ $\ (x=t,\phi,x^{A})$, the projections (components) of the divergence of the Faraday tensor $\nabla_{\beta}F^{\alpha\beta}$ become

[TABLE]

where we have used $\nabla_{\alpha}t^{\alpha}=0$ and $\nabla_{\alpha}\phi^{\alpha}=0$ for Killing vectors $t^{\alpha}$ and $\phi^{\alpha}$.

Then, the components of Maxwell’s equations $\nabla_{\beta}F^{\alpha\beta}=4\pi j^{\alpha}$ become as follows:

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

$x^{A}$-component:

[TABLE]

Since $F^{tA}$, $F^{\phi A}$, $F^{AB}$ are components, these equations are also written

[TABLE]

where we substitute Eq. (98) to $F^{AB}$ in Eq. (122).

II.6.6 MHD-Euler equations

MHD-Euler equations (43) are also written in index free notation,

[TABLE]

Substituting the current vector

[TABLE]

and using the Cartan identity, and the $t,\phi$ symmetries, each component becomes as follows:

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

$x^{A}$-component:

[TABLE]

Analogously to the ideal-MHD conditions, after a somewhat lengthy calculation, integrability conditions and a set of first integrals of MHD-Euler equations are derived as follows:

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

These are combined and written using another function of $\Upsilon$, $\Lambda(\Upsilon)$,

[TABLE]

$x^{A}$-components:

[TABLE]

Derivations of the above Eqs. (128)-(131) are detailed in Appendix D.

For the case without meridional flow fields, $u^{A}=0$, a set of first integrals for the stationary and axisymmetric system can be derived from the above set of equations by taking a limit of the stream function, $[\sqrt{-g}\Psi](\Upsilon)\rightarrow$ constant. In this limit, the right hand side of the first integral (130) becomes finite as shown in Gourgoulhon:2011gz . In Appendix E, we present a direct proof for the same case with pure rotational flows, since our formulation is slightly different from Gourgoulhon:2011gz .

III Formulation for numerical computation and numerical method

In Sec. II, a set of elliptic PDEs for computing gravitational fields $\{\psi,\tilde{\beta}_{a},\alpha\psi,h_{ab}\}$, and electromagnetic fields $\{\alpha\Phi_{\Sigma},A_{a}\}$ of stationary and axisymmetric systems are derived from Einstein’s and Maxwell’s equations. The number of variables for gravitational fields is 11, as it is augmented with the conformal factor $\psi$ and a condition $\tilde{\gamma}=f$ is added to determine it. The number of electromagnetic potentials are 4, and 4 PDEs are derived from Maxwell’s equations (39). The apparent number of variables and equations matches, though there are four additional coordinate conditions (26) to be imposed on the metric, and a gauge condition (76) on the electromagnetic potentials. For a set of matter and electric currents, $\{h,T,s,\rho,u^{\alpha},j^{\alpha}\}$, 12 variables in total appear in a system of 10 equations in Sec. II.4 which are MHD-Euler equations (43), ideal MHD condition (44) (3 components), the continuity equations for the rest mass conservation and the current conservation, and normalization of the 4-velocity $u\cdot u=-1$. Instead of solving the equation for local thermal energy conservation and the 2-parameter EOS (17) (that is, $h=h(\rho,s)$) simultaneously with the above system, we assume 1-parameter (barotropic) EOS.999The component along $u^{\alpha}$ of the ideal MHD condition $F\cdot u=0$ is trivial. The component along $u^{\alpha}$ of the relativistic MHD-Euler equation constrains the flow to be isentropic. The apparent number of the variables and equations for the matter and current also match. This is also the case for the derived system of algebraic equations for the first integrals and integrability conditions.

As shown in previous sections, for stationary and axisymmetric ideal MHD, the system of equations for matter and currents are integrable analytically when the $t$ and $\phi$ components of the electromagnetic potential 1-form and the densitized stream function are homologous. We rewrite these equations to be solved iteratively in our numerical code.

III.1 Formulation for solving Maxwell’s equations in ideal MHD

Our formulation is to solve the 3+1 decomposed Maxwell’s equations in the form of elliptic equations for the electromagnetic 1-form. Those are Eqs. (70) and (72) and can be integrated by prescribing the current $j^{\alpha}$. The ideal MHD condition constrains the 4 components of $j^{\alpha}$. In particular, the $t$ and $\phi$ components appearing in Eq. (131) are an inseparable combination, which is a consequence of the integrability conditions (110) requiring the potentials $A_{t}$ and $A_{\phi}$ to be functions of a single master potential $\Upsilon$. There are several ways to rearrange the system of equations derived in Sec. II.6 to compute $j^{\alpha}$. In the rearrangement used in this paper, we choose that the master potential $\Upsilon$ to be equal to $A_{\phi}$. This choice is general enough to generate interesting solutions for rotating compact stars associated with mixed toroidal and poloidal magnetic fields.

III.1.1 Form of the currents

As shown in Eq. (122), the Maxwell’s equations relate $F^{AB}=\epsilon^{AB}B$ (98) with the $x^{A}$ components of the current $j^{A}$,

[TABLE]

Multiplying Eq. (132) by $A^{\prime}_{t}A^{\prime}_{\phi}$ and using the first integral of the $t$ and $\phi$ components of MHD-Euler equations (128) and (129), we have

[TABLE]

The combination of the $t$ and $\phi$-components of the current $j^{\alpha}$ has a similar form as above; from the first integral of the MHD-Euler equations (131), we have

[TABLE]

The above expressions for $j^{\alpha}$ are symmetric in $t$ and $\phi$, and they can be used for taking the limit as either $A_{t}$ or $A_{\phi}$ approaches to a constant. When $A_{\phi}$ is not constant, one can, without loss of generality, choose $\Upsilon=A_{\phi}$ because $A_{\phi}$ is an arbitrary function of the master potential $\Upsilon$. Since $A^{\prime}_{\phi}=1$ and $A^{\prime\prime}_{\phi}=0$ for this case, we can derive simpler expressions for $j^{\alpha}$ with the help of Eqs. (114) and (129);

[TABLE]

where $\delta^{AB}$ is the Kronecker delta, and $B_{A}$ and $\omega_{A}$ defined in Eqs. (96) and (103) are substituted.

III.1.2 Calculation of $j^{t}$

Among the four components of the current $j^{\alpha}$, there are three independent components; the $t$ and $\phi$ components appear to be a combination as in Eq. (136). Therefore we propose to use the $t$-component of Maxwell’s equations to determine the $t$-component of the current $j^{t}$. From Eq. (55), using the relations $\rho_{\Sigma}=-j^{\alpha}n_{\alpha}=\alpha j^{t}$, and D_{a}F^{a}=\psi^{-6}{\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}(\psi^{6}F^{a}), $j^{t}$ is written

[TABLE]

Then, we move the $j^{t}$ term in Eq. (136) to the right-hand side to isolate $j^{\phi}$, and use this $j^{\phi}$ as source of the spatial components of Maxwell’s equations (72) for evaluating $A_{a}$. This method works because $j^{t}$ is related to $A_{t}$ under a choice of $\Upsilon=A_{\phi}$ and hence $A_{t}$ is a prescribed function of $A_{\phi}$, $A_{t}=A_{t}(A_{\phi})$ on the support of ideal MHD fluid, that is $F_{a}$ in Eq. (137) is related to $A_{t}$ through Eq. (65) as

[TABLE]

since $A_{t}=-\alpha\Phi_{\Sigma}+A_{a}\beta^{a}$.

In actual computations, we also solve the $n^{\alpha}$ component of Maxwell’s equations (70) to cross-check the consistency of this method by comparing a solution and a prescribed function $A_{t}(A_{\phi})$. It is also necessary to solve (70) in the electro-vacuum spacetime outside of the compact stars because the above argument is valid only on the support of ideal MHD fluid.

III.2 Elliptic PDE solver

As discussed in a series of papers Uryu:2016dqr ; Tsokaros:2018zlf ; cocal , one of the basic concepts of the cocal code is to develop a simple and straightforward numerical method for computing data sets on a 3D slice $\Sigma$. Our idea is to formulate vectorial or tensorial elliptic PDEs in terms of Cartesian components, and apply the same elliptic PDE solver as that for the elliptic PDE of a scalar function on spherical coordinates $(r,\theta,\phi)$. Our scalar elliptic PDE solver, for example for Eq. (25), uses a multipole expansion of Green’s function

[TABLE]

where ${\bf x}$ and ${\bf x}^{\prime}$ are points in $V$, ${\bf x},{\bf x}^{\prime}\in V\subset\Sigma$,

[TABLE]

where $r_{>}:=\sup\{r,r^{\prime}\}$, $r_{<}:=\inf\{r,r^{\prime}\}$, $\epsilon_{m}=1$ for $m=0$, and $\epsilon_{m}=2$ for $m\geq 1$, and $P_{\ell}^{~{}m}(\cos\theta)$ are the associated Legendre polynomials. We will truncate the expansion in $\ell$ at a certain positive integer $L$ so that $0\leq\ell\leq L$.101010 Obviously, in the present aim for computing axisymmetric configurations, it is not necessary to expand in the azimuthal angle $\phi$. The Cartesian components of vector or tensor variables have trivial dependencies on $\phi$, which may be easily integrated analytically. We, however, keep $\phi$ dependencies in the formulation and the $\phi$ integrals as Eq. (139) in the numerical code for future extensions.

The function $\chi({\bf x})$ in Eq. (139) is a homogeneous solution, {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\Delta\chi({\bf x})=0, to be used for imposing boundary conditions on $\Phi({\bf x})$. The function $\chi({\bf x})$ may be included in the Green’s function, if the boundary is a concentric sphere on the spherical coordinate cocal . For this particular problem, that is, computations of compact stars that have flat asymptotics, among all elliptic PDEs, Eqs. (21)-(24), (34)–(37), (70)-(73) and (80), all of them except for one can be integrated setting $\chi({\bf x})$ to be constant, since errors introduced to the potential is negligible if the boundary of the computational domain is taken far enough from the source. An exception for the choice for $\chi({\bf x})$ is Eq. (70) to determine $\alpha\Phi_{\Sigma}$.

As mentioned in the previous subsection, $\Phi_{\Sigma}$ is related to $A_{t}$, and $A_{t}$ is determined from the integrability condition $A_{t}=A_{t}(A_{\phi})$ on the support of the ideal MHD fluids, and from Maxwell’s equations on electro-vacuum spacetime outside of the fluids. Because we also assume that $A_{\phi}$ as well as its derivatives are continuous across the stellar surface, while $A_{t}$ and hence $\alpha\Phi_{\Sigma}$ are continuous across the surface but their derivatives are not. Therefore, in solving $\alpha\Phi_{\Sigma}$, we impose a boundary condition not only at the boundary of the computational domain, but also at the stellar surface. Our idea to impose the boundary condition at the stellar surface is essentially the same as the one described in Sec 3.1 of Bocquet:1995je .

We assume that the stellar surface is a single valued function of spherical coordinates $\theta$ and $\phi$ as $r=R(\theta,\phi)$, whose origin $r=0$ is placed inside of the star. The homogenous solution $\chi({\bf x})$ outside of fluid support is regular at $r\rightarrow\infty$, that is, $\chi\propto r^{-\ell-1}$ for $r\geq R(\theta,\phi)$, and $\chi({\bf x})$ is determined so that Eq. (139) (in which $\Phi$ is replaced by $\alpha\Phi_{\Sigma}$) satisfies the boundary value

[TABLE]

To achieve this, we expand $\chi({\bf x})$ with coefficients $(a_{\ell m},b_{\ell m})$ as,

[TABLE]

where ${\cal Y}_{\ell}^{~{}m}(\theta)$ is defined by

[TABLE]

and the expansion in $\ell$ is also truncated at $L$ here. To determine $(a_{\ell m},b_{\ell m})$ from imposing boundary conditions, we apply the method of least squares. Writing the boundary value of the volume integral term in Eq. (139),

[TABLE]

and \chi^{\rm B}\,=\,\chi({\bf x})\big{|}_{r=R(\theta,\phi)}, we define the squared residuals,

[TABLE]

and apply the method of least squares to minimize $I$, that is, we solve a linear system,

[TABLE]

to determine a set of coefficients $(a_{\ell m},b_{\ell m})$. In the above Eq. (145), a summation is taken over all grid points $(\theta_{j},\phi_{k})$ which will be introduce in later section, and $\ell$ is truncated also here $0\leq\ell\leq L$.

III.3 Equations for the matter variables

Finally, we introduce final forms of the first integrals and integrability conditions which are suitable for the self-consistent field (SCF) iteration scheme used in our numerical method.

We assume the one-parameter EOS to have a single independent thermodynamic variable for simplicity, and assume homentropic flow, $s(\Upsilon)={\rm constant}$. Because of these choices, we have 5 independent variables for the matter $\{h,u^{\alpha}\}$. In the following, the 4-velocity $u^{\alpha}$ may also be written in 3+1 form,

[TABLE]

where $v^{\alpha}$ is the spatial component of the velocity that satisfies $v^{\alpha}\nabla_{\alpha}t=0$, and both expressions $u^{\alpha}$ and $v^{\alpha}$ (or $v^{a}$) are mixedly used. The meridional components $u^{A}$ are written $u^{A}=u^{t}v^{A}$.

Assuming a choice $\Upsilon=A_{\phi}$, a possible arrangement of equations for the SCF iteration of hydrodynamic variables becomes as follows.

For the meridional velocity $u^{A}$, the rest mass conservation (89) is used;

[TABLE]

For the $t$-component of the 4-velocity, $u^{t}$, the norm $u\cdot u=-1$ is used;

[TABLE]

For the $\phi$-component of the 4-velocity, $u^{\phi}$, $x^{A}$ component of ideal MHD condition (114) is used;

[TABLE]

For a thermodynamic variable, the enthalpy $h$, the combination of $t$ and $\phi$-components of MHD-Euler equations (130) is used;

[TABLE]

As mentioned earlier, even in the case of no meridional flows (purely rotational flows), $u^{A}=0$ and $[\sqrt{-g}\Psi](\Upsilon)={\rm constant}$, the above set of equations for matter are valid, and Eq. (151) can be used as the same form (see, Appendix D).

III.4 Assumptions for arbitrary functions

To specify a model of a rotating star, a concrete form of each arbitrary function that appears in the integrability conditions (110) and in the first integrals (128)-(130) has to be prescribed. We partly follow a choice made in YYE06 for these functions, which are used in our previous paper Uryu:2014tda , but we also introduce new functional forms below. As mentioned in Sec.III.1.1, we choose $\Upsilon=A_{\phi}$ as the independent variable instead of the master potential $\Upsilon$.

III.4.1 A smoothed step function

We introduce a class of a two parameter sigmoid function $\Xi^{\prime}(x;b,c)$ that varies from [math] to $1$ in a region $x\in[0,1]$ whose transition width is $b$ $(0<b<1)$, and transition center $c$ $(0<c<1)$,

[TABLE]

and its integral $\Xi(x;b,c)$,

[TABLE]

We make use of these functions in a region where $A_{\phi}$ varies on the fluid support and its contour is closed as will be explained later. Functions $\Xi^{\prime}(A_{\phi})$ and $\Xi(A_{\phi})$ are defined by

[TABLE]

and

[TABLE]

The smoothed step function $\Xi^{\prime}(A_{\phi})$ varies on a region $A_{\phi}\in[A_{\phi,\rm S}^{\rm max},A_{\phi}^{\rm max}]$, where $A_{\phi}^{\rm max}$ and $A_{\phi,\rm S}^{\rm max}$ are the maximum values of $A_{\phi}$ on the fluid support, and that on the stellar surface, respectively, Here, $A_{\phi}^{\rm max}>A_{\phi,\rm S}^{\rm max}$ is assumed. Note that Eqs. (154) and (155) are not a mere substitution of

[TABLE]

to Eqs. (152) and (153), since the prime of Eq. (154) is with respect to $A_{\phi}$, not $x$, and a constant of integration in Eq. (155) is chosen appropriately.

The parameters $(b,c)$ determine the width and position of the transition of $\Xi^{\prime}$, and are set to be $(b,c)=(0.2,0.5)$ in the following applications. Other types of smooth step functions such as those made from Hermite interpolation polynomials could be used in the same manner.

III.4.2 Models

For the function $\Lambda(A_{\phi})$, we choose

[TABLE]

where $\Lambda_{0}$, $\Lambda_{1}$, and ${\cal E}$ are constants. $\Lambda_{0}$ and $\Lambda_{1}$ are set by hand, while ${\cal E}$ is calculated from a condition to be imposed on a solution. In case of pure rotational flows without magnetic fields, the constant ${\cal E}$ agrees with the injection energy RNS .

For the function $A_{t}(A_{\phi})$, we choose

[TABLE]

where $C_{e}$ is a constant that relates to the net electric charge of the star, and $\Omega_{\rm c}$ is a constant. As discussed in YYE06 , the choice of the first term corresponds to the rigid rotation in the limits of $[\sqrt{-g}\Psi]^{\prime}\rightarrow 0$ or $B_{\phi}\rightarrow 0$, because of relation (115).

The current (135) and (136) involve terms with a derivative of function $[\sqrt{-g}\Lambda_{\phi}](A_{\phi})$ coupled with the magnetic fields. Since we assume no electric current outside of the star, $[\sqrt{-g}\Lambda_{\phi}]^{\prime}(A_{\phi})$ has to vanish outside of the fluid support. This is why we prepare a smooth function that varies between $[0,1]$ in a region $A_{\phi}\in[A_{\phi,\rm S}^{\rm max},A_{\phi}^{\rm max}]$ as in Sec. III.4.1. We choose a smooth function,

[TABLE]

where the parameter $\Lambda_{\phi 0}$ is the range of the function, $[\sqrt{-g}\Lambda_{\phi}]^{\prime}(A_{\phi})\in[0,\Lambda_{\phi 0}]$ set by hand.

In the later sections, we only present solutions without the meridional circulation flows, hence for $[\sqrt{-g}\Psi](A_{\phi})$, we set

[TABLE]

We have also tested a few models for $[\sqrt{-g}\Psi](A_{\phi})$, and successfully computed solutions with meridional flows, although so far we have calculated solutions whose meridional flows do not affect equilibrium of the stars. For example, we may choose the same form as Eqs. (159) and (160),

[TABLE]

where $a_{\Psi}$ is a parameter to be set by hand.

III.4.3 Differentially rotating models

When magnetic fields and meridional flows exist inside of compact stars, Eq. (150) implies that the stellar rotation $\Omega:=u^{\phi}/u^{t}$ becomes inevitably differential in general, because a combination $B_{\phi}/\rho\sqrt{-g}$ is not a function of $A_{\phi}$ or $\Upsilon$. When there is no meridional flow $u^{A}=0$, $[\sqrt{-g}\Psi]^{\prime}=0$, on the other hand, the form of the function $A_{t}$ (158) results in a uniform rotation as mentioned.

It seems that the latter case with no meridional flows $[\sqrt{-g}\Psi]^{\prime}=0$ is sometimes misinterpreted in the literature, as stated in Bocquet:1995je , that only the uniform rotation $\Omega={\rm constant}$ is allowed in this case. This statement seems to have been made because a distribution of $A_{\phi}$ usually becomes toroidal and hence such a toroidal differential rotation $\Omega(A_{\phi})$ was considered to be unnatural. Such differential rotation laws in which $\Omega$ depends on $\Upsilon$ (or $A_{\phi}$) are, however, allowed mathematically, and may not necessarily be too unrealistic to be rejected. For example, one can assume moderate, or weak, differential rotations,

[TABLE]

by setting $\max|\delta\Omega(\Upsilon)|$ to be a few tens of %, or less, of $\Omega_{\rm c}$. Various rotation laws can also be used for the case with meridional flow $u^{A}\neq 0$ (that is, $[\sqrt{-g}\Psi]^{\prime}\neq 0$), but it is more likely that some kind of instability such as the magnetorotational instability may be induced.

III.4.4 Other models

In our previous paper Uryu:2014tda , the functional form for $A_{t}(A_{\phi})$ was chosen the same as Eq. (158), and for the function $\Lambda(A_{\phi})$, $\Lambda_{0}$ was set to be zero in Eq. (157). Our previous choices for $[\sqrt{-g}\Lambda_{\phi}](A_{\phi})$ and $[\sqrt{-g}\Psi](A_{\phi})$ in Uryu:2014tda were taken from those used in YYE06 . For $[\sqrt{-g}\Lambda_{\phi}]$, we have chosen

[TABLE]

where values of the constant coefficient $a$ and index $k$ are set by hand, and $\Theta(x)$ is the Heaviside function. In YYE06 , it was found that the solutions have comparable strength in poloidal and toroidal components of magnetic fields when the index is about $k=0.1$. We replace Eqs. (166) and (167) with Eqs. (159) and (160) to bring a smoothness as well as better control in the behavior of the functional forms, although either choice of the function reproduce qualitatively the same set of solutions. Also in Uryu:2014tda , for $[\sqrt{-g}\Psi](A_{\phi})$, we have chosen

[TABLE]

where the values of constant coefficient $a_{\Psi}$ and index $p$ are given by hand, for which we set $p=1$ following YYE06 .

III.5 Alternative form of the first integral of

MHD-Euler equations under pure rotational flows

When the flow fields are pure rotational $u^{A}=0$, the 4-velocity $u^{\alpha}$ becomes

[TABLE]

For the stationary and axisymmetric perfect fluid spacetime without magnetic fields, the first integral of the relativistic Euler equations can be derived as a consequence of the Cartan identity (2). For the simplest uniformly rotating case, $\Omega={\rm constant}$, the helical vector $k=t+\Omega\phi$ becomes a Killing vector. Then, the Euler equations are written, with the help of Eq. (2),

[TABLE]

and hence the first integral is derived as $h{\underline{u}}\cdot k=\mbox{constant}$. This relation is used for determining a thermodynamic variable, the enthalpy $h$ in this case, of uniformly rotating non-magnetized stars.

In the presence of magnetic fields, the corresponding first integral (130) (or (265)) for determining the relativistic enthalpy $h$, in place of the above relation $h{\underline{u}}\cdot k=\mbox{constant}$, was not derived from the Cartan identity (2) as discussed in previous sections. We show an alternative derivation of the first integral of ideal MHD flow using the Cartan identity (2), which is valid only for the case of pure rotational flows (170).

The canonical momentum, $\pi_{\alpha}=hu_{\alpha}$, respects the symmetries $\mbox{\pounds}_{t}\pi_{\alpha}=\mbox{\pounds}_{\phi}\pi_{\alpha}=0$. Although the angular velocity $\Omega$ is a certain function which also respects the symmetries $\mbox{\pounds}_{t}\Omega=\mbox{\pounds}_{\phi}\Omega=0$, $\Omega\phi^{\alpha}$ is not a Killing vector. Hence for a certain p-form $Q$, a relation,

[TABLE]

is satisfied. Then, the first term of the MHD-Euler equations (43) divided by the enthalpy $h$ becomes

[TABLE]

Hence, Eq. (43) is rewritten,

[TABLE]

Substituting the current (124), the $t$ and $\phi$-components of Eq. (174) become,

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

Above we used $t\cdot dQ=0$ and $\phi\cdot dQ=0$ for a scalar $Q$. For the $x^{A}$-components, combining Eq. (174) dotted with the basis $e_{A}$

[TABLE]

and

[TABLE]

we have

$x^{A}$-component:

[TABLE]

where $e_{A}\cdot dQ=\partial_{A}Q$ for a scalar $Q$.

Substituting the $x^{A}$-components of Maxwell’s equations (122), to the $t$ and $\phi$ components of MHD-Euler equations (175) and (176),

[TABLE]

the integrability conditions,

[TABLE]

are derived. Hence, using Eqs. (122) and (97) the current term of $x^{A}$-components (179) becomes

[TABLE]

Substituting Eq. (183), and the integrability conditions (182) for $A_{t}$ and $A_{\phi}$ and (115) for $\Omega$ into Eq. (179), we have

$x^{A}$-component:

[TABLE]

To derive a first integral of this Eq. (184), we have a few choices to reduce the first two terms: we may assume any of (i) $s=\mbox{constant}$, (ii) $T/h=[T/h](s)$, (iii) $\rho h=[\rho h](p)$ with $dh-Tds=dp/\rho$, or (iv) $s=s(\Upsilon)$. Then the following argument goes analogously. Here we assume a homentropic fluid, $s=\mbox{constant}$, for simplicity, and rewrite Eq. (184),

[TABLE]

or we separate the contribution of differential rotation by defining $j(\Omega):=u^{t}u_{\phi}$, and rewrite

[TABLE]

From the converse of the Poincaré lemma, the integrability condition becomes

[TABLE]

For the latter Eq. (186), the arbitrary function $\lambda(\Upsilon)$ is related to the current as

[TABLE]

and the first integral is written

[TABLE]

where ${\cal E}$ is a constant.

We have implemented this formulation, assuming $\Upsilon=A_{\phi}$, and $\Omega={\rm constant}$,

[TABLE]

and, although we do not show the result in this paper, we have computed a few solutions which agree well with those calculated from Eq. (151).

III.6 Remarks on numerical method

III.6.1 Finite difference and iteration

Given the forms of functions presented in Sec III.4, a set of integral equations of the field equations (139) and algebraic equations arranged from the first integrals and integrability conditions (148)-(151) derived in Sec. III.2 and III.3 are a full system of equations for computing magnetized rotating compact stars. These equations are discretized on spherical coordinates $(r,\theta,\phi)\in[r_{a},r_{b}]\times[0,\pi]\times[0,2\pi]$ that cover a star and an asymptotic region, where the origin $r_{a}=0$ is located at the center of the star. Then a self-consistent field iteration method is applied to calculate a converged solution. The numerical code is developed on cocal code extending previously developed rotating compact star code in which the waveless formulation is used Shibata:2004qz ; Uryu:2016dqr ; cocal . A numerical method used in the present code for magnetized compact stars, including setups for coordinate grid points and grid spacing, finite difference schemes for derivatives and integrals, and the self-consistent iteration scheme, are common with the above mentioned rotating compact star code on cocal. Readers who are interested in the details of the code are advised to consult papers Uryu:2016dqr ; cocal .

In Table 1, we reproduce a list of relevant parameters for the coordinate grids presented in previous papers Uryu:2016dqr ; cocal , and, in Table 2, present the numbers of grid points and other grid parameters used in actual computations shown in the later sections. An important difference from the previous calculations for non-magnetized rotating stars is the inclusion of higher multipoles and higher resolution in the $\theta$ direction. As we will see below, stronger toroidal magnetic fields tend to concentrate near the equatorial plane, hence it is necessary to increase the number of terms in the multipole expansion in (140) and (142) to as high as $L\gtrsim 30$, and accordingly the grid points in the $\theta$ direction to $N_{\theta}\gtrsim 144$.

Typically, with the grid setup SE3 $L=40$ in Table 2, it takes 6 minutes per 1 iteration using 1 CPU thread of Xeon E5-2687W v3 3.1GHz, and for a convergence about 500-1000 iterations are required.

III.6.2 Parameters

In our formulation, parameters to specify a magnetized rotating model appear in the integrability conditions shown in Sec. III.4.2. For the case without meridional flows, those are $\Lambda_{0}$, $\Lambda_{1}$, and $\cal E$ in Eq. (157), $\Omega_{\rm c}$ and $C_{e}$ in Eq. (158), and $\Lambda_{\phi 0}$ in Eq. (159). A set of parameters $b$ and $c$ contained in smooth step functions $\Xi(A_{\phi})$ in Eqs. (157) and (159) may be distinct in general but are set to have the same value in both equations. In addition to these parameters, we augment the number of parameters by introducing an equatorial radius $R_{0}$ in coordinate length for rescaling the radial coordinate $r$ Uryu:1999uu .

Another set of parameters is introduced from the EOS, which is also one of the integrability conditions. In cocal, a piecewise polytropic EOS and a variant of such piecewise EOS to model, for example, quark matter, are implemented CocalQuark . In this paper, we simply use a (single segment) polytropic EOS

[TABLE]

which introduces two parameters, $K$ and $\Gamma$, the adiabatic constant and the adiabatic index.

The values of the parameters $\{\Lambda_{0},\Lambda_{1},\Lambda_{\phi 0},b,c\}$ are prescribed, which control the strength of electromagnetic fields. As in the computations of non-magnetized rotating compact stars, three parameters, $\{{\cal E},\Omega_{\rm c},R_{0}\}$, are determined from three conditions, which are a given value of the maximum rest mass density $\rho_{\rm c}$, at (or near) the center of the star, the normalization of the equatorial radius, $r_{\rm eq}$, as $r_{\rm eq}/R_{0}=1$, and the given value for the deformation $r_{p}/r_{\rm eq}$ at the north pole, $r_{p}$. These conditions are imposed on Eq. (151) and the resulting set of three algebraic equations is simultaneously solved to determine $\{{\cal E},\Omega_{\rm c},R_{0}\}$ during iteration.

Finally, the parameter $C_{e}$ in Eq. (158) is left to be determined. We fix this value from the condition that the asymptotic (net) electric charge $Q$ vanishes:

[TABLE]

where $\int_{\infty}$ is the surface integral over a sphere $S_{r}$ with radius $r$, in the limit, $\int_{\infty}:={\displaystyle\lim_{r\rightarrow\infty}}\int_{S_{r}}$. This integral is evaluated at a large radius of our computational region, typically $r\sim 10^{4}R_{0}$, at every 30 iterations, then the secant method is applied to find a solution of $C_{e}$ to have $Q=0$. One can also compute a charged solution with setting a finite value to $Q$.

IV Results

In Uryu:2014tda , we have presented preliminary results for relativistic rotating star solutions associated with mixed poloidal and toroidal magnetic fields. As mentioned in Sec. III.4, we have modified the form of arbitrary functions from those used in Uryu:2014tda . We have also improved numerical codes to maintain expected accuracy; for example, the virial relation is satisfied in higher precision. The numerical computations presented below are performed using smaller to larger numbers of grid points and multipoles as shown in Tables 1 and 2 for studying the convergence of the solutions.

In the following computations, we choose the adiabatic EOS (191) with indices $\Gamma=2$ or $3$ and the adiabatic constant $K$ so that the compactness of a spherical solution having rest mass $M_{0}=1.5M_{\odot}$ becomes $M/R=0.2$. Reference quantities for the Tolman-Oppenheimer-Volkov (TOV) solutions for these EOSs are tabulated in Table 3. For the model parameters to determine the strength of electromagnetic fields, we choose three sets listed in Table 4. To our knowledge, since our first paper Uryu:2014tda was published, cocal is the only code that can calculate fully relativistic rotating compact stars associated with mixed poloidal and toroidal magnetic fields without any approximation in the formulation other than assumptions of stationarity and axisymmetry.

IV.1 Extremely magnetized solutions

We present three solutions of magnetized rotating compact stars in Fig. 1, and corresponding physical quantities in Tables 5 and 6. Definitions of these quantities are summarized in Appendix F. The model parameter of each solution is P1, P2, and P3, respectively in Table 4, where the model P1 is a normal mass solution with $\Gamma=2$ EOS, P2 is a supramassive solution with $\Gamma=2$ and is rotaing near the Kepler limit, and P3 is a normal mass solution with $\Gamma=3$.

As shown in Table 6, these solutions are associated with extremely strong poloidal and toroidal magnetic fields about an order of $\sim 10^{17}$–$10^{18}$[G], while the mass and radius of these compact stars are close to those of common neutron stars. For the models P1 and P3, the maximum values of the toroidal and poloidal components, ${B^{\rm max}_{\rm pol}}$ and ${B^{\rm max}_{\rm tor}}$ respectively, are comparable, and even for P2, ${B^{\rm max}_{\rm tor}}$ is about $1/3$ of ${B^{\rm max}_{\rm pol}}$. As reported also in other works, however, the bulk energy of the toroidal magnetic fields ${\cal M}_{\rm tor}$ is much smaller than that of the poloidal fields ${\cal M}_{\rm pol}$; as shown in Table 6, the energy of the poloidal fields accounts for more than 90% of the total electromagnetic energy $\cal M$.

In the top to bottom left panels of Fig. 1, the contours of $p/\rho$ and the poloidal and toroidal magnetic fields are presented. Although the toroidal magnetic field component ${B_{\rm tor}}$ is not dominating in the whole electromagnetic energy, ${B_{\rm tor}}$ is concentrated near the equatorial surface so that its maximum value is comparable to that of poloidal component ${B_{\rm pol}}$. This feature has been often observed in the other Newtonian MRNSNewtonian or approximate calculations PolTorMS .

A new feature can be seen in these panels for models P1 and P2. When the toroidal field ${B_{\rm tor}}$ is extremely strong, the magnetic energy density locally dominates over the mass energy density, and hence expel the matter from the region of extremely strong toroidal magnetic fields. In the middle left panel for the model P2, we can observe that the $p/\rho$ contours are deformed around the ${B^{\rm max}_{\rm tor}}$, and in the top left panel for the model P1, there are small closed circles of the density contours near the equatorial surface. For the model P1, a profile of $p/\rho$ along the equatorial radius near the surface (and hence $\rho$ or $\varepsilon$) almost drops to zero. Hence, we expect that, with a little stronger magnetic fields, which can be easily achieved by changing the parameters in Table 4, the matter will be completely expelled from this region, and hence a toroidal electro-vacuum tunnel will be formed inside of the compact star.

Roughly speaking, this happens because the pressure/energy density of the electromagnetic fields dominates over those of the matter in this toroidal region near the surface. To see this, in Fig. 2, we show the plots of spatial trace part of stress-energy tensor $T_{a}{}^{a}=T^{\alpha\beta}\gamma_{{\alpha\beta}}=(T_{M}^{\alpha\beta}+T_{F}^{\alpha\beta})\gamma_{{\alpha\beta}}$ separating contributions from the matter $T_{M}^{\alpha\beta}$ (14) ($T_{a}{}^{a}=T_{M}^{\alpha\beta}\gamma_{{\alpha\beta}}$) and the electromagnetic fields $T_{F}^{\alpha\beta}$ (15) ($T_{a}{}^{a}=T_{F}^{\alpha\beta}\gamma_{{\alpha\beta}}$). As can be seen in the left panel for the model P1, the dominance of the matter to the electromagnetic fields exchange in this region. In the right panel for model P2, there is a sizable amount of contribution from $T_{F}^{\alpha\beta}$ but it does not dominate over $T_{M}^{\alpha\beta}$.

In the middle panel of each row of Fig. 1, contours of metric potentials around the compact stars are plotted. Using the waveless formulation, we are able to compute non-conformal flat components of the metric such as $h_{xz}$ as shown in these panels.

In the right panel of each row of Fig. 1, contours of $t$ and $\phi$ components of the electromagnetic 1-form $A_{\alpha}$ are shown. As mentioned in Sec. III.4, integrability of ideal MHD equations require $A_{t}$ to be a function of a master potential for which can be seen correctly imposed on the fluid support of compact stars. Since we assume electro-vacuum spacetime outside of the star, the $A_{t}$ component is continuously but not smoothly connected at the stellar surface. Since we also assume that the net charge at infinity $Q$ (evaluated at a larger radius from the source in actual computations) vanishes, the contours of $A_{t}$ become positive (red dashed curves) near the poles, and negative (blue dashed curves) near the equator. These panels show that the method of solving for $A_{t}$ as described in Sec. III.2 is working consistently in these computations.

IV.2 Convergence test

In Fig. 3, the convergence of integrated quantities are plotted for the models P1, P2, and P3. Those are the convergence of $I_{\rm vir}/|{\cal W}|$ in the left, and that of $|1-M_{\rm K}/M_{\rm ADM}|$ in the right panels. The relativistic virial relation $I_{\rm vir}$ is defined as the volume integral of the spatial trace of Einstein’s equations as (283) in the Appendix, and its residuals shown in the left panel decrease as ${O}(\Delta r^{2})$ as expected. Strictly speaking, the numerically evaluated volume integral $I_{\rm vir}$ does not approach to zero as the resolution goes much higher, because it is evaluated on a large but finite computational domain $r\in[0,10^{6}R_{0}]$, and also evaluated from a solution in which a large but a finite number of multipoles are used for approximation. Hence what we can conclude here is the fact that the actual value of $I_{\rm vir}/|{\cal W}|$ in our setup is smaller than the finite difference error and hence can not be probed with the present highest resolutions such as SD3 or SE3.

The asymptotic Komar and ADM mass, $M_{\rm K}$ and $M_{\rm ADM}$ are known to agree in the framework of the waveless formulation under the gauge choice Eq. (26) Shibata:2004qz . The residual $|1-M_{\rm K}/M_{\rm ADM}|$ of each model decrease slower than ${O}(\Delta r^{2})$ as shown in the right panel, which is the same behavior as that of non-magnetized rotating compact stars Uryu:2016dqr . Hence we may conclude that the strongly magnetized solutions are calculated with comparable precision as the non-magnetized solutions in the cocal code.

In Fig. 4, we show the profile of $p/\rho$ along the x-axis (the radial coordinate in the equatorial plane) near the surface for the model P1. As mentioned in the above, extremely strong toroidal magnetic fields expel the matter from the toroidal region. To resolve such a relatively small scale toroidal structure, it is necessary to increase the numbers of grid points and multipoles. We performed convergence tests to examine the profile of $p/\rho$ with systematically increasing the resolution under a fixed number of multipoles, and also increasing the number of multipoles under a fixed resolution. In the top panel of Fig. 4, the number of multipoles is fixed to $L=40$, and the resolution is increased from SE12 to SE3 in Table 2. It can be observed that the largest error appears at the bottom of the $p/\rho$ profile around $x=0.89R_{0}$, and that the profile converges at the levels of resolutions around SE23 to SE3.

In the bottom panel of Fig. 4, plotted are a convergence of $p/\rho$ profiles near the equatorial surface of the model P1 with respect to the number of multipoles used in the elliptic solver (139). In this test, resolution SE3 and modified resolutions of it (SE3p, SE3t, SE3tp in Table 2) are used. It is confirmed that those modified resolutions do not affect the profile. For example, we compute the case with $L=50$ with SE3p and SE3t with increasing $N_{\theta}$, the number of grid points in $\theta$ coordinate, but those profiles overlap as seen in the plot. The profiles are also the same when the number of $N_{\phi}$ is increased, which is not related to an accuracy of a solution but is necessary for computating solutions with larger $L$. As seen in the panel, the profiles gradually change as increasing $L$ from 20 to 60. Because of the limit of computational resources, we do not perform computations higher than $L=60$ with a resolution SE3tp. The solutions still slightly changes from $L=50$ with SE3t resolution to the $L=60$ with SE3tp resolution, but the overall difference of the profile is getting smaller with increasing $L$.

V Discussion

In this paper, we have presented the full details of a formulation and a numerical method for computing stationary and axisymmetric equilibriums of fully relativistic rotating compact stars associated with mixed poloidal and toroidal magnetic fields. We have successfully calculated solutions associated with extremely strong poloidal and toroidal magnetic fields, and found solutions whose mass energy density is expelled by the energy density of toroidal magnetic fields. As presented in Fig. 4, the structure of this toroidal low density region can be calculated accurately using a large number of multipoles in our Poisson solver (Sec. III.2). From the top left panel of Fig. 1, the size of the toroidal region in the $\theta$ direction is around $\sim 0.03$ radian and hence requires a resolution $N_{\theta}>\pi/0.03\sim 100$. The number $N_{\theta}=384$ of the resolution SE3tp is sufficient to resolve this structure. On the other hand, a Legendre polynomial $P_{60}^{0}(\cos\theta)$ has only 60 nodes, which may resolve roughly $\pi/60\sim 0.05$ radians in the $\theta$ direction, and hence this is a reason for slow convergence in the number of multipoles $L$.

One of the new features of our method is to solve all components of Maxwell’s equations to determine all components of the electromagnetic potential 1-form $A_{\alpha}$. This allows us to compute electromagnetic configurations under various circumstances. In this paper, we assume an electro-vacuum spacetime outside of the compact star, which results in a surface charge distribution when we compute an asymptotically charge neutral solution. This is the same assumption used in the first relativistic computation of magnetized rotating equilibrium by the Meudon group Bocquet:1995je . It would be also possible to compute a solution which continues to force free magnetosphere outside of the star. The cocal code is also available for computing black holes. Hence an extension to a black hole associated with electromagnetic fields, that is, a black hole magnetosphere, is also a part of our future project.

Acknowledgements.

This work was supported by JSPS Grant-in-Aid for Scientific Research(C) 18K03624, 15K05085, 18K03606, 17K05447, NSF Grant PHY-1662211, NASA Grant 80NSSC17K0070, and the Marie Sklodowska-Curie grant agreement No.753115.

Appendix A Notation and relations

A.1 Orthogonal curvilinear coordinates and basis

When a system of coordinate functions $x^{A}$ is introduced for a set of points $x^{\alpha}$ in a space, the basis and dual basis are respectively $\{e_{A}^{\alpha}\}$ and $\{\nabla_{\alpha}x^{A}\}$, where $e_{A}^{\alpha}:=\frac{\partial x^{\alpha}}{\partial x^{A}}$. The coordinates are called orthogonal when

[TABLE]

for any pair of indices $A$ and $B$. Derivatives of this expression give

[TABLE]

since $\nabla_{\alpha}\delta_{A}{}^{B}=0$. Projecting these to the basis $e_{B}^{\alpha}$ we have

[TABLE]

for any dual basis $\nabla_{\alpha}x^{C}$. Therefore, we have

[TABLE]

A.2 Index free notation for differential forms and vectors

In some manipulations of equations in Sec. II and III, we use index free notations for differential forms and vectors for convenience. We summarize correspondences between index and index free notations in this subsection. For more details see e.g. Gourgoulhon:2011gz ; Gourg2012

A.2.1 Differential forms

For a 1-form $w$, we write an exterior derivative $dw$ which corresponds to the index notation as

[TABLE]

This notation is often used for an electromagnetic potential 1-form $A$, and for a canonical momentum 1-form $h{\underline{u}}$, where $\underline{u}$ is a dual 1-form of the 4-velocity vector $u$ (which is $u^{\alpha}$ in the abstract index notation). Faraday 2-form $F=dA$ is written,

[TABLE]

$F$ is a closed 2-form, $dF=0$, which is written in index notation,

[TABLE]

A.2.2 Inner product and Cartan identity

Inner product of a $p$-form $\omega$ and a vector $u$ is denoted with a dot in index free notation,

[TABLE]

Using this, the Cartan identity for a $p$-form $\omega$ is written

[TABLE]

in index free notation. As a rule for the inner product between a vector and a $p$-form, we assume that when the vector is operated to $p$-form from left (right), the vector is contracted with the left (right) most index of the $p$-form, for example, $u\cdot F=-F\cdot u$ for a 2-form $F$.

A.2.3 Ideal MHD condition

Ideal MHD condition is written $F\cdot u=0$. This implies $\mbox{\pounds}_{u}F=0$. This is shown using $dF=0$ as

[TABLE]

Or, using a potential 1-form $A$, ($F=dA$)

[TABLE]

where $u\cdot dA=u\cdot F=0$, $d^{2}=0$ and Cartan identity are used.

Also for any vector proportional to $u$, that is, with an arbitrary scalar function $\lambda$, $F\cdot(\lambda u)=0$ holds , which implies $\mbox{\pounds}_{\lambda u}F=0$. This guarantees that a flux of $F$ over any surface along a given family of flow lines is conserved Carter79 .

A.2.4 Integrability condition

When smooth functions $A$ and $f$ satisfy $fdA=C={\rm const}$, a relation is derived:

[TABLE]

Hence $f=f(A)$. As $f(A)dA=dF(A)=C={\rm const}$, $dF(A)/dA=f$. If the constant $C=0$, then $f=dF/dA=0$.

A.3 4-velocity

We decompose the 4-velocity $u^{\alpha}$ with respect to $t^{\alpha}$ as

[TABLE]

where $v^{\alpha}$ is spatial vector $v^{\alpha}\nabla_{\alpha}t=0$. In the first integrals and in the currents, $t$ and $\phi$ components $u_{\alpha}$ appears, which are calculated as follows;

[TABLE]

[TABLE]

where $t^{\alpha}=\alpha n^{\alpha}+\beta^{\alpha}$ is used.

The coordinate basis for vector $\phi^{\alpha}$ is related to the Cartesian basis $\hat{x}^{\alpha}$ and $\hat{y}^{\alpha}$ as

[TABLE]

and the basis for 1-form

[TABLE]

Hence, the $\phi$-components of the 4-velocity are related to those of Cartesian coordinates as

[TABLE]

and these becomes on the $\phi=0$ (meridional) plane,

[TABLE]

The same relations between $\phi$-components and the Cartesian components are used for the electric currents $j^{\alpha}$ ($j^{\phi}$ and $j_{\phi}$).

Appendix B 3+1 decomposition of Faraday tensor

In this Appendix, we derive the 3+1 form of Faraday tensor and its divergence. The spatial projection of $F_{\alpha\beta}=(dA)_{\alpha\beta}$ can be derived explicitly as follows; For $\bar{F}_{\alpha}$, we use the Cartan identity $n\cdot dA=\mbox{\pounds}_{n}A-d(n\cdot A)$,

[TABLE]

where the relation, $\mbox{\pounds}_{n}n_{\alpha}=D_{\alpha}\ln\alpha$, is used. For $\bar{F}_{\alpha\beta}$, since $F_{\alpha\beta}=(dA)_{\alpha\beta}$ is independent of the geometry of the manifold, its spatial projection becomes its spatial part,

[TABLE]

which can be shown more explicitly as

[TABLE]

as $K_{\alpha\beta}$ is a symmetric tensor. The divergence $\nabla_{\beta}F^{\alpha\beta}$ is also decomposed with respect to $\Sigma_{t}$. The projection of $\nabla_{\beta}F^{\alpha\beta}$ to the hypersurface normal $n_{\alpha}$ becomes

[TABLE]

The projection of $\nabla_{\beta}F^{\alpha\beta}$ to the hypersurface $\Sigma_{t}$ becomes

[TABLE]

Hence, on $\Sigma_{t}$, we have

[TABLE]

Appendix C Derivation of equations for electromagnetic potentials

In this Appendix, we derive the final form of Maxwell’s equations implemented in the cocal code for computing electromagnetic potentials. Since we introduce a conformal decomposition of the spatial metric Eq. (6) as in Sec.II.2.1, the divergence with respect to the conformal metric $\tilde{\gamma}_{ab}$ is simplified to that of flat metric $f_{ab}$, \tilde{D}_{a}A^{a}={\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}A^{a}.

Projecting along $n^{\alpha}$, Eq. (65) becomes

[TABLE]

Separating the flat Laplacian from the first term,

[TABLE]

an elliptic equation for $\alpha\Phi_{\Sigma}$ is derived

[TABLE]

where the source $S$ is written

[TABLE]

The second term of the source (226) vanishes under the Dirac gauge condition {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}\tilde{\gamma}^{ab}=0. Also note that, the fourth and fifth terms are derived as below since $\tilde{\gamma}=f$ is satisfied,

[TABLE]

Projecting to $\Sigma_{t}$, Eq. (67) becomes

[TABLE]

The first term, from which an elliptic operator is separated as below, is rewritten,

[TABLE]

Using an identity,

[TABLE]

where ${\tilde{v}}^{a}=\tilde{\gamma}^{ab}v_{b}$, we have

[TABLE]

Hence,

[TABLE]

From the first term of the right hand side of Eq. (232), the flat Laplacian \displaystyle-{\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\Delta\tilde{A}_{a} is isolated,

[TABLE]

We keep $\tilde{D}_{a}$ instead of replacing it by {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a} and a connection $C^{c}_{ab}$ in a couple of terms in the Eq. (233), to shorten the equation. Then, a set of elliptic equations for $A_{a}$ is derived,

[TABLE]

where the source $S_{a}$ is written

[TABLE]

Appendix D Derivation of first integrals of MHD-Euler equations

In this Appendix, we derive a set of integrability conditions and first integrals (240)-(246) of the relativistic MHD-Euler equations (43).

For the $t$ and $\phi$ components of the MHD-Euler equations (125) and (126), substituting Eq. (89) to $u^{A}$ and Eq. (122) to that of the current $j^{A}$ in the above set of equations, and multiplying by $\rho\sqrt{-g}$, we have

[TABLE]

Substituting the integrability conditions (110), these are rewritten

[TABLE]

These relations imply that the terms in parenthesis are a function of $\Upsilon$. Hence, introducing the densitized scalars $[\sqrt{-g}\Lambda_{t}](\Upsilon)$, and $[\sqrt{-g}\Lambda_{\phi}](\Upsilon)$, for each component, the sufficient conditions for the $t$ and $\phi$ components are written,

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

These are combined and written using another function of $\Upsilon$, $\Lambda(\Upsilon)$,

[TABLE]

For the $x^{A}$-component (127), multiplying $\rho\sqrt{-g}$, and, substituting Eq. (89) and Eq. (122) as well as definitions $(dA)_{AB}=\epsilon_{AB}B_{\phi}$ and $d(h{\underline{u}})_{AB}=-\epsilon_{AB}\omega_{\phi}$, we have

[TABLE]

In this Eq. (243), the first two terms and the sixth term multiplied by $A^{\prime}_{t}A^{\prime}_{\phi}$ become as follows; substituting the integrals of $t$ and $\phi$ components of MHD-Euler equations (240) and (241),

[TABLE]

The terms in the first parenthesis of r.h.s. vanish because of the first integral of the $x^{A}$-components of the ideal MHD condition (114), and all other terms are proportional to $\partial_{A}\Upsilon$, as $A_{t}$, $A_{\phi}$, $\sqrt{-g}\Psi$, $\sqrt{-g}\Lambda_{t}$, $\sqrt{-g}\Lambda_{\phi}$, and $\Lambda$, are functions of $\Upsilon$. Hence, with an assumption to the thermodynamic variable, that is, the entropy $s$ to be a function of $\Upsilon$, $s=s(\Upsilon)$, Eq. (243) multiplied by $A^{\prime}_{t}A^{\prime}_{\phi}$ is rewritten

[TABLE]

Therefore, we obtain the first integral of the $x^{A}$-components of the MHD-Euler equations (243) as

[TABLE]

Appendix E First integral and integrability conditions for the case

of pure rotational flow

For the case without meridional flow fields, $u^{A}=0$, the system of first integrals and Maxwell’s equations under stationarity and axisymmetry can be recast into a single equation to be solved for a single independent variable which may be called the Grad-Shafranov equation with a toroidal flow fields. However as we have noted, we do not reduce the number of variables in our formulation, but rather solve the hydrostationary equation and Maxwell’s equations simultaneously. Although the derivation presented in Sec. II.6 can be applied for the case with pure rotational flow, we repeat the derivation below for clarity.

E.1 Ideal MHD condition for purely rotational flow

We assume the 4-velocity $u^{\alpha}$ of the flow field in the absence of meridional flow $u^{A}=0$ as

[TABLE]

The ideal MHD condition $F_{\alpha\beta}u^{\beta}=0$ in this case becomes

$t$-component

[TABLE]

$\phi$-component

[TABLE]

$x^{A}$-component

[TABLE]

For the case with meridional flow, integrability conditions can be found in the above ideal MHD condition alone. The absence of the meridional stream function in this case trivializes the $t$ and $\phi$ components of ideal MHD conditions and hence the integrability conditions are not derived from these equations.

E.2 MHD-Euler equations for pure rotational flow

Substituting $u^{A}=0$ and $j^{\alpha}=j^{t}t^{\alpha}+j^{\phi}\phi^{\alpha}+j^{A}e_{A}^{\alpha}$, the MHD-Euler equations becomes

$t$-component:

[TABLE]

$\phi$-component:

[TABLE]

$x^{A}$-component:

[TABLE]

E.3 Integrability conditions for the case of purely rotational flow

Substituting the $x^{A}$-components of Maxwell’s equations (122) to the meridional current $j^{A}$ appearing in the $t$ and $\phi$ components of the MHD-Euler equations (251) and (252), we have

[TABLE]

These relations require integrability conditions for consistency, namely, a master potential $\Upsilon$ is introduced as follows,

[TABLE]

The $x^{A}$-components of ideal MHD conditions (250), and the above integrability conditions for $A_{t}$ and $A_{\phi}$ imply

[TABLE]

and hence,

[TABLE]

or, introducing the angular velocity $\Omega$,

[TABLE]

Therefore, $\Omega$ should be a function of $\Upsilon$ as well,

[TABLE]

E.4 First integral of meridional components of MHD-Euler equations

Derivation of the integrability of $x^{A}$-component of MHD-Euler equations proceeds analogous to the case with non-zero meridional flow. A difference is an absence of the stream function $\sqrt{-g}\Psi$. In Eq. (130), the stream function $[\sqrt{-g}\Psi]^{\prime}(\Upsilon)$ appears in the denominator of the definition of an arbitrary function $[\sqrt{-g}\Lambda](\Upsilon)$. In paper Gourgoulhon:2011gz , we have proved that such a combination always become finite, and hence the relations derived in previous sections are valid also for the case of pure rotational flow under simply taking a limit $[\sqrt{-g}\Psi](\Upsilon)\rightarrow{\rm constant}$. In this section, we prove this fact by repeating the derivation of the previous section, and derive (130) directly as a part of integrability conditions.

We recast the $x^{A}$-components of MHD-Euler equations in the same way as the case for generic flow in the previous section. To proceed, we use a relation,

[TABLE]

derived from normalization of the 4-velocity $u\cdot u=-1$, and the integrability condition (258). Multiplying by the factor $\frac{2A^{\prime}_{t}A^{\prime}_{\phi}}{A^{\prime}_{t}u^{t}-A^{\prime}_{\phi}u^{\phi}}$, the kinetic term of the $x^{A}$-components of MHD-Euler equations for purely rotational flow (253), is rewritten

[TABLE]

where a consistency Eq. (258) is used. Also, the Lorenz force term of Eq. (253) becomes

[TABLE]

with Eqs. (97) and (122).

Because $A_{t}$, $A_{\phi}$ and $B\sqrt{-g}$ are functions of the master potential $\Upsilon$ as shown in Eq. (256), the $x^{A}$-components of MHD-Euler equations (253) multiplied by the factor $\frac{2A^{\prime}_{t}A^{\prime}_{\phi}}{A^{\prime}_{t}u^{t}-A^{\prime}_{\phi}u^{\phi}}$, is rewritten, with an assumption of $s=s(\Upsilon)$,

[TABLE]

The above relation suggests that because of the converse of the Poincaré lemma, the first term is a function of $\Upsilon$,

[TABLE]

This is compared with Eq. (130).

Since $\sqrt{-g}B=[\sqrt{-g}B](\Upsilon)$, we introduce the following functions of $\Upsilon$,

[TABLE]

Taking derivatives of these with respect of $\Upsilon$, and combining them, we have a relation,

[TABLE]

Finally, substituting Eqs. (265) and (268) to (264), the consistency of the $x^{A}$ components yields

[TABLE]

This relation is compared with the result for the generic flow (131); Eq. (269) agrees with Eq. (131) in the limit $[\sqrt{-g}\Psi](\Upsilon)\rightarrow\mbox{constant}$.

Appendix F Definitions of mass, angular momentum, and virial relation

For the reader’s convenience, we summarize definitions of tabulated quantities in Tables 5 and 6 (see also e.g., RNS ; Gourg2012 ). Those include the rest mass $M_{0}$, ADM mass $M_{\rm ADM}$, Komar mass $M_{\rm K}$, total angular momentum $J$, the virial relation $I_{\rm vir}$ and other related quantities including electromagnetic energy $\cal M$ and its decomposition.

In Table 5, $\bar{R}_{0}$ and $\bar{R}_{z}$ are the equatorial and polar radius of a compact star in the proper length, respectively. The proper equatorial radus $\bar{R}_{0}$ is defined by

[TABLE]

and for $\bar{R}_{p}$, the integral is taken along $z$-axis. The unbarred $R_{0}$ and $R_{p}$ are the equatorial and polar radii in coordinate length.

The rest mass $M_{0}$ is written

[TABLE]

where $dS_{\alpha}=\nabla_{\alpha}t\sqrt{-g}d^{3}x$, and $d^{3}x=r^{2}\sin\theta drd\theta d\phi$ on the spherical coordinates. This is conserved irrespective of the choice of a slice $\Sigma$ for the rest mass conservation Eq. (42). The above volume integral over the hypersurface $\Sigma$ is non zero only on the fluid support.

The ADM mass $M_{\rm ADM}$ is defined and calculated by

[TABLE]

where $\rho_{H}:=T_{\alpha\beta}n^{\alpha}n^{\beta}$ is a component of the stress energy tensor $T_{\alpha\beta}$ normal to the hypersurface $\Sigma$. To check the consistency of the solution, both of the surface integral (273) and the volume integral (273) are evaluated. Relative errors between those values are $\sim 0.01\%$ for the solutions presented in Sec. IV. The surface integral is calculated on a sphere with radius around $r\sim 10^{4}R_{0}$, while the volume integral is taken within a computational domain within a radius around $r\sim 10^{6}R_{0}$. For the surface integral (273), we replace $\tilde{\gamma}^{ab}\rightarrow f^{ab}$, \tilde{D}_{a}\rightarrow{\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D_{a}, and $d\tilde{S}_{a}=\nabla_{a}r\sqrt{\tilde{\gamma}}d^{2}x=\nabla_{a}r\sqrt{f}d^{2}x=dS_{a}$. These are exact at spatial infinity, and they introduce only a negligible numerical error at the above radius where the surface integral (273) is evaluated.

The Komar mass $M_{\rm K}$ associated with the global timelike Killing field $t^{\alpha}$ is defined by

[TABLE]

and the asymptotic Komar mass whose $t^{\alpha}$ is a symmetry of an asymptotically flat spacetime by

[TABLE]

where the source terms $j_{a}:=-T_{\alpha\beta}\gamma^{\alpha}{}_{a}n^{\beta}$, and $S:=T_{\alpha\beta}\gamma^{\alpha\beta}$ are the components of the 3+1 decomposed stress-energy tensor $T_{\alpha\beta}$. In deriving (277), a relation, $(G_{\alpha\beta}-8\pi T_{\alpha\beta})(\gamma^{{\alpha\beta}}+n^{\alpha}n^{\beta})=0$ is used. Since $T_{\alpha\beta}=T^{M}_{\alpha\beta}+T^{F}_{\alpha\beta}$ contains electromagnetic contributions, the support of the volume integral (275) is non-compact. All integrals (275), (277), and (277) should reproduce the same value, when the waveless condition (27) and the coordinate conditions (26) are imposed at least asymptotically Shibata:2004qz . We computed Eqs. (275), (277), and (277) to check the consistency of the solutions, and found that they agree in the same order as $M_{\rm ADM}$ mentioned above.

For the total angular momentum $J$, the surface and volume integrals are evaluated,

[TABLE]

and the difference between the values from Eqs. (279) and from (279) is typically ${O}(0.1)\%$. Also for $J$, a term including $j_{a}$ in the volume integral contains contributions from fluid as well as electromagnetic fields, and hence it is integrated over a non-compact support. The values of $M_{\rm ADM}$, $J$, and $|1-M_{\rm K}/M_{\rm ADM}|$ listed in Table 5 are those of the volume integrals, (273), (275) and (279).111111 In previous papers for non-magnetized rotating stars Uryu:2016dqr ; RNScocal , the ratio of the kinetic energy and the gravitational energy, $T/|W|$ was defined following RNS ,

$\displaystyle W$ $\displaystyle:=$ $\displaystyle M_{\rm ADM}-M_{\rm P}-T,$

(280)

$\displaystyle T$ $\displaystyle:=$ $\displaystyle\textstyle\frac{1}{2}\int_{\Sigma}\Omega\,dJ.$

(281) where the proper mass $M_{\rm P}$ was defined by

$\textstyle M_{\rm P}\,:=\,\int_{\Sigma}T^{\alpha}\!{}_{\beta}u^{\beta}dS_{\alpha}\,=\,\int_{\Sigma}\epsilon u^{\alpha}dS_{\alpha}.$

(282) In the solutions presented in Sec. IV, $T^{\alpha}\!{}_{\beta}$ includes the electromagnetic $T^{F}_{\alpha\beta}$, while $u^{\alpha}$ is defined only on the fluid support. Likewise $\Omega$ in (281) is undefined outside of a star, although $dJ$ has a non zero electromagnetic contribution there. Because of this ambiguity, we do not calculate the values of $M_{\rm P}$ and $T/|W|$.

The relativistic virial theorem for an Einstein-Maxwell spacetime coupled with charged and magnetized perfect fluid Gourgoulhon & Bonazzola (1994) is computed to determine the accuracy of solutions. It is a vanishing integral of the spatial trace of Einstein’s equations over a hypersurface $\Sigma$,

[TABLE]

The equality of the ADM mass and the Komar mass $M_{\rm ADM}=M_{\rm K}$ has been proved for stationary spacetimes Beig (1978), and for the waveless formulation Shibata:2004qz . The integrals $\cal T$, $\Pi$, ${\cal M}$ and $\cal W$ are defined by

[TABLE]

which become the kinetic, internal, electromagnetic, and gravitational energies, in the Newtonian limit. In the integrand of $\cal W$ (F), ${}^{3}\!\tilde{R}$ is a scalar curvature of a conformally related spacelike hypersurface associated with a conformal 3-metric $\tilde{\gamma}_{ab}$. We define a virial integral $I_{\rm vir}$ as

[TABLE]

whose values for the selected solutions are presented in Table 6. The magnetic energy term $\cal M$ (286) is decomposed into contributions from the electric fields as well as the poloidal and toroidal magnetic fields, for which we define, respectively,

[TABLE]

which are also listed in Table 6.

Finally, the electric charge $Q$ defined in Eq. (192) becomes

[TABLE]

where $dS_{\alpha\beta}=\frac{1}{2}(\nabla_{\alpha}t\nabla_{\beta}r-\nabla_{\alpha}r\nabla_{\beta}t)\sqrt{-g}\,d^{2}x$, and $dS_{a}=\nabla_{a}r\sqrt{\gamma}\,d^{2}x$, which is evaluated on a large sphere $S$ in the asymptotics of $\Sigma$. Rewriting the charge $Q$ in the form of volume integral,

[TABLE]

the volume integral over the MHD fluid support $Q_{M}$ and the surface charge $Q_{S}$ at the stellar surface should contribute to the total charge $Q$. In our formulation, the form of $Q_{S}$ is not given, and the values of $Q_{M}$ are listed in Table 6.

Appendix G Imposition of symmetry of the electromagnetic

vector potential

When an exact 2-form $F=dA$ respects the symmetry $\mbox{\pounds}_{t}F=0$, a gauge potential $f$ exists such that $A$ transformed by $A\rightarrow A+df$ satisfies $\mbox{\pounds}_{t}(A+df)=0$.

proof) The Cartan identity $t\cdot dF=\mbox{\pounds}_{t}F-d(t\cdot F)$ implies $d(t\cdot F)=0$ when $\mbox{\pounds}_{t}F=0$. Hence, because of the Poincaré lemma, a function $\Phi$ exists such that $t\cdot F=d\Phi$ on a simply connected manifold. This implies

[TABLE]

Under a gauge transformation with a potential $f$, $A\rightarrow A+df$, $(F\rightarrow F)$, $\mbox{\pounds}_{t}A$ is transformed as

[TABLE]

Hence for an $f$ that satisfies

[TABLE]

$A+df$ satisfies the symmetry $\mbox{\pounds}_{t}(A+df)=0$.

We have freedom to choose another gauge potential ${\widehat{f}}$ that respects the symmetry $\mbox{\pounds}_{t}\widehat{f}=0$. This gauge potential does not affect the above transformation to impose the symmetry on the potential $A$, namely $A\rightarrow A+df+d\widehat{f}$ respects the symmetry. With this gauge freedom we may, for example, impose Coulomb gauge (vanishing spatial divergence) {\raise 4.30554pt\hbox{{}^{\ \circ}}}\!\!\!\!\!D^{a}A_{a}=0 where $a$ is a spatial 3D index.