Initial data for general relativistic simulations of multiple electrically charged black holes with linear and angular momenta

TL;DR

This paper develops a method to generate initial data for simulating charged binary black holes with linear and angular momenta in general relativity, advancing the study of such systems in strong gravitational fields.

Contribution

The authors introduce a new formalism within the conformal transverse-traceless framework to model charged black hole binaries with arbitrary momenta, extending existing numerical tools.

Findings

Code successfully recovers known solutions.

Demonstrates excellent self-convergence for generic configurations.

Enables future dynamical simulations of charged black hole systems.

Abstract

A general relativistic, stationary and axisymmetric black hole in a four-dimensional asymptotically-flat spacetime is fully determined by its mass, angular momentum and electric charge. The expectation that astrophysically relevant black holes do not posses charge has resulted in a limited number of investigations of moving and charged black holes in the dynamical, strong-field gravitational (and electromagnetic) regime, where numerical studies are necessary. Apart from having a theoretical interest, the advent of multimessenger astronomy with gravitational waves offers new ways to think about charged black holes. In this work, we initiate an exploration of charged binary black holes by generating valid initial data for general relativistic simulations of black hole systems that have generic electric charge, linear and angular momenta. We develop our initial data formalism within the framework of the conformal transverse-traceless (Bowen-York) technique using the puncture approach, and apply the theory of isolated horizons to attribute physical parameters (mass, charge and angular momentum) to each hole. We implemented our formalism in the case of a binary system by modifying the publicly available TwoPunctures and QuasiLocalMeasures codes. We demonstrate that our code can recover existing solutions and that it has excellent self-convergence properties for a generic configuration of two black holes.

Equations

$${}^{(4)}{\mathcal{R}}_{ab}-\frac{1}{2}{g}_{ab}{}^{(4)}\mathcal{R}$$

$${\nabla }_a{}^{(4)}{F}^{ab}$$

$${\nabla }_a{}^{(4)}{}^{\star }{F}^{ab}$$

$${}^{(4)}{}^{\star }{F}^{ab}=\frac{1}{2}{ϵ}^{abcd}{}^{(4)}{F}_{cd}.$$

$$4\pi {}^{(4)}{T}_{ab}^{\mathrm{EM}}={}^{(4)}{F}_{ac}{}^{(4)}{F}_{bd}{g}^{cd}-\frac{1}{4}{g}_{ab}{}^{(4)}{F}_{cd}{}^{(4)}{F}^{cd}.$$

$${\gamma }_b^a={\delta }_b^a+n^a{n}_b.$$

$${\gamma }_{ab}=g_{ab}+n_a{n}_b.$$

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

$$d{s}^2=-{\alpha }^{2}d{t}^2+{\gamma }_{ij}(d{x}^i+{\beta }^{i}dt)(d{x}^j+{\beta }^{j}dt).$$

$$K_{ab}=-{\gamma }_a^c{\gamma }_b^d{\nabla }_c{n}_d.$$

$$\mathcal{E}$$

$$S^i$$

$$\mathcal{R}+K^2-K_{ij}{K}^{ij}$$

$$D_j(K^{ij}-{\gamma }^{ij}K)$$

$$E^a$$

$$B^a$$

$${}^{(4)}{F}_{ab}=n_a{E}_b-n_b{E}_a+{ϵ}_{abcd}{B}^c{n}^d,$$

$${}^{(4)}{}^{\star }{F}_{ab}=n_a{B}_b-n_b{B}_a-{ϵ}_{abcd}{E}^c{n}^d.$$

$$D_a{E}^a$$

$$D_a{B}^a$$

$$\begin{gathered} 4\pi T_{\mathrm{EM}}^{ab}=\frac{1}{2}(n^a{n}^b+{\gamma }^{ab})(E_c{E}^c+B_c{B}^c) \\ +2n^{(a}{ϵ}^{b)cd}{E}_c{B}_d-(E^a{E}^b+B^a{B}^b), \end{gathered}$$

$$4\pi \mathcal{E}$$

$$4\pi S^i$$

$$\begin{gathered} d{s}^2=-\left(1-\frac{2m}{r}+\frac{q}{r^2}\right)d{t}^2+ \\ {\left(1-\frac{2m}{r}+\frac{q}{r^2}\right)}^{-1}d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2), \end{gathered}$$

$${}^{(4)}A=-\frac{q}{r}dt.$$

$$r=R\left(1+\frac{m}{R}+\frac{R_H^2}{R^2}\right),$$

$$d{s}^2=-{\Psi }^{-4}d{t}^2+{\Psi }^4{\delta }_{lk}d{x}^{l}d{x}^k,$$

$$\Psi =\sqrt{1+\frac{m}{R}+\frac{R_H^2}{R^2}}=\sqrt{{\left(1+\frac{m}{2R}\right)}^2-{\left(\frac{q}{2R}\right)}^2}.$$

$$E^R={\Psi }^{-6}\frac{q}{R^2}.$$

$$R_{\mathcal{S}}={\left(\frac{1}{4\pi }{\int }_{\mathcal{S}}ϵ\right)}^{\frac{1}{2}},$$

$$J_{\mathcal{S}}=-\frac{1}{8\pi }{\int }_{\mathcal{S}}(\phi \cdot \omega )ϵ+2(\phi \cdot {}^{(4)}A){}^{(4)}{}^{\star }F,$$

$$Q_{\mathcal{S}}=\frac{1}{4\pi }{\int }_{\mathcal{S}}{}^{(4)}{}^{\star }F,$$

$$M_{\mathcal{S}}=\frac{1}{2R_{\mathcal{S}}}{\left[(R_{\mathcal{S}}^2+Q_{\mathcal{S}}^2{)}^2+4J_{\mathcal{S}}^2\right]}^{\frac{1}{2}}.$$

$$J_{\mathcal{S}}^{\text{GR}}=-\frac{1}{8\pi }{\int }_{\mathcal{S}}(\phi \cdot \omega )ϵ=\frac{1}{8\pi }{\int }_{\mathcal{S}}{\phi }^a{R}^b{K}_{ab}ϵ,$$

Full text

I Introduction

The successful detection of gravitational waves from the inspiral and merger of binary black holes by the LIGO-Virgo interferometers Abbott et al. (2016a, b, 2017a, 2017b); The LIGO Scientific Collaboration and the Virgo Collaboration (2018) was made possible not only by technological advancements in instrumentation but also by substantial improvements in theoretical modeling that furnished the gravitational wave templates necessary for performing matched filtering Flanagan and Hughes (1998a, b); Aylott et al. (2009); Ajith et al. (2012); Hinder et al. (2014). To generate a bank of complete template signals, the equations of general relativity have to be solved during the late compact binary inspiral, merger and post-merger phases, because these events involve extreme gravitational fields, whose description with post-Newtonian methods is not accurate. Obtaining an analytic solution to describe these systems during these dynamic stages is not possible. Therefore, numerical integration of the full Einstein equations provides the only viable avenue for understanding such relativistic astrophysical systems from first principles and for helping to build gravitational wave templates during the most dynamical phases of their evolution.

Assuming that general relativity is the correct theory of gravity, the problem of two black holes is solved by integrating Einstein’s equations in vacuum. Despite the simpler description of black hole spacetimes compared to spacetimes with matter, it took decades for the field of numerical relativity to mature enough to be able to stably evolve two black holes until merger Pretorius (2005a); Campanelli et al. (2006); Baker et al. (2006). Some of the issues that hindered the development were due to the highly non-linear character of the Einstein equations, the coordinate freedom of general relativity, and the intrinsically singular nature of black holes. However, since the 2005 breakthrough, numerical relativity has advanced considerably with state-of-the-art codes that can simulate the inspiral and merger of uncharged binary black holes, and extract gravitational waves with high precision (see, e.g., Gopakumar et al. (2008); Lovelace et al. (2012); Chu et al. (2016); Lovelace et al. (2016); Jani et al. (2016); Healy et al. (2017, 2019) and references therein). Numerical relativity furnishes invaluable information for gravitational-wave detection and analysis, which includes the development of templates (see, e.g., Baumgarte et al. (2008); Scheel et al. (2009); Ajith et al. (2012); Hinder et al. (2014)), and the accurate parameter estimation of already detected events Lange et al. (2017).

Apart from binary black holes, binary neutron stars and binary black hole–neutron stars are also most the promising gravitational wave sources for currently operating interferometers Baumgarte and Shapiro (2010a). In fact, among the eleven confirmed detections of gravitational waves so far The LIGO Scientific Collaboration and the Virgo Collaboration (2018), event GW170817 is attributed to the inspiral and merger of a binary neutron star Abbott et al. (2017) (although a binary black hole–neutron star cannot be ruled out Yang et al. (2018); Hinderer et al. (2018); Foucart et al. (2018); Coughlin and Dietrich (2019) as a possibility). A complete simulation of compact binaries with matter requires the evolution of the spacetime coupled to matter, radiation/neutrinos, and electromagnetic fields in conjunction with detailed microphysics. A full solution including radiation/neutrinos without approximation is impossible at this time, and even with approximation, evolution of perfect fluids with existing numerical schemes involves density floors and other ad hoc prescriptions that are necessary to stabilize the calculations (see, e.g., Font (2008); Etienne et al. (2012a)), but are designed such that their impact on the global solution is minimal. However, this means that in a sense, simulations involving perfect fluids are not as “clean” as the ones in vacuum, which do not require ad hoc prescriptions. Nevertheless, many important results have been obtained through binary neutron star and binary black hole–neutron star simulations in full general relativity, see Shibata and Taniguchi (2011); Faber and Rasio (2012); Lehner and Pretorius (2014); Paschalidis (2017); Paschalidis and Stergioulas (2017); Baiotti and Rezzolla (2017); Duez and Zlochower (2019) for reviews (see also Cardoso et al. (2012) for other applications of numerical relativity).

Interesting spacetimes that are as “clean” as vacuum spacetimes, but have received little attention in numerical relativity, are those described by Einstein-Maxwell’s theory. This theory involves only gravitational and electromagnetic fields, and the corresponding spacetimes are referred to as electrovacuums or electrovacs. However, force-free electrodynamics has received some attention Palenzuela et al. (2010); Lehner et al. (2012); Moesta et al. (2012); Alic et al. (2012); Paschalidis and Shapiro (2013); Paschalidis et al. (2013a); Ponce et al. (2014); East and Yang (2018), but those simulations are not “clean”, in the sense that when the force-free conditions are violated during the evolution (typically in current sheets), one must interfere and enforce them to continue the calculations. On the other hand, electrovacuum spacetimes can be solved without physical approximations or ad hoc prescriptions, as the only assumption here is that electromagnetism and gravitation are described by the source-free Einstein-Maxwell equations. This simplification is the reason why these spacetimes have attracted numerous theoretical and analytic investigations for a long time, including the celebrated Kaluza-Klein theory Kaluza (1921); Klein (1926) unifying gravity and electromagnetism.

Examples of interesting electrovacuum spacetimes are those with electrically charged black holes.111It is also possible to include magnetic charges. This will not be done in the study presented in this paper, so we always take the term charge to mean electric charge. We note the extension of the work to include magnetic charges would be straightforward. The case of a single charged non-rotating black hole is analytically solved by the Reissner-Nordström metric Reissner (1916); Nordström (1918). This solution has been extended to non-vanishing angular momentum in the Kerr-Newman spacetime Newman et al. (1965), which generalizes the uncharged rotating black hole solution found by Kerr Kerr (1963). Another interesting class of solutions with multiple black holes is the static Majumdar-Papapetrou solution Majumdar (1947); Papapetrou (1945) that describes non-spinning black holes whose electric repulsion and gravitational attraction balance, producing a zero net force condition, and thus equilibrium. The hypothesis of staticity was relaxed to simple stationarity by Perjés (1971); Israel and Wilson (1972); Hartle and Hawking (1972). This list summarizes the known analytical solutions of the source-free Einstein-Maxwell equations in four-dimensional asymptotically-flat spacetimes.

A reason why the source-free Einstein-Maxwell theory has been primarily confined to the realm of theoretical explorations is the fact that astrophysically relevant black holes are not believed to be electrically charged, as the charge would be neutralized by the surrounding plasma Wald (1984) or as result of a pair-production through a Schwinger-like process Gibbons (1975). Nonetheless, there are some viable mechanisms that result in a black hole with non-zero charge. One example is the model proposed by Wald (1974), where the charge is retained due to the presence of an external magnetic field. This is known as the “Wald mechanism”. It was shown in Wald (1974) that if an asymptotically uniform magnetic field $B_{0}$ can be sustained, a black hole with mass $M$ spinning with angular momentum $J$ would acquire an electric charge $Q=2B_{0}J$ (measured in geometrized units222The conversion factor between our units and the International System of Units (SI) is $c^{2}G^{-\frac{1}{2}}(4\pi\varepsilon_{0})^{\frac{1}{2}}=$1.16\text{\times}{10}^{20}\text{,}\mathrm{C}\text{,}{\mathrm{km}}^{-1}$, so$1\text{,}\mathrm{\text{$M_{\odot}$}}$=$1.71\text{\times}{10}^{20}\text{,}\mathrm{C}, with $c$ speed of light in vacuum, $G$ gravitational constant, $\varepsilon_{0}$ vacuum permittivity and $M_{\odot}$ solar mass.), which we can rewrite as $Q/M=2B_{0}\chi M$ with $\chi=J/M^{2}$ the black hole dimensionless spin parameter. Since for black holes $\chi^{2}\leq 1$, there exists a maximum possible charge-to-mass ratio in the Wald mechanism: $(Q/M)\leq(Q/M)_{\rm max}\equiv 2B_{0}M$ Wald ([1974](#bib.bib61)). In the case of a solar mass black hole in the galactic magnetic field Eckart *et al.* ([2012](#bib.bib62)); Eatough *et al.* ([2013](#bib.bib63)) the ratio has to be $Q/M\leq${10}^{-24}. The charge-to-mass ratio quantifies the deformation of the spacetime due to electromagnetism, so if it is very small it means that the spacetime is well-described by a vacuum (uncharged) black hole. Black holes with mass $M\gtrsim${10}^{9}\text{,}\mathrm{\text{$M_{\odot}$}}$immersed in a magnetic field of order${10}^{11}$~{}\mathrm{G}$ would be needed to reach values of $Q/M$ large enough to be relevant for the spacetime structure. Fields of such strength are expected to be found only in neutron stars. Based on the Wald mechanism, it has been recently proposed that a binary black hole – neutron star could provide a suitable environment to charge the black hole itself Levin et al. (2018). A second case in which charged black holes might occur in the Universe is immediately after the collapse of a compact star when the resulting hole might briefly retain some charge Ray et al. (2003). A similar scenario is the collapse of magnetized stars Nathanail et al. (2017), which was also considered as a candidate for fast-radio bursts Liu et al. (2016). Finally, charged black holes can emerge in more exotic theories associated with “hidden” gauge fields and elementary particles whose charge is a fraction of the electron charge Cardoso et al. (2016).

In spite of the apparently compelling reasons to believe that astrophysical black holes have practically zero net charge compared to their mass, it is still worth studying the source-free Einstein-Maxwell system to advance our comprehension of strong-field gravitation and electromagnetism in this largely unexplored territory. The interplay between electromagnetism and gravity in a highly dynamical spacetime, which can be probed only with numerical investigations, can offer a unique laboratory for both theoretical and more exotic astrophysical studies. For example, the inclusion of charge in highly relativistic collisions of black holes (see, e.g., Sperhake et al. (2008, 2009); Berti et al. (2010); Sperhake et al. (2013) for such studies with zero charge) would advance our understanding in a new direction never explored before. Another interesting application of dynamical electrovacuums is related to cosmic censorship. In a recent series of papers, it was argued that strong cosmic censorship can be violated by electrovacuums with a positive cosmological constant Mo et al. (2018); Dias et al. (2018a, b). In contrast, the case without cosmological constant is not settled yet Cardoso et al. (2018).

The coalescence and merger of charged black holes may present new interesting phenomenology. For instance, Zhang (2016) proposed that the faint potential electromagnetic counterpart to GW150914 Abbott et al. (2016); Connaughton et al. (2016) might have been the result of merger of charged black holes. Another hypothesized mechanism along similar lines invokes magnetic reconnection Fraschetti (2018). Subsequently, Liebling and Palenzuela (2016) tested the idea of Zhang (2016) with relativistic simulations. However, the set-up considered by the authors had some limitations: only equal-mass, equal-charge, non-rotating black holes were studied, and the initial data did not satisfy the constraints of the field equations. A more systematic study of this part of the parameter space of charged black holes requires that one starts with constraint-satisfying initial data for black hole configurations with arbitrary charge, mass ratio, linear and angular momenta.

The most common avenue for generating constraint-satisfying initial data is provided by the $3+1$ decomposition of spacetime Arnowitt et al. (2008a); York (1971). In this approach, one casts the Einstein-Maxwell equations to an initial value problem in which the four-dimensional spacetime is foliated by successive timeslices obtained via the dynamical evolution of the system.333It is worth mentioning that another common approach to building spacetimes in the computer is the generalized harmonic formalism Pretorius (2005b, a). When performing this decomposition, both Maxwell’s and Einstein’s equations are split in two sets: the evolution and the constraint equations. The former move the system forward in time, whereas the latter must be satisfied at all times and must be used to generate the initial data for the evolution. In this paper, we primarily focus on the constraint equations.

Einstein-Maxwell’s theory was first cast in a $3+1$ decomposition by Thorne and MacDonald (1982) and more than 25 years later, Alcubierre et al. (2009) proved that the evolution equations are symmetric hyperbolic, and hence admit a well-posed initial value problem. Moreover, the authors extended the work of Bowen (1985) to generate initial data for electrically charged black holes at a moment of time-symmetry – the spacetime is invariant with respect to time reversal. Recent applications of this formalism are the head-on collisions by Zilhão et al. (2012, 2014a) (the interested reader can find several cogent additional reasons motivating the numerical study of charged black holes in these references). In these works, the authors evolved initial data generated with the same formalism described by Alcubierre et al. (2009) and were mostly interested in comparing the electromagnetic and gravitational emissions. Finally, the same group also investigated numerically the non-linear stability of a Kerr-Newman black hole Zilhão et al. (2014b).

In this paper, we extend the work of Alcubierre et al. (2009) and Zilhão et al. (2014b) to generate initial data for charged, rotating, and moving black holes in a self-consistent way.444We note that the formalism outlined in this paper applies not only to electromagnetism but to any U(1) charge (such as the one described in Cardoso et al. (2016)). We adopt the conformal transverse-traceless formalism Bowen and York (1980) treating the black holes as punctures to solve for the metric, and take advantage of the Reissner-Nordström solution in isotropic coordinates to solve for the electromagnetic fields. This strategy involves two major challenges. The first is that non-linear partial differential equations have to be solved. This can be done only numerically for generic binary black hole configurations. To address this issue we modify the TwoPuncture code Ansorg et al. (2004) to solve the resulting elliptic differential equations. The single-domain pseudo-spectral character of the code results in an accurate solution and it is quickly convergent. The second challenge is that the physical interpretation of the results is not transparent. The parameters given as input for the algorithm (the bare parameters, such as mass and charge) in general are not actual physical quantities of the resulting black holes. Hence, we apply the theory of isolated horizons Ashtekar et al. (2000), which provides a quasi-local machinery for linking the bare black hole parameters with the physical ones and is suitable for simulations. We implement this numerically by modifying the QuasiLocalMeasures code Dreyer et al. (2003).

We structure the paper as follows. In Section II we review the mathematical tools necessary for generating initial data for charged black holes. In particular, we present the $3+1$ decomposition of Einstein-Maxwell’s theory and review the Reissner-Nordström solution in isotropic coordinates and the formalism of isolated horizons. In Section III we solve the constraints with the conformal transverse-traceless technique. Our numerical implementation and tests are detailed in Section IV. Finally, Section V summarizes our findings and describes possible future research directions.

In Appendix A, we prepared a summary of the important equations and steps needed to generate initial data for generic systems of charged black holes. The Appendix provides a distilled overview of the analytic content of this paper. For the reader who is interested only in the gist of the algorithm/equations and the results of our work, we suggest they skip to Appendix A, and then read Sections IV.1 and IV.2, where we present our results.

Notation and conventions

We assume that gravity and electromagnetism are described by Einstein-Maxwell’s theory Wald (1984) and we follow the same notation as in Misner et al. (1973). In particular, we use Einstein’s summation convention and the signature of the metric is $(-,+,+,+)$. We use geometrized units with $G=c=1$, where $c$ is the speed of light in vacuum and $G$ is the gravitational constant. The unit of charge is defined so that the proportionality constant in Coulomb’s law is 1 (for more details, see Jackson (1975)). Indices $a$, $b$, $c$, and $d$ run in the set $\{0,1,2,3\}$, whereas the other Latin letters, such as $i$, $j$ or $k$, run in the set $\{1,2,3\}$ and are referred to as spatial components. Parentheses and square brackets in the indices mean symmetrization and anti-symmetrization, respectively. We also use the abstract index notation Wald (1984). We reserve the symbol $\nabla$ for the four-dimensional covariant derivative associated with the spacetime metric $g_{ab}$, and $D$ for the three-dimensional covariant derivative, compatible with the spatial metric $\gamma_{ij}$. We denote the determinant of these metrics as $g=\det{g_{ab}}$ and $\gamma=\det{\gamma_{ij}}$. We prepend the symbol “${(4)}$” to all the four-dimensional tensors, with exception of the metric $g_{ab}$. For the completely antisymmetric Levi-Civita tensor we use the convention that $\epsilon_{1230}=\sqrt{-g}$, and $\epsilon_{123}=\sqrt{\gamma}$, and denote the Levi-Civita symbol with $\bar{\epsilon}_{ijk}$ or $\bar{\epsilon}^{ijk}$.

II Formalism

In this Section we describe the theoretical tools that we use later to generate initial data for arbitrary configurations of charged black holes. Specifically, in Section II.1 we survey the $3+1$ decomposition of Einstein-Maxwell’s equations, focusing on the constraint equations. Section II.2 reviews the Reissner-Nordström solution for a single charged stationary black hole in isotropic coordinates. Section II.3 summarizes the theory of isolated horizons, which we employ to assign the black hole physical properties: mass, charge and angular momentum.

II.1 $3+1$ decomposition of Einstein-Maxwell

In this paper we study systems described by the source-free Einstein-Maxwell equations Wald (1984)

[TABLE]

where ${}^{(4)}\mathcal{R}_{ab}$ is the Ricci tensor associated with the metric $g_{ab}$, ${}^{(4)}\mathcal{R}={}^{(4)}\mathcal{R}^{a}_{\;a}$, ${}^{(4)}F_{ab}=2\,{}^{(4)}A_{[a,b]}$ is the Maxwell field-strength tensor, with ${}^{(4)}A_{a}$ the four-potential, and ${}^{(4)}{}^{\star}F_{ab}$ is its Hodge dual, defined by

[TABLE]

The electromagnetic stress-energy tensor is

[TABLE]

Solving the coupled Einstein-Maxwell equations in four dimensions is a challenging task. In particular, the form of Equations (1) is not suitable for a numerical solution. Therefore, we adopt the standard $3+1$ decomposition to express the equations as a Cauchy problem, and cast them in a form amenable for numerical integration Baumgarte and Shapiro (2010b).

Assuming that the spacetime is described by a globally hyperbolic Lorentzian manifold $\mathcal{M}$ with metric tensor $g_{ab}$, $\mathcal{M}$ can be foliated by a family of spacelike non-intersecting hypersurfaces $\Sigma_{t}$, taken as level surfaces of a time function $t$. Let $n^{a}$ be the future-directed, timelike unit vector normal to $\Sigma_{t}$. The projection operator along this vector is $n^{a}n_{b}$, whereas the one onto $\Sigma_{t}$ is

[TABLE]

The induced metric on $\Sigma_{t}$, is derived by applying twice the projection operator on $g_{ab}$, which yields

[TABLE]

The induced metric is purely spatial ($\gamma_{ab}n^{b}=0$), it encodes the intrinsic curvature of the hypersurfaces $\Sigma_{t}$ and can be used to define a spatial covariant derivative $D_{i}$ on $\Sigma_{t}$.

Instead of working with the normal vector $n^{a}$, it is convenient to use the normalized time vector

[TABLE]

where $\alpha$ and $\beta^{a}$ are the lapse function and shift vector. With these quantities, the spacetime metric assumes the Arnowitt-Deser-Misner (ADM) form Arnowitt et al. (2008a, b)

[TABLE]

The spatial metric is not sufficient to fully describe the curvature properties of the four-dimensional spacetime. The extrinsic curvature $K_{ab}$ supplies the missing information by expressing how $\Sigma_{t}$ is embedded in $\mathcal{M}$, and is defined as

[TABLE]

Just like the induced metric (which we will also refer to as the three-metric throughout), the extrinsic curvature is purely spatial. The Riemann tensor can be expressed in terms of $\gamma_{ij}$ and $K_{ij}$, and therefore Einstein’s equations can be rewritten in terms of $3+1$ quantities. The resulting 3+1 ADM (à la York) formalism Arnowitt et al. (2008a); York (1971) of general relativity consists of four constraints and twelve evolution equations. The constraints are the direct consequence of the integrability conditions that $\gamma_{ij}$ and $K_{ij}$ have to satisfy to have $\Sigma_{t}$ properly embedded in $\mathcal{M}$. On the other hand, the evolution equations provide a prescription to move from one timeslice to the next provided a gauge choice is made. The evolution equations preserve the constraints: if the constraints are initially satisfied, they will always be satisfied. However, when they are not satisfied, the simulated system is not a solution of the Einstein equations. The same split into evolution equations and constraint equations holds for Maxwell’s theory, too. In complete analogy to Einstein’s theory, Maxwell’s evolution equations preserve the Maxwell constraints, if the constraints are initially satisfied. For this reason, it is important to start with valid, constraint-satisfying initial data. In this work, we focus only on the constraint equations, precisely because our goal is the generation of valid initial data for general relativistic simulations in Einstein-Maxwell theory.

Let ${}^{(4)}T_{ab}^{\text{EM}}$ be the stress-energy tensor, and define

[TABLE]

The Einstein constraints then become Baumgarte and Shapiro (2010b)

[TABLE]

with $\mathcal{R}$ being three-dimensional Ricci scalar associated with $\gamma_{ij}$, and $K$ the trace of the extrinsic curvature. Equation (10a) is known as the Hamiltonian constraint, Equations (10b) as the momentum constraints.

Equations (10) are not the only constraints in Einstein-Maxwell’s theory. As for Einstein’s equations, a $3+1$ split of Maxwell’s equations must be performed.555A more detailed derivation of the three-dimensional Maxwell equations from the four-dimensional ones can be found in the Appendix of Alcubierre et al. (2009) (see also Baumgarte and Shapiro (2003)). First, we introduce the electric and magnetic fields as measured by normal observers with four-velocity $n^{a}$,

[TABLE]

which are both purely spatial ($n_{a}E^{a}=n_{a}B^{a}=0$). The electromagnetic tensor becomes

[TABLE]

and its dual is

[TABLE]

With these decompositions, Maxwell’s equations can be expressed in terms of $3+1$ quantities. As in the case of the Einstein equations, the $3+1$ split leads to evolution and constraint equations. In particular, the electromagnetic constraints are

[TABLE]

The electromagnetic sector couples with the spacetime through the stress-energy tensor $T^{ab}_{\mathrm{EM}}$ which is re-written in terms of the 3+1 variables as

[TABLE]

where $\epsilon^{bcd}=n_{a}\epsilon^{abcd}$. Plugging Equation (15) into the source terms of Equations (9), we find

[TABLE]

which are the familiar electromagnetic energy density and Poynting vector.

II.2 The Reissner-Nordström spacetime

The Reissner-Nordström spacetime Reissner (1916); Nordström (1918) describes an isolated non-rotating black hole with electric charge $q$ and mass $m$ Wald (1984). This solution will be the base of our generalization to charged black hole systems. In Boyer-Lindquist coordinates ($t,r,\theta,\phi$), the Reissner-Nordström metric is given by

[TABLE]

and the electromagnetic potential of the solution is

[TABLE]

In the following Sections we will adopt the puncture approach, so we transform the Boyer-Lindquist coordinates to isotropic ones. In order to do so, we define a new radial coordinate $R$ satisfying

[TABLE]

with $R_{H}=\frac{1}{2}\sqrt{m^{2}-q^{2}}$ the radius of the black hole horizon in isotropic coordinates. The metric then assumes the following form

[TABLE]

with $\delta_{lk}$ the flat Euclidean metric, and $\Psi$ the conformal factor defined as

[TABLE]

As is clear from Equation (20), the spatial metric is manifestly conformally flat in isotropic coordinates. Moreover, there is no magnetic field and the electric field has only an $R$ component

[TABLE]

As a result, the Poynting vector defined in Equation (16b) is identically zero everywhere.

II.3 Isolated horizons

Once the constraint equations are solved, it is important to interpret the physical configuration to which the initial data correspond. This can be achieved by locating the black hole apparent horizons and applying the theory of isolated horizons Ashtekar et al. (2000) (see Ashtekar and Krishnan (2004) for a review). Isolated horizons provide a quasi-local notion of the black hole physical properties. In this Section we review basic identities of the formalism, including, in particular, the electric charge of the horizon, and the electromagnetic field contribution to angular momentum, elements that have not received much attention in numerical relativity applications Dreyer et al. (2003); Schnetter et al. (2006).

Isolated horizons have several desirable features. For instance, they always lie inside the event horizon, to which they reduce for stationary spacetimes, and they imply the existence of a future singularity Penrose (1965); Hawking and Penrose (1970). Most relevant for our purpose, they provide well-defined notions of mass, charge and angular momentum. For spacetimes with suitable symmetries, these quasi-local physical quantities coincide with the global ones defined from conservation laws (for example via ADM integrals), as we verify this explicitly for the Reissner-Nordström case in Appendix B. However, in general, the quasi-local definitions and those at infinity differ Ashtekar et al. (2000). Furthermore, the formalism does not provide a quasi-local definition of linear momentum due to the lack of a meaningful notion of space-translational symmetry in curved spacetime Ashtekar and Krishnan (2004); Krishnan (2002).

Here, we follow closely Dreyer et al. (2003) in using isolated horizons to assign black hole physical parameters. Given a spatial section $\mathcal{S}$ of an isolated horizon, the variables we are interested in are defined as follows. First, the areal radius is given by

[TABLE]

where $\epsilon$ is the area two-form on the 2-surface, given by $\epsilon=\frac{1}{2}\sqrt{\mathbbm{q}}\bar{\epsilon}_{ab}\operatorname{d}\!{x}^{a}\wedge\operatorname{d}\!{x}^{b}$, where $\mathbbm{q}_{ab}$ is the induced metric on the horizon, $\mathbbm{q}=\det\mathbbm{q}_{ab}$, and $\bar{\epsilon}_{ab}$ is the two-dimensional antisymmetric symbol. $\int_{\mathcal{S}}\epsilon$ is the surface area of the horizon.

Next, the definition of the angular momentum is based on an approximate rotational killing vector field $\varphi^{a}$ on the 2-surface Ashtekar et al. (2000)

[TABLE]

where $\omega$ is the form that satisfies the condition $t^{a}\nabla_{a}k^{b}=t^{a}\omega_{a}k^{b}$ for any vector $t^{a}$ tangent to $\mathcal{S}$, with $k^{b}$ being the outgoing future-directed vector normal to $\mathcal{S}$ . By construction of $k^{b}$, $\omega$ always exists Ashtekar et al. (2000). The two terms in the right-hand-side of Equation (24) are the gravitational and electromagnetic contribution to the horizon angular momentum.

The charge is defined by means of Gauss’s law

[TABLE]

and finally, the gravitational mass of the isolated horizon is given by

[TABLE]

For Kerr-Newman black holes, this formula perfectly reduces to the equation that relates total mass, irreducible mass, charge and angular momentum Ashtekar et al. (2001).

The definitions of angular momentum and charge involve four-dimensional quantities, but during simulations with the $3+1$ formalism, it is more convenient to use 3+1 variables. In Dreyer et al. (2003), it was shown that the gravitational contribution to the horizon angular momentum can be computed using an ADM-like formula

[TABLE]

where $R^{a}$ is the spatial unit vector normal to $\mathcal{S}$. The electromagnetic component of the angular momentum $J_{\mathcal{S}}^{\text{EM}}=J_{\mathcal{S}}-J_{\mathcal{S}}^{\text{GR}}$ depends on both ${}^{(4)}A$ and ${}^{(4)}{}^{\star}F$. The first is directly accessible if instead of the electric and magnetic fields one evolves the vector potential Del Zanna et al. (2003); Giacomazzo et al. (2011); Etienne et al. (2012b, 2015); Fragile et al. (2018), whereas the second has components

[TABLE]

When integrated over a spatial 2-surface the term $2n_{[a}B_{b]}$ does not contribute because $n_{a}=(-\alpha,0,0,0)$. Therefore, the electromagnetic contribution to the horizon angular momentum becomes

[TABLE]

where $\epsilon_{abc}=n^{d}\epsilon_{abcd}$. By use of Equation (28), Equation (25) for the charge becomes

[TABLE]

These definitions provide a complete characterization of black holes during a general relativistic simulation with the $3+1$ decomposition. An example of how the integrations above are performed is in Appendix C.

III Solving the constraint equations

To solve the constraint equations we adopt the conformal transverse-traceless approach, also referred to as Bowen-York technique Bowen and York (1980). The goal of this method is to expose and specify degrees of freedom containing physical information about the system by applying conformal transformations on the spatial quantities, and working directly on the conformal variables instead of the physical ones.

The first step in the method is to conformally decompose $\gamma_{ij}$ by introducing the conformal factor $\psi$ and metric $\bar{\gamma}_{ij}$

[TABLE]

In the following we use an overbar to indicate conformal quantities.

A common assumption when generating multiple black hole initial data is that the spatial metric is conformally flat Cook (2000); Alcubierre et al. (2009); Baumgarte and Shapiro (2010b); Lousto et al. (2012). In other words, we fix the conformal three-dimensional metric $\bar{\gamma}_{ij}$ to be the flat Euclidean metric $\delta_{ij}$ (in Cartesian coordinates). This choice greatly simplifies computations and it is a good approximation for the systems we are interested in studying, in spite of the fact that conformally flat spatial slices of the Kerr metric do not exist Garat and Price (2000). Conformal flatness limits the maximum equilibrium value that the black hole dimensionless spin can attain Lovelace et al. (2012); Lousto et al. (2012), but values of order $0.9$ are completely achievable. Thus, we do not anticipate this approximation to impose severe constraints on the equilibrium values of the black hole spin and charge. Considering what happens in the uncharged case Brandt and Seidel (1995); Gleiser et al. (1998), we expect that conformal flatness will generate initial data with spurious gravitational (and electromagnetic) radiation in the charged black hole cases, too. Nonetheless, this is not a major concern since in dynamical simulations the system is evolved until this “junk” radiation propagates away, and the fields relax to their quasi-equilibrium values.

In addition to the conformal decomposition of the metric, it is also useful to transform the extrinsic curvature $K_{ij}$ by separating it into its traceless $A_{ij}$ and trace ($K=K^{i}_{\;i}$) parts

[TABLE]

Following standard practice, we adopt the maximal slicing condition $K=0$ Smarr and York (1978), and introduce a conformal, traceless extrinsic curvature $\bar{A}_{ij}$ as

[TABLE]

Then, $\bar{A}_{ij}$ can be split into a transverse-traceless and a longitudinal part

[TABLE]

We set $\bar{A}^{ij}_{\text{TT}}=0$, which corresponds to suppressing the radiative degrees of freedom, so

[TABLE]

The longitudinal part can always be expressed in terms of a vector $V$ as

[TABLE]

where Cartesian coordinates are adopted. Going back to Equation (33), the extrinsic curvature is given by

[TABLE]

We already exploited much of the freedom we had in specifying variables during the previous steps. Under these assumptions, we just need the vector $V^{i}$ and the conformal factor $\psi$ to fully determine $\gamma_{ij}$ and $K_{ij}$, and the constraint Equations (10) take the form

[TABLE]

where $\nabla^{2}=\partial_{k}\partial^{k}$.

Next, we turn to the electromagnetic sector of the problem. We rescale the electromagnetic fields as in Alcubierre et al. (2009)

[TABLE]

The factor $\psi^{6}$ is chosen in order to have $D_{i}E^{i}=\psi^{-6}\partial_{i}\bar{E}^{i}$, where we used the fact that for any vector $v^{i}$ it holds true that $D_{i}v^{i}=\gamma^{-1/2}\partial_{i}\left(\sqrt{\gamma}v^{i}\right)$. The Maxwell constraints (14) read

[TABLE]

These equations do not depend on the conformal factor $\psi$, so the electromagnetic constraints can be solved independently from the spacetime ones. Moreover, the equations are linear; hence we can superpose solutions.

Having fixed the conformal scalings of the $E^{i}$ and $B^{i}$ fields, the source terms $\mathcal{E}$ and $S^{i}$ of the Einstein constraints conformally transform as

[TABLE]

where

[TABLE]

With these redefinitions, the Einstein constraints become

[TABLE]

The problem is now greatly simplified because the momentum constraints do not depend on $\psi$, are linear in $V^{i}$, and along with the Hamiltonian constraint have decoupled from the Maxwell constraints.

Next, we exploit the linearity of Equation (43b) by decomposing $V^{i}$ as

[TABLE]

where $V^{i}_{0,\text{GR}}$ solves the homogeneous Equation (10b) (when $\bar{S}^{i}=0$), and ${V}^{i}_{\text{EM}}$ the inhomogeneous one.666The subscript [math] does not indicate any component, but it reminder that the field is a solution of the homogeneous equation. The first term does not contain any reference to the electromagnetic sector of the problem. Thus, as in Ansorg et al. (2004), we choose

[TABLE]

with $R_{n}=\left\lvert\mathbf{x}-\mathbf{x}_{n}\right\rvert$ the Euclidean coordinate distance from puncture $n$, where $\mathbf{x}_{n}$ is the location of the $n$-th puncture, and $P_{n}^{i}$ and $S^{k}_{n}$ are its linear and angular momenta, respectively. Equation (45) solves the homogeneous version of Equation (43b), and it is known that for suitable single black hole solutions $P_{\text{ADM}}^{i}=P^{i}$ and $J_{\text{ADM}}^{i}=S^{i}$, with $P_{\text{ADM}}$ and $J_{\text{ADM}}$ being the ADM linear and angular momenta evaluated at infinity Bowen and York (1980); Ansorg et al. (2004), respectively. By use of the decomposition (44), the momentum constraints further reduce to three decoupled linear equations for $V^{i}_{\text{EM}}$, effectively replacing Equation (43b) with

[TABLE]

We also manipulate the Hamiltonian constraint (43a) further by separating the singular part of the conformal factor from the finite one $u$, motivating our ansatz based on the conformal factor of the Reissner-Nordström spacetime in Equation (21),

[TABLE]

We introduce the following abbreviations for compactness

[TABLE]

Therefore, the conformal factor becomes

[TABLE]

Equation (47) is essentially an ansatz that states that our solution is a superposition of Reissner-Nordström black holes plus corrections (in $u$), which parallels what is performed in the uncharged case Brandt and Brügmann (1997).

Expanding Equation (43a), we reach

[TABLE]

In deriving the last expression, we used the fact that the Laplacian of $\eta$ and $\varphi$ is zero. Equation (50) is a second order, non-linear elliptic partial differential equation in $u$ that depends on $V_{\text{EM}}^{i}$ through the term $\bar{A}_{ij}\bar{A}^{ij}$. Now, the momentum and Hamiltonian constraints (10) have been re-expressed as elliptic equations (46), (50) for $V_{\text{EM}}^{i}$ and $u$. The associated boundary conditions are found from the assumption of asymptotic flatness so that $u$ and $V_{\text{EM}}^{i}$ have to go to zero at spatial infinity. In this paper we assume that $u$ and $V_{\text{EM}}$ are regular everywhere, and thus they can be found with standard numerical methods that can solve Equations (46), (50).

The problem of generating valid initial data for multiple charged black holes is now reduced to solving Equations (46) and (50), which is done once Maxwell-compliant electromagnetic fields are found. In this paper, we assume that each puncture is endowed with a Reissner-Nordström electromagnetic field, and hence the total conformal electric field is a superposition of Reissner-Nordström electromagnetic fields in isotropic coordinates, i.e.,

[TABLE]

nwhere $\hat{R}_{n}$ is the radial unit vector centered on the $n-$th puncture. In the case of a single, non-rotating black hole with zero linear momentum, our choice of Reissner-Nordström fields exactly produces a spatial slice of that solution, since the constraints are solved by $V_{\text{GR}}=V_{\text{EM}}=0$, and $u=0$ [so $\psi=\Psi$, where $\Psi$ is given in Equation (21)]. For systems of spinning black holes with linear momenta, the superposition of Reissner-Nordström fields is a first approximation to the equilibrium electromagnetic field generated by these configurations. As for the gravitational fields generated in the puncture approach (and the gauge fields), we expect that the time evolution will relax our electromagnetic-field initial data to their quasi-equilibrium values on a light-crossing timescale. An advantage of choosing Reissner-Nordström electromagnetic fields is that they allow for a clear description of each black hole in the system with a specific charge, whose isolated horizon value $Q_{\mathcal{S}}$ equals the “bare” charge entering Equation (51). In addition, since there is no magnetic field, the source term of Equation (46) vanishes, and so $V_{\text{EM}}=0$ (even for multiple black holes with linear and angular momenta). Thereby, this choice ensures that there are no electromagnetic contributions to the extrinsic curvature, implying that the parameters entering Equation (45) can be still interpreted as $P_{\text{ADM}}^{i}=P^{i}$ and $J_{\text{ADM}}^{i}=S^{i}$.

The choice of Reissner-Nordström electromagnetic fields is by no means unique. Another possibility is Kerr-Newman fields in quasi-isotropic coordinates. We present a detailed discussion of this case and the complexities associated with it in Appendix D.

IV Numerical implementation

We implement the formalism outlined in the previous Sections by modifying the TwoPunctures Ansorg et al. (2004) and QuasiLocalMeasures open-source codes Dreyer et al. (2003). The software is run within the Cactus infrastructure Goodale et al. (2003) and all physical variables are interpolated on a Carpet grid Schnetter et al. (2004); Paschalidis et al. (2013b). Black hole apparent horizons are found with AHFinderDirect Thornburg (2004).

The main component in our software stack is TwoChargedPunctures, which is used to generate initial data for two punctures located at $(\pm b,0,0)$ given the bare black hole properties ($M_{n}$, $Q_{n}$, $P^{i}_{n}$, $S^{i}_{n}$). This code implements a pseudo-spectral collocation method that solves the constraint equations (46) and (50) to find $u$ and $V_{\text{EM}}^{i}$.

In what follows, we adopt Reissner-Nordström electromagnetic fields. Since there is only an electric field, $\bar{S}^{i}=0$ in Equation (46), and the momentum constraint is trivially satisfied by $V_{\text{EM}}^{i}=0$. Hence, we only need to solve the Hamiltonian constraint (50).

TwoChargedPunctures implements a single domain pseudo-spectral method that covers all $\mathbb{R}^{3}$ with spatial infinity on the grid. This region is parametrized by the coordinates $(A,B,\phi)$, with $A,B\in[-1,1]$ and $\phi\in[0,2\pi]$. To be more specific, the code uses a system of bispherical coordinates that transform to the usual Cartesian ones with the law777This parametrization is slightly different compared to what is done in Ansorg et al. (2004). The spectral expansion used here treats $A$ and $B$ on equal footing, i.e., the spectral decomposition in $A$ and $B$ uses the same Chebyshev polynomial basis, unlike what is reported in Ansorg et al. (2004).

[TABLE]

where the $x$ axis is along the line connecting the two punctures. Equations (52) describe a set of cylindrical-like coordinates around the $x$ axis with a radius that depends on both $A$ and $B$.

The coordinates $(A,B,\phi)$ live on a compact grid where spatial infinity corresponds to $A=1$, which makes it straightforward to impose the desired outer boundary conditions ($u\to 0$ at infinity). This condition is enforced by solving the equations for an auxiliary variable $\mathpzc{U}$ defined as $u=(A-1)\mathpzc{U}$. The code expands $\mathpzc{U}$ in Chebyshev polynomials along $A$ and $B$, and adopts a Fourier basis along $\phi$. The coordinates are discretized with $n_{A}$, $n_{B}$ and $n_{\phi}$ grid points chosen as the zeros of Chebyshev polynomials $T_{n_{A}}(x)$, $T_{n_{B}}(x)$ and of the sine function $\sin{(n_{\phi}\phi)}$. The coefficients of the spectral expansion are found by evaluating the relevant equation on the collocations points and solving the corresponding multidimensional non-linear system with a modified Newton-Raphson method Barrett et al. (1994) (more details on how this is done can be found in Section II of the original paper Ansorg et al. (2004)). We consider the equations to be solved, when the residuals are smaller than a threshold value. To choose this threshold value, we solve for increasingly smaller values of this threshold and compute the ADM and the horizon masses. When these masses have converged to within one part in $10^{6}$, we consider the solution converged.

With the equations solved and $u$ known, TwoChargedPunctures reverts back to the physical fields using Equations (37), (39) (44), (45) and (47). We then spectrally interpolate the physical fields on a Carpet grid where AHFinderDirect is subsequently run to locate the apparent horizons. Once the horizons are found, we compute mass, charge and angular momentum of each black hole with our version of QuasiLocalMeasures, which we call QuasiLocalMeasuresEM, and which implements the formalism of isolated horizons for the full Einstein-Maxwell theory as reviewed in Section II.3. Moreover, having the spectral expansion of the fields we can interpolate them at a very large radius to compute the ADM mass, the linear and angular momenta.

IV.1 Code validation

We validate our approach and numerical implementation with a series of tests that are presented in this section.

We report our results in terms of the input bare mass $M$ of the punctures, which is the only mass known a priori. In all the runs, we confine the black hole in a region where the Carpet grid resolution is $\Delta x_{i}=$0.0078$~{}M$, which usually guarantees that the diameter of the horizon is resolved by about $100$ points, making it easily found by AHFinderDirect. We also fix the resolution of the AHFinderDirect grid to be $79$ points in the azimuthal direction and $39$ in the meridional direction. We have confirmed that the resolution on the AHFinderDirect grid has negligible impact in our results. In the cases presented here, doubling the AHFinderDirect grid resolution introduces a variation in the computed parameters of order $0.01\text{\,}\mathrm{\char 37\relax}$. We compute ADM integrals by spectrally interpolating our fields on a sphere of radius $10000~{}M$, and discretized with $256$ points in both the meridional and azimuthal directions.

As a first test, we made sure that our modified code with zero charge, TwoChargedPunctures, produces the same output as the standard open-source TwoPunctures code. This is not a trivial test because the equations used in our code and in the original one are different, having different numerical properties, even though they are mathematically equivalent. In particular, our formulation is more susceptible to numerical instabilities due to the finite-arithmetic error in regions close to the puncture. The reason for this is that our equations have terms that are not present in the original code, but that should perfectly cancel out when $Q=0$. Such a numerical cancelation near the punctures is not trivial. However, the result of the test with different spectral resolutions shows that the two implementations agree at the round-off-error level for punctures with no charge.

Another key test that our code successfully passes consists in recovering the only conformally flat analytical solutions known: the Reissner-Nordström and the case of two black holes with the same charge-to-mass ratio (see Appendix E for more details), both of which are found with $u=0$. We find the solution $u=0$ is recovered to machine precision everywhere outside the horizons, and it is non-identically zero only very close to the punctures, again due to numerical precision.

The next test for TwoChargedPunctures is reproducing the numerical solution found by Alcubierre et al. (2009) for two non-rotating black holes with opposite charge-to-mass ratio starting at rest. Figure 1 reports the value of $u$ along the $x$, $y$ and $z$ axes for a system of two punctures with the same mass but opposite charge ($Q_{1}=-Q_{2}=$0.5$~{}M$). We graphically superposed our plot with Figure 1 in Alcubierre et al. (2009), finding perfect agreement.

Continuing the progression of complexity in the considered systems, we generate a single puncture with angular momentum, but no linear momentum, and one with linear momentum but no angular momentum (Figures 2 and 3, respectively). In these single-black hole cases, we compare the horizon mass with the ADM mass measured at infinity and we find agreement of order $0.1\text{\,}\mathrm{\char 37\relax}$ even with resolution as low as $n=16$. The same is true for the ADM angular momentum and the horizon spin, as computed with QuasiLocalMeasuresEM. We repeated these two tests by aligning the linear and angular momentum vectors once along the $x$ direction and once along the $z$ direction to ensure that the built-in asymmetry in the coordinates [Equations (52)] does not spoil expected symmetries in symmetric configurations. By doing this, we find that the solutions are rotationally invariant to better than one part in $10^{6}$ for a resolution $n=32$ or higher.

IV.2 Convergence

Finally, we considered the generic system shown in Figure 4. This is formed by two equal-mass black holes with charge $Q_{1}=-$0.3$~{}M$ and $Q_{2}=$0.5$~{}M$. Both black holes are spinning with angular momentum $S^{z}_{1}=S^{z}_{2}=$0.5$~{}M^{2}$. The black holes also have linear momentum $P_{1}^{x}=P_{2}^{z}=-$0.5$~{}M$. The solution for $u$ for this system is depicted in Figure 5. With QuasiLocalMeasuresEM, we find that the quasi-local angular momenta (charges) agree with their bare counterparts to within one part in $10^{4}$ ($10^{8}$). We find that the mass of the first horizon is 1.187$~{}M$ and the second is 1.202$~{}M$. The total (ADM) mass of the system is $$2.337$~{}M$, and the difference between this value and the sum of the individual masses is the binding energy plus contribution from “junk” radiation.

This system is used to study the self-convergence properties of the code. In particular, we consider the maximum relative error of $u$ with respect to a reference solution at high resolution $N$. For this, we sampled $u$ on a set of points $\mathcal{T}$ and computed the infinity norm

[TABLE]

We choose $\mathcal{T}$ as the set of points where spheres of radii $1~{}M$, $2~{}M$, $5~{}M$, $10~{}M$, $100~{}M$, and $1000~{}M$ intersect the coordinate axes for $x>0,\,y>0,\,z>0$.

We set as a reference solution ($N$) one obtained at high-resolution with $n_{A}=n_{B}=n_{\phi}=n=64$, which is between $n=50$ and $n=70$ that were used for self-convergence tests in the original TwoPunctures code Ansorg et al. (2004). Here, we simply choose resolutions which are multiples of 4, but our results do not depend on this choice. Our convergence test (Figure 6) shows that the algorithm is robust; $u$ quickly converges to its high-resolution value. The code converges approximately at sixth-order. We also verified that the code exhibits the same convergence properties when we repeat the convergence test with $Q_{1}=Q_{2}=0$, which also agree with the convergence properties of the original TwoPunctures code Ansorg et al. (2004). The convergence of $u$ also results in excellent convergent behavior for both the ADM mass and momenta and the horizon properties as computed by QuasiLocalMeasuresEM.

V Conclusions and future work

Gravitational waves offer new opportunities to study the Universe that are not accessible with electromagnetic or neutrino astronomy. In this landscape, numerical-relativity simulations are a powerful tool to gain insight into the properties and the characteristics of both the waves and their sources. The majority of numerical-relativity simulations of black holes to-date do not treat the electric charge. This is because it is believed that astrophysically relevant black holes should have a charge which is negligibly small compared to the mass. For this reason, there are no studies of highly dynamical electrovacuum spacetimes that involve the inspiral and merger of binary black holes with charge and spin. Nevertheless, electrovacuum spacetimes are of great interest, having both a theoretical appeal and exotic astrophysical applications.

In this paper, we initiated an effort toward solving the coupled Einstein-Maxwell equations in a dynamical and fully general relativistic regime. The first step to perform this type of simulations is the generation of valid initial data. Here, we employed the conformal transverse-traceless approach to build a formalism for generating initial data for multiple black holes with charge, angular and linear momenta. Moreover, we applied the theory of isolated horizons to attribute the physical mass, charge and angular momentum to the horizon, providing a solid understanding of the physical content of our initial data. We implemented the formalism in a software based on the TwoPuncture and the QuasiLocalMeasures open-source codes, verifying our implementation with a series of tests involving analytical or previously-known results. The algorithm was found to recover the expected solutions and showed excellent convergence properties.

With the valid initial data for charged, rotating and moving punctures it is now possible to simulate dynamical evolution of several systems that have never been taken in consideration, such as ultra-relativistic head-on collision, and the quasi-circular or eccentric inspiral and merger of two black holes. As a first application of the formalism outlined in this paper we plan to study in the near-future the case of charged and spinning black holes in quasi-circular orbit. Some of these simulations are already underway, and will be presented in forthcoming work.

Acknowledgements.

We are indebted to the authors of the open-source software that we used: Cactus, Carpet, TwoPunctures, AHFinderDirect and QuasiLocalMeasures. We thank V. Cardoso, S. Gralla, L. Lehner, U. Sperhake, J. R. Westernacher-Schneider, and M. Zilhão for useful discussions. Computations were performed on the Ocelote cluster at The University of Arizona.

Appendix A Algorithm and important equations

In this Appendix, we sketch the algorithm and summarize the important equations to generate initial data for $3+1$ evolutions of arbitrary systems of $\mathcal{N}$ black holes with electric charge, linear and angular momenta using the conformal transverse-traceless decomposition. In the following $n$ is used to index the n-th black hole in the system that is $n\in\{1,\dots,\mathcal{N}\}$. Unless otherwise specified sums in this Appendix are over all punctures. We also assume that each black hole is endowed with Reissner-Nordström electromagnetic fields $(\bar{E}^{i},\bar{B}^{i})$ associated with electric charge $Q_{n}$. The steps in generating the initial data are as follows:

Choose the bare parameters $M_{n},Q_{n},S_{n}^{i},P_{n}^{i},\mathbf{x}_{n}$ for each black hole, respectively representing mass, charge, angular momentum, linear momentum, and position. 2. 2.

Compute the conformal electromagnetic fields $(\bar{E}^{j}_{n},\bar{B}^{j}_{n})$ for each black hole. Under the assumption of Reissner-Nordström fields, we obtain

[TABLE]

with $R_{n}=\left\lvert\mathbf{x}-\mathbf{x}_{n}\right\rvert$ the Euclidean coordinate distance from puncture $n$, and $\hat{R}_{n}^{i}$ the corresponding unit vector. Then, superpose the conformal electromagnetic fields of all black holes,

[TABLE] 3. 3.

Solve the inhomogeneous momentum constraint for $V_{\text{EM}}^{i}$

[TABLE]

with

[TABLE]

and imposing as a boundary condition that $V_{\text{EM}}^{i}\to 0$ at spatial infinity.

Given our choice for the electromagnetic fields [Equation (54)], $\bar{S}^{i}=0$, so $V_{\text{EM}}=0$ is a solution of the momentum constraint (58). 4. 4.

Compute the total auxiliary vector $V^{i}$,

[TABLE]

with

[TABLE] 5. 5.

Solve the Hamiltonian constraint for $u$, imposing $u\to 0$ at spatial infinity

[TABLE]

with

[TABLE] 6. 6.

With $\psi$ now known, compute the physical fields that are necessary for the evolution

[TABLE] 7. 7.

Find the isolated horizons ${\mathcal{S}}_{n}$ and compute the associated physical properties

[TABLE]

where ${}^{(4)}{}^{\star}F$ is the dual of the electromagnetic tensor, $\epsilon$ is the horizon surface $2$-form, ${}^{(4)}A$ is the electromagnetic vector potential, $\varphi$ is the approximate rotational Killing vector on $\mathcal{S}_{n}$, and $\omega$ is defined in the main text [see Equation (24)]. $Q_{\mathcal{S}_{n}},R_{\mathcal{S}_{n}},J_{\mathcal{S}_{n}}$ and $M_{\mathcal{S}_{n}}$ are respectively the charge, radius, angular momentum and mass of the n-th horizon.

Appendix B Isolated horizon in the Reissner-Nordström solution

The goal of this Appendix is to show that the formalism of isolated horizons produces the expected black hole properties in the case of the Reissner-Nordström solution. This can be proven starting from metric (17), which we rewrite here for convenience

[TABLE]

with electromagnetic potential

[TABLE]

In this case, a spherical surface with coordinate radius $r_{+}=m+\sqrt{m^{2}-q^{2}}$ is a Killing horizon, which implies that it is an isolated horizon. This is because every Killing horizon which is topologically $S^{2}\times\mathbb{R}$ is an isolated horizon Ashtekar et al. (2000). Therefore, the metric $\mathbbm{q}_{ab}$ induced on the spatial section of the horizon is simply the metric on a spherical surface [$\operatorname{d}\!{s}^{2}=\mathbbm{q}_{ab}\operatorname{d}\!{x}^{a}\operatorname{d}\!{x}^{b}=r^{2}(\operatorname{d}\!{\theta}^{2}+\sin^{2}\theta\operatorname{d}\!{\phi}^{2})$], and the value of $R_{\mathcal{S}}$ defined by Equation (23) coincides with $r_{+}$ itself, since the radial coordinate in Equation (77) is the areal radius. In this case, the rotational vector $\varphi$ in (24) is taken to be the generator of the azimuthal symmetry on the sphere, which is also a Killing vector of the entire spacetime. Hence, we find that $\varphi\cdot{}^{(4)}A=0$ as ${}^{(4)}A$ has only a temporal components and $\varphi$ only spatial. Moreover, since the future-directed vector $k^{a}$ orthogonal to $\mathcal{S}$ has only radial and temporal component, and any $t^{a}$ tangent to $\mathcal{S}$ has only azimuthal and meridional components, $t^{a}\nabla_{a}k^{b}=0$. By construction, we also have $t^{a}\nabla_{a}k^{b}=t^{a}\omega_{a}k^{b}=0$, which implies that $\omega_{a}=0$, because the equation is zero for each $t^{a}$. Hence, by use of Equation (24), we conclude that $J_{\mathcal{S}}=0$.

To compute charge and mass, we need the electromagnetic tensor, which is given by

[TABLE]

and its dual

[TABLE]

The integration of ${}^{(4)}{}^{\star}F/4\pi$ over any sphere of coordinate radius $r$ results in exactly $q$, so Equation (25) implies $Q_{\mathcal{S}}=q$.

Finally, from Equation (26) the horizon mass is

[TABLE]

For a Reissner-Nordström black hole $m$, $q$ are interpreted as the spacetime total energy and electric charge, respectively Wald (1984). Therefore, in this case, the bare mass (charge), the isolated horizon mass (charge), and the physical mass (charge) all coincide.

Appendix C Computing the charge of an isolated horizon

In this Appendix we discuss how we perform the computation of the horizon charge. To compute the charge of the horizon, we need to perform the following integration (see Section II.3):

[TABLE]

This quantity is coordinate-independent, so choosing Cartesian coordinates $(x^{a})=(x,y,z)$, we can write

[TABLE]

with $\gamma$ determinant of the spatial metric. We introduce a parametrization of $\mathcal{S}$ with polar coordinates $(\theta,\phi)$ around the origin $(x_{0},y_{0},z_{0})$,

[TABLE]

with $s(\theta,\phi)$ suitable smooth function. This is always possible since by hypothesis $\mathcal{S}$ has spherical topology and by construction $Q_{\mathcal{S}}$ does not depend on the parametrization. Then, the first term in Equation (83) can be written as

[TABLE]

where $J_{xy}(\theta,\phi)$ is the Jacobian of the transformation (84) involving the coordinates $x$ and $y$

[TABLE]

The remaining terms in Equation (83) are dealt with accordingly.

In QuasiLocalMeasuresEM, we use the parametrization $s(\theta,\phi)$ provided by AHFinderDirect, and we compute the derivatives in the Jacobians using a centered, second-order accurate finite-difference scheme.

Appendix D Kerr-Newman spacetime

In this Appendix we review the Kerr-Newman spacetime and discuss challenges associated with using the Kerr-Newman electromagnetic fields as source terms in the Hamiltonian and momentum constraints.

The Kerr-Newman black hole with mass $m$, electric charge $q$, and angular momentum $am$ in Boyer-Lindquist coordinates $(t,r,\theta,\phi)$ is Wald (1984)

[TABLE]

with

[TABLE]

The electromagnetic vector potential is

[TABLE]

Following the usual procedure for generating puncture initial data, we transform to quasi-isotropic coordinates by introducing a new radial coordinate $R$ as in Zilhão et al. (2014b)

[TABLE]

with $R_{H}=\frac{1}{2}\sqrt{m^{2}-a^{2}-q^{2}}$ radius of the black hole horizon in the new coordinate system. The metric takes now the form

[TABLE]

where

[TABLE]

where $\Psi$ is the conformal factor, $l,k\in\{R,\theta,\phi\}$, and $\alpha$, $\beta$, $\gamma$ and $h$ functions of $(R,\theta,\phi)$, with

[TABLE]

The non-zero components of the electromagnetic fields are888 Our expression for $E^{\theta}$ differs from the corresponding one in Equation (3.5) of Zilhão et al. (2014b) by a factor of $r/R$. We find that the electric field components listed in Zilhão et al. (2014b) do not satisfy Maxwell’s equations, and that Gauss’s law yields a value for the charge that is correct for spherical surfaces, but the value is different on non-spherical surfaces, e.g. ellipsoidal ones. We have checked that our electric fields satisfy Maxwell’s equations, and, as a result, Gauss’s law yields the correct electric charge even on non-spherical surfaces. We conclude that $E^{\theta}$ in Zilhão et al. (2014b) has a typographical error.

[TABLE]

The conformal fields are obtained by scaling by $\sqrt{\gamma}=\Upsilon\rho^{-4}R^{-1}\sin\theta$

[TABLE]

In these coordinates, the conformal fields are regular for $R\to 0$ (in this limit $\rho\sim r\sim 1/R$).

However, Equations (46) and (50) are in Cartesian coordinates. Transforming to Cartesian coordinates as in flat spacetime, the conformal fields are obtained as

[TABLE]

where here $i\in\{x,y,z\}$ and the factor of $R^{2}\sin\theta$ is the determinant of the Jacobian of the transformation and ensures that the resulting fields $\bar{E}^{i}$ and $\bar{B}^{i}$ satisfy the Maxwell constraints

[TABLE]

In these coordinates, the fields are singular when $x,y,z\to 0$. Given this singular behavior, $V_{\text{EM}}^{i}$ is expected to be singular as well near the punctures because the source of the momentum constraint (46) diverges with a high power of R. This is precisely what we find when we implement our algorithm with the Kerr-Newman electromagnetic fields. In particular, for a single Kerr-Newman black hole without linear momentum, the singular source terms are $\bar{S}^{x}$ and $\bar{S}^{y}$, which at leading order for $x,y,z\to 0$ scale as

[TABLE]

A possible approach to dealing with the singular source would be to separate the singular part of the solution from the regular one, as is done for the Hamiltonian constraint. However, this approach typically requires a known analytic solution, and this does not seem possible within the conformal flatness approximation, because the Kerr-Newman solution does not admit conformally flat spatial slices. In future work, we will explore potential solutions to these challenges by lifting the conformal flatness approximation.

Appendix E Generalized Majumdar-Papapetrou

In this Appendix we show that our formalism recovers spatial slices of a generalized Majumdar-Papapetrou’s solution found by Alcubierre et al. (2009) when each black hole is at rest, non-spinning and all black holes have the same charge-to-mass ratio. This happens because under these assumptions the momentum constraint is trivially satisfied, and the Hamiltonian one is solved by $u=0$, as we verify in what follows.

Given our definitions of $\eta$ and $\varphi$ [Equations (48)], if the charge-to-mass ratio is fixed to $\lambda$ for every black hole, then $\varphi=\lambda\eta$. Moreover, with our choice of Reissner-Nordström fields, there are no magnetic fields, so the electromagnetic energy is $8\pi\bar{\mathcal{E}}=4\partial_{a}\varphi\partial^{a}\varphi$, where the factor of $4$ arises from the fact that $\varphi$ is not the electrostatic potential but it is half of it. Plugging the ansatz $u=0$ into the Hamiltonian constraint [Equation (50)] yields

[TABLE]

But, $\psi=\sqrt{\kappa^{2}-\phi^{2}}$, thus, multiplying the last equation by $\psi^{2}$, and expressing the derivatives of $\psi$ in terms of $\kappa$, $\phi$ and their derivatives yields

[TABLE]

Plugging $\kappa=1+\eta=1+\varphi/\lambda$, and $\partial_{a}\kappa=\partial_{a}\varphi/\lambda$ into this last expression, after some algebra we find that the Hamiltonian constraint is satisfied. If we choose $\lambda=1$ we find

[TABLE]

which describes a spatial slice of the Majumdar-Papapetrou spacetime with $N$ extremal black holes Majumdar (1947); Papapetrou (1945).