# Accurate Simulations of Binary Black Hole Mergers in Force-free Electrodynamics

Abstract: We provide additional information on our recent study of the electromagnetic emission produced during the inspiral and merger of supermassive black holes when these are immersed in a force-free plasma threaded by a uniform magnetic field. As anticipated in a recent letter, our results show that although a dual-jet structure is present, the associated luminosity is ~100 times smaller than the total one, which is predominantly quadrupolar. Here we discuss the details of our implementation of the equations in which the force-free condition is not implemented at a discrete level, but rather obtained via a damping scheme which drives the solution to satisfy the correct condition. We show that this is important for a correct and accurate description of the current sheets that can develop in the course of the simulation. We also study in greater detail the three-dimensional charge distribution produced as a consequence of the inspiral and show that during the inspiral it possesses a complex but ordered structure which traces the motion of the two black holes. Finally, we provide quantitative estimates of the scaling of the electromagnetic emission with frequency, with the diffused part having a dependence that is the same as the gravitational-wave one and that scales as L^(non-coll)_(EM) ≈ Ω^((10/3)–(8/3)), while the collimated one scales as L^(coll)_(EM) ≈ Ω^((5/3)–(6/3)), thus with a steeper dependence than previously estimated. We discuss the impact of these results on the potential detectability of dual jets from supermassive black holes and the steps necessary for more accurate estimates.

## Summary (5 min read)

### 1. INTRODUCTION

- In contrast, magnetic fields generated by the circumbinary accretion disk could play an important role and the dynamics of the plasma in the inner region can then be described within the force-free (FF) approximation.
- These physical conditions are indeed similar to those considered in the seminal investigations of BH electrodynamics of Blandford and Znajek (Blandford & Znajek 1977), who addressed the question of whether the rotational energy of an isolated BH can be extracted efficiently by a magnetic field.

### 2. EVOLUTION EQUATIONS

- The authors solve the combined system defined by the Einstein and Maxwell equations and model either an isolated rotating BH or a BH binary inspiralling in quasi-circular orbits.
- More specifically, the authors solve the Einstein equations Rμν − 1 2 Rgμν = 8πTμν, (1) where Rμν , gμν , and Tμν are the Ricci, the metric, and the stress-energy tensors, respectively.
- Such scalar fields are initialized to zero, but are driven into evolution as soon as violations of the EM constraints are produced.
- In the rest of their discussion the authors will use the expression “electrovacuum” to denote the case when currents and charges of the Maxwell equations are zero.

### 2.1. The Einstein Equations

- For the solution of the Einstein equations the authors make use of a three-dimensional finite-differencing code that adopts a conformal-traceless “3 + 1” BSSNOK formulation of the equations (see Pollney et al. 2007 for the full expressions in vacuum and Baiotti et al. 2008 for the case of a spacetime with matter).
- The code is based on the CactusComputational Toolkit (Allen et al. 2000) and employs adaptive mesh-refinement techniques via the Carpet-driver (Schnetter et al. 2004).
- For compactness the authors will not report here the details regarding the adopted formulation of the Einstein equations and the gauge conditions used, which can, however, be found in Pollney et al. (2007, 2011).
- The authors also note that recent developments, such as the use of eighth-order finite-difference operators or the adoption of a multiblock structure to extend the size of the wave zone, have been recently presented in Pollney et al. (2009, 2011).
- Here, however, in order to limit the computational costs and because a very high accuracy in the waveforms is not needed, the multiblock structure was not used and the authors have used a fourth-order finite-difference operator with a third-order Implicit-Explicit Runge–Kutta integration in time (see Section 2.3).

### 2.2. The Maxwell Equations

- More specifically, these evolution equations describe damped wave equations and have the effect of dynamically controlling the possible growth of the violations of the constraints and of propagating them away from the problematic regions of the computational domain where they are produced.
- The charge density q can be computed either through the evolution (Equation (13)) or by inverting the constraint (Equation (7)).

### 2.3. Numerical Treatment of the Force-free Conditions

- As noted before, within an FF approximation the stressenergy tensor is dominated by the EM part and the contribution coming from the matter can be considered zero.
- From the numerical point of view, specific strategies must be adopted in order to enforce the FF constraints expressed by Equations (25) and (26).
- The authors strategy, however, differs from both the previous ones and follows the same philosophy behind the choice of the driver defined by Equation (29).

### 3. ANALYSIS OF RADIATED QUANTITIES

- The calculation of the EM and gravitational radiation generated during the inspiral, merger, and ringdown is an important aspect of this work as it allows us to measure the amount correlation between the two forms of radiation.
- Any measure of these quantities in the strong-field region is therefore subject to ambiguity and risks producing misleading results.
- The term Φ0 in Equation (36) has been maintained (it disappears at null infinity) to account for the possible presence of an ingoing component in the radiation at finite distances.
- In particular, Equation (36) shows that the net flux is obtained by adding (with the appropriate sign) the respective contributions of the outgoing and ingoing fluxes.
- The second approach that the authors have followed for the computation of the emitted luminosity is the evaluation of the flux of the Poynting vector across a 2-sphere at large distances in terms of the more familiar 3+1 fields Ei and Bi in Equation (19).

### 4. ASTROPHYSICAL SETUP AND INITIAL DATA

- More specifically, the authors consider the astrophysical conditions during and after the merger of two supermassive BHs, each of which is surrounded by an accretion disk.
- During this phase, the binary evolves on the dynamical viscous timescale τd of the circumbinary accretion disk, which is regulated by the ability of the disk to transport its angular momentum outward (either via shear viscosity or magnetically mediated instabilities).
- Because τGW and τd have a very different scaling with D, more specifically, τGW ∼ D4, while τd ∼ D2, at a certain time the timescale τGW becomes smaller than τd.
- This represents the astrophysical scenario in which their simple model is then built.

### 4.1. Initial Data and Grid Setup

- The authors construct consistent BH initial data via the “puncture” method as described in Ansorg et al. (2004).
- The authors use these two configurations to best isolate the effects due to the binary orbital motion from those related to the spins of the two BHs.
- The authors note that similar initial data were considered by Koppitz et al. (2007), Pollney et al. (2007), and Rezzolla et al. (2008a, 2008b, 2008c) but they have recalculated them here using a higher resolution and improved initial orbital parameters.
- Shorter, higher-resolution simulations have also been carried out to perform consistency checks.
- Note that the initial Arnowitt, Deser, Misner (ADM) mass of the spacetime is not exactly 1 due to the binding energy of the BHs.

### 5. ACCURATE FORCE-FREE ENFORCEMENT

- As mentioned in Section 2.3, several different approaches are possible to enforce the FF conditions (25) and (26) in the plasma.
- In fact, since this approach acts “by hand” on the EM fields and converts them to values which would yield an FF regime, one is guaranteed that the constraints (25), (26), and (30) are satisfied.
- Denotes the second step of the “discrete” approach of Palenzuela et al. (2010a), which amounts to the modification of the electric field according to Equation (32), also known as 3. discrete2.
- On the other hand, the corresponding currents when the prescriptions discrete1 and discrete2 are used (right column) do not show evident signs of descending currents and, rather, they show unphysical features around the BH and discontinuities along the ∼ ±45◦ diagonals when seen in the (x, z) plane.
- As a final remark the authors also note that their prescriptions (29) and (33) also provide a saving in computational costs.

### OF BBH MERGERS

- After having discussed the details of their implementation of the FF conditions and having shown its higher accuracy with respect to alternative suggestions in the literature, in the following the authors concentrate their discussion on the FF electrodynamics accompanying the inspiral and merger of BH binaries.
- In particular, the authors will discuss the subtleties that emerge with the subtraction of the background radiation, the spatial distribution of the charge density, the EM and GW zones, and the scaling of the EM luminosity with frequency.

### 6.1. Subtraction of Background Radiation

- Hence, a proper identification of this background radiation is essential for the correct measure of the emitted luminosity and to characterize its properties.
- The generic expression (37) for the EM luminosity can be evaluated in combination with Equation (39), that is, by setting as background values those of the Newman-Penrose scalars Φ2 and Φ0 at the initial time.
- As a result, the background choice (39) represents by far the most convenient one.
- Note that the only modes that have a regular time modulation, and are therefore radiative, are (Φ2)22 and (Φ0)22, while the real parts of the (Φ2)20 and (Φ0)20 are essentially constant in time, indicating that these are not radiative modes, and could represent a way to measure the background radiation.

### 6.2. Properties of the EM Luminosity

- Evolution, measured in hours before the merger, of the luminosities as computed with expression (37) and either the prescriptions (39) or (40) for the background subtraction.
- More specifically, the thick lines refer to the total luminosity, while the thin ones to the luminosity in a polar cap of 5◦ semi-opening angle, measured using either expression (39) (red solid line), expression (40) (blue dotted line), or through the expression in terms of the Poynting vector (41) (black dashed line).
- A few comments should be reserved about the different spatial distributions of the EM fluxes that come with the different prescriptions for the subtraction of the background radiation and that are erased when computing the luminosities as integral quantities.
- The differences introduced by the spin are reported in the right panel of Figure 4, which refers to the binary s6 and thus with BHs having a dimensionless spin of J/M2 0.6.
- Less obvious, however, is the fact that the wave zones can be different whether one is considering the gravitational or the EM radiation, with the latter starting at considerably larger distances than the former.

### 6.3. Frequency Scaling

- As remarked already in Paper I, an accurate measure of the evolution of the collimated and non-collimated contributions of the emitted energies is crucial to predict the properties of the system when the two BHs are widely separated.
- The scaling ∼Ω2/3 is clearly incompatible with their data and the authors suspect the accelerated motion of the BHs to be behind this difference and longer simulations will be useful to draw robust conclusions.
- The right panel of Figure 6, which is the same as the left one but where the authors extrapolate the scaling back in frequency.
- The authors rough estimate is therefore that the collimated emission will be larger than the diffused one at an orbital frequency Ω = (1/2)ΩGW 3.2 × 10−5.
- While this is an exciting possibility, the authors should also bear in mind that, when extrapolated back to the time when it becomes dominant, the collimated emission has also decreased by almost one order of magnitude and to luminosities that are only of the order of ∼1042 erg s−1.

### 6.4. Charge-density Distribution

- Providing information which is complementary to the one already presented by Palenzuela et al. (2010b, 2010c) and Neilsen et al. (2011).the authors.
- To further limit the amount of information that can be extracted directly from their simulation is the fact that an FF code does not allow for an unambiguous calculation of the plasma velocity, which can only be estimated a posteriori based on a certain number of assumptions.
- Be appreciated from the first two columns of Figure 7, while the second contribution is the only one responsible for the charge distribution in the last column.
- Note that since they both refer to isolated spinning BHs (although with different spins), the right column of Figure 7 should be compared with the right column of Figure 2, which shows instead the electric currents.
- Within an FF approach, the consequences of this regular and alternate distribution of positive and negative charges, it is clear that it can lead to rather intriguing particle acceleration processes along the surfaces separating regions of different charges.

### 7. PROSPECTS AND CONCLUSIONS

- Assessing the detectability of the EM emission from merging BH binaries is much more than an academic exercise.
- Nevertheless, relying on a number of assumptions with varying degree of realism, several investigations have been recently carried out to investigate the properties of these EM counterparts either during the stages that precede the merger or in those following it.
- The authors have therefore provided the first quantitative estimates of the scaling of the EM emission with frequency and shown that the diffused part has a dependence that is very close to the one exhibited by the GW luminosity and therefore of the type Lnon-coll EM ≈ Ω10/3−8/3.
- Unfortunately, however, their use of an FF condition (and their ability to maintain it essentially to machine precision) prevents us from producing such electric fields and hence the corresponding accelerations.
- Even in the optimistic case in which the majority of the Poynting flux is converted into radio emission via synchrotron processes, the EM radiation (either collimated or diffused) will eventually exit the evacuated central region around the binary and penetrate in the ambient medium.

Did you find this useful? Give us your feedback

...read more

##### Citations

150 citations

112 citations

### Cites background from "Accurate Simulations of Binary Blac..."

...…last few years has seen work on force-free magnetospheres of binary black hole and neutron star systems (e.g., Palenzuela, Lehner & Liebling 2010; Alic et al. 2012; Palenzuela et al. 2013; Paschalidis, Etienne & Shapiro 2013), motivated in part by the possibility of observing electromagnetic…...

[...]

^{1}, University of Göttingen

^{2}, Harvard University

^{3}, Aarhus University

^{4}, University of Birmingham

^{5}, University of Sydney

^{6}, Search for extraterrestrial intelligence

^{7}, University of Pennsylvania

^{8}, NASA Exoplanet Science Institute

^{9}, University of Copenhagen

^{10}, Max Planck Society

^{11}, Ames Research Center

^{12}, Yale University

^{13}, Iowa State University

^{14}

74 citations

68 citations

### Cites background or methods from "Accurate Simulations of Binary Blac..."

...In a similar way, the method for linear relaxation terms described in subsection 4.2.1 can be generically used for non-linear algebraic currents with the condition that the non-linear terms are evaluated at previous time steps, as it was considered in (Alic et al. 2012)....

[...]

...…limit can also be achieved by considering an effective anisotropic conductivity with a generic form given by (Komissarov 2004; Moesta et al. 2012; Alic et al. 2012) Ji = qvid + σ‖ B2 [ (EkBk)B i +χ(E2−B2)E i ] , (56) where σ‖ is the (anisotropic) conductivity along the magnetic field lines....

[...]

63 citations

### Cites background from "Accurate Simulations of Binary Blac..."

...…al. 2011; Yunes et al. 2011a) or with the aid of electromagnetic observations (Kocsis et al. 2006; Kocsis & Loeb 2008; Kocsis et al. 2012; Giacomazzo et al. 2012; Farris et al. 2012; Noble et al. 2012; Bode et al. 2012; Alic et al. 2012; Gold et al. 2014b; McKernan et al. 2013; Farris et al. 2015)....

[...]

##### References

4,007 citations

### "Accurate Simulations of Binary Blac..." refers background in this paper

...Finally, the local EM flux from the jets can in principle be enhanced if the BHs are spinning and, indeed, within a Blandford–Znajek process one expects that the luminosity from the jets increases quadratically with the spin of the BH (Blandford & Znajek 1977; Palenzuela et al. 2010b)....

[...]

...2009, Astrophysics and Space Science Library, Vol. 357, Neutron Stars and Pulsars Beskin, V. S. 1997, Soviet Physics Uspekhi, 40, 659 Binétruy, P., Bohé, A., Caprini, C., & Dufaux, J.-F. 2012, arXiv:1201.0983 Blandford, R. D., & Znajek, R. L. 1977, Mon....

[...]

...However, a series of more recent numerical simulations in which currents and charges are taken into account, have suggested the intriguing possibility that a mechanism similar to the original one proposed by Blandford and Znajek may be activated in the case of binaries (Palenzuela et al. 2009, 2010; Palenzuela et al. 2010b,a; Moesta et al. 2012); note that Palenzuela et al. (2010b,a); Moesta et al. (2012) also make use of a ar X iv :1 20 4....

[...]

...These physical conditions are indeed similar to those considered in the seminal investigations of BH electrodynamics of Blandford and Znajek (Blandford & Znajek 1977), who addressed the question of whether the rotational energy of an isolated BH can be extracted efficiently by a magnetic field....

[...]

...In particular, the Blandford–Znajek mechanism is likely to be valid under rather general conditions, namely even if stationarity and axisymmetry are relaxed and even if a non-spinning BH is simply boosted through a uniform magnetic field....

[...]

1,568 citations

1,349 citations

### "Accurate Simulations of Binary Blac..." refers methods in this paper

...The first one uses the Newman-Penrose scalars Φ0 (for the ingoing EM radiation) and Φ2 (for the outgoing EM radiation), defined using the same tetrad (Teukolsky 1973): Φ0 ≡ Fµν lνmµ , Φ2 ≡ Fµνmµnν ....

[...]

...Following Teukolsky (1973), we compute the total EM luminosity as a surface integral across a 2-sphere at a large distance: L EM = lim r→∞ 1 2π ∫ r2 ( |Φ2|2 − |Φ0|2 ) dΩ , (36) which results straightforwardly from the integration of the component of EM stress-energy tensor (4) along the timelike…...

[...]

1,053 citations

### "Accurate Simulations of Binary Blac..." refers background in this paper

...…Maxwell equations (2) and (3) is referred to as “extended” because it incorporates the so-called divergencecleaning approach, originally presented in Dedner et al. (2002) in flat spacetime, and which amounts to introducing two additional scalar fields, Ψ and Φ, that propagate away the deviations…...

[...]

913 citations