Charged and rotating multi-black holes in an external gravitational field

We construct analytical and regular solutions in four-dimensional General Relativity which represent multi-black hole systems immersed in external gravitational field configurations. The external field background is composed by an infinite multipolar expansion, which allows to regularise the conical singularities of an array of collinear static black holes. A stationary rotating generalisation is achieved by adding independent angular momenta and NUT parameters to each source of the binary configuration. Moreover, a charged extension of the binary black hole system at equilibrium is generated. Finally, we show that the binary Majumdar-Papapetrou solution is consistently recovered in the vanishing external field limit. All of these solutions reach an equilibrium state due to the external gravitational field only, avoiding in this way the presence of any string or strut defect.


Introduction
Multi-black holes solutions are intriguing astrophysical objects both from the theoretical and from the phenomenological point of view. On the theoretical side, these solutions disclose the non-linear nature of General Relativity and represent an important playground in which testing the laws of black hole mechanics. On the experimental side, the recent remarkable observations of gravitational waves [1] heavily rely on the interactions between two black holes in a binary system: thus an analytical description of such a spacetime is of utmost relevance for the interpretation of the measurements.
Of course one of the main obstacle in modelling a stationary multi-gravitational sources system is to provide a mechanism to balance the gravitational attraction of the bodies. Otherwise the system naturally tends to collapse. Usually the equilibrium is granted by the introduction of cosmological struts which prop up the gravitating bodies, but these one-dimensional objects must be constituted by matter that violates physically reasonable energy conditions. Alternatively in some cases cosmic strings of infinite length can be used to support the gravitational collapse. In both cases these objects are symptomatic for conical singularities which plague the spacetime. Other objects, such as the Misner string, which occur in the presence of the gravimagnetic parameter, also known as NUT charge, are sometimes used to prevent the merging [2,3], but these bring in also harmful issues such as closed timelike curves. Instead our main objective is to remain as close to phenomenology as possible, therefore we will try to avoid physical structures which are outside realistic experimental or observational range.
We are interested in stationary and axisymmetric models basically for mathematical reasons: the Einstein equations enjoy some integrability properties which can be exploited to generate highly non-trivial solutions without resolving the equation of motion. To take advantage of the symmetries underling axisymmetric and stationary General Relativity we will mainly use the inverse scattering technique invented by Belinski and Zakharov [4], which will be briefly summarised in section 2. Also the coupling of the Einstein theory to the Maxwell electromagnetic field preserves the integrability of the axisymmetric and stationary system, therefore we will explore the possibility of including electromagnetic charge to some multi-black hole configuration, even though it seems quite difficult that these charged objects can exist spontaneously in Nature, because the distribution of matter in the Universe is neutral on large scale.
Basically the only regular multi-black holes configurations known so far were the extremal ones where the gravitational attraction is balanced by the electromagnetic force, such as in the Majumdar-Papapetrou solution [5,6] or the magnetised dihole [7] 1 . Only very recently a new physical mechanism able to maintain the equilibrium was presented in the realm of vacuum General Relativity [9]. It consists in the introduction of an external gravitational field endowed with a number of integration constants related to the field multipolar expansion. In [9] a simple model was studied, which represents a system of two black holes immersed in a dipole-quadrupole external gravitational field. Here we extend such solution in several ways: firstly, we consider an arbitrary number of black holes immersed in an external field described by an infinite number of multipole momenta. Then, in section 5 and 6, we build the rotating and charged generalisation of [9].
External gravitational fields represent a natural setting for multi-black holes systems, as recent gravitational waves detection proceeding from the center of galaxies confirms. In particular it has been shown in [10] that the multipolar gravitational field, we will deem here, can be produced by a distribution of matter such as thin disks or rings, typical shape of gravitational objects such as galaxies or nebulae. Anyway the solutions considered in this article will be pure vacuum solution, without any energy-momentum tensor. In principle a distribution of matter might be possibly considered very far away from the black holes. In this sense our solutions can be interpreted as local models for binary or multi-black hole configurations. In a certain sense the metrics presented here have to be considered as the gravitational analogous of the stationary 1 C-metrics embedded in an external magnetic field, found by Ernst [8], can represent also a couple of black holes, when considering their maximal extension. However these objects are two copies of the same black hole, always causally disconnected, and not interacting with each other, thus not really physically realistic. black holes in Melvin magnetic universe. In fact, also in this latter case, the solution is a pure electrovacuum solution with no definite sources for the electromagnetic external field, therefore their feasibility remains in the proximity domain of the black bodies.
Single black holes in external multipolar gravitational field have been pioneered by Doroshkevich, Zel'dovich and Novikov [11], later studied by Chandrasekhar [12], Geroch and Hartle [13]. These solutions are known in the literature as deformed black holes. The novelty of our proposal is to take advantage of the external field to sustain the black bodies and prevent their collision. From a mathematical point of view this means that the solution can be regularised from conical singularities that usually affect multi-black hole metrics.
The plan of the paper is the following: we begin by reviewing the inverse scattering method, which is the leading method for the generation of our solutions, in section 2. Then, we describe the main properties of the background metric which produces the "distortion" on the black holes sources in section 3: there we discuss both internal and external deformations from the point of view of the gravitational multipoles, even if in the following we make use of the external deformations only. Section 4 is devoted to the array of collinear and static black holes immersed in the external field, while in sections 5 and 6 we present the rotating and charged generalisations, respectively, of the binary system found in [9]. Finally, we summarise our work and present some conclusions.

Inverse scattering method
We make use of the inverse scattering method [4,14,15] to superimpose black holes on top of the background spacetime containing an external gravitational field. Thus, we begin by discussing the key features of the solution generation technique and subsequently by constructing the quantities apt to generate the desired spacetime.
The inverse scattering method relies on the integrability of the Einstein equations for the class of stationary and axisymmetric spacetimes, which can be described by the general metric in Weyl coordinates [16] where a, b = 0, 1 and x 0 = t, x 1 = φ. The spacetime (2.1) possesses two commuting Killing vectors proportional to ∂ t and ∂ φ . We assume that the coordinate ρ is chosen such that det g = −ρ 2 . The vacuum Einstein equations R µν = 0 can be equivalently written as where U = ρg ,ρ g −1 and V = ρg ,z g −1 are two 2 × 2 matrices. We see that, by solving equation (2.2a) for g, one is able to find the function f in quadratures from equations (2.2b) and (2.2c). Thus, the problem of solving the vacuum Einstein equations is reduced to the problem of finding the matrix g.
One can show that the integrability condition for the Einstein equations (2.2) is equivalent to the linear equations for the generating matrix Ψ(ρ, z, λ), where the commuting differential operators D 1 and D 2 are given by and λ is a complex spectral parameter.
In the inverse scattering method we prepare a seed solution (g 0 , f 0 ), and then find a generating matrix Ψ 0 which satisfies the linear eigenvalue equations (2.3). Given such a Ψ 0 , we introduce the functions where w k are arbitrary (complex) constants, called poles. µ k andμ k are called solitons and anti-solitons, respectively, and they satisfy µ kμk = −ρ 2 . We will make us of the solitons only, in the following. We associate a 2-dimensional vector (called BZ vector) to each (anti-)soliton with arbitrary constants m (k) 0 a . Given the matrix a new metric is obtained by adding 2N solitons to (g 0 , f 0 ) as c (g 0 ) ca and C f is an arbitrary constant. The new metric (2.8) satisfies by construction the Einstein equations (2.2) and is such that det g = −ρ 2 .

Background metric: internal and external multipoles
We discuss the general solution to the Einstein equations in vacuum, which contemplates both the internal deformations of the source and the contributions which come from matter far outside the source. This general solution finds its roots in the pioneeristic work of Erez and Rosen [17], and it was lately discussed and expanded in [11,13,12] to include the deformations due to an external gravitational field.
The general solution for the Weyl metric is given, following the conventions in [18], by ψ = ∞ n=1 a n r n+1 + b n r n P n , (n + 1)(p + 1)a n a p (n + p + 2)r n+p+2 (P n+1 P p+1 − P n P p ) + npb n b p r n+p n + p (P n P p − P n−1 P p−1 ) , (3.2b) where r := ρ 2 + z 2 defines the asymptotic radial coordinate and P n = P n (z/r) is the n-th Legendre polynomial. The real constants a n describe the deformations of the source, while the real parameters b n describe the external static gravitational field.
We observe that the "internal" part a n , which is related to the deformations of the source, is asymptotically flat: this seems to contradict Israel's theorem [19], which states that the only regular and static spacetime in vacuum is the Schwarzschild black hole. Actually, the internal deformations lead to curvature singularities not covered by a horizon, in agreement with the theorem. Because of this feature, in the following sections we will discard the internal contributions and focus on the external ones only.
On the converse, the "external" part b n is not asymptotically flat. This is in agreement with the physical interpretation: this part of the metric represents an external gravitational field generated by a distribution of matter located at infinity. This interpretation parallels the Melvin spacetime one, and in fact there is no stress-energy tensor here to model the matter responsible of the field. The asymptotia is not flat because at infinity there is the source matter and the spacetime is not empty. Actually it can be shown that the curvature invariants, such as the Kretschmann scalar, can grow indefinitely at large distances 2 in some directions. In this sense, a spacetime containing such an external gravitational field should be considered local, in the sense that the description is meaningful in the neighborhood of the black holes that one embeds in this background. In this regard these metric are not different with respect to the usual single distorted black hole studied in the literature [18], [20]. To have a completely physical solution, one should match the black holes immersed in the external field with an appropriate distribution of matter (such as galaxies), which possibly asymptotes Minkowski spacetime at infinity. A model for matter content consistent with the multipolar background expansion treated here and based on a thin ring distribution can be found in [10] 3 .
Accelerating generalisations of the multipolar gravitational background studied in this section naturally are endowed with a Killing horizon of Rindler type. These accelerating horizon are typically located in between the black hole sources and the far region [21].
The constants a n and b n are related to the multipole momenta of the spacetime Q n . The relativistic definition of the multipole momenta was given by Geroch [22], and lately refined by Hansen [23]. That definition applies for a stationary, axisymmetric and asymptotically flat spacetime: the modified Ernst potential associated to the spacetime is expanded at infinity, and the first coefficients of the expansion correspond to the multipole momenta [24]. Clearly, the internal deformations asymptote a flat spacetime, while the external ones do not.
In the following, we compute the multipole momenta associated to the above deformations, in order to clarify the intepretation of a n and b n . In particular, we calculate the internal multipoles contributions by means of the standard definition, and then we propose a new approach for determining the multipoles associated to external deformations. We shall see that the proposed definition is analogous to the usual one and gives a consistent interpretation of the external field parameters.

Internal multipoles
Let us start with the internal deformations. If we turn off the b n , then we are left with We define ξ as a function of the Ernst potential (cf. Appendix B) that in our case reads We want to expand the above expression around infinity: following [24], we bring infinity to a finite point by definingρ = ρ/(ρ 2 + z 2 ),z = z/(ρ 2 + z 2 ), and by conformally rescaling ξ int : (3.6) One can prove thatξ is uniquely determined by its value on the z-axis: since infinity corresponds toρ =z = 0 in the new coordinates, then we expandξ aroundz for ρ = 0 where M j are the expansion coefficients. The multipole momenta are completely determined by the coefficients M j , and in particular one can show (see [24] and references therein) that the first four coefficients M 0 , . . . , M 3 are exactly equal to the first four multipole momenta Q int 0 , . . . , Q int 3 . Thus, for the internal deformations, we find We see that there is a direct correspondence between the coefficients a n and the multipole momenta, at least at the first orders. Q int 0 is the monopole term, which is zero since no source (e.g. no black hole) is present. Q int 1 is the dipole, Q int 2 is the quadrupole and Q int 3 is the octupole moment. The subsequent multipole momenta can be still computed from the expansion (3.7) by means of a recursive algorithm, but now they will be non-trivial combinations of the coefficients M j [24]. For instance, the 16-pole is given by

External multipoles
The above construction can not work for the external deformations: the asymptotia is not flat, and infinity is the place where the sources of the external field are thought to be, hence it does not make sense to expand there. On the converse, it is meaningful to detect the effects of the deformation in the origin of the Weyl coordinates: thus, by paralleling the Geroch-Hansen treatise, we propose to expand the modified Ernst potential in the origin of the cylindrical coordinates. Now we consider only the external deformations with modified Ernst potential According to the above discussion, we assume that ξ ext is completely determined on the z-axis, so we expand around the origin for ρ = 0 11) where N j are the expansion coefficients. We define the first four multipole momenta Q ext 0 , . . . , Q ext 3 as the coefficients N 0 , . . . , N 3 , which are equal to Again, the monopole moment is zero because of the absence of a source. We observe, contrary to (3.8), that the octupole Q ext 3 is given by a non-trivial mixing of the constants b 1 and b 3 . This definition occurs only for the first momenta: a definition which takes into account higher momenta can be achieved by generalising the procedure in [24]. For the time being, we content ourselves by proposing the following interpretation: the coefficients b n are related to the multipole momenta generated by the external gravitational field, similarly to what happens for the internal deformations. Then, the presence of the external gravitational field affects the momenta of a black hole source immersed in it.

Seed for the inverse scattering construction
We are interested in the external deformations only, hence we consider a n = 0 hereafter, and focus only on the contributions from b n . We express the external gravitational field metric in a form which is suitable for the inverse scattering procedure, i.e.
The parameters b n are related the multipole momenta of the external field, as explained above. Metric (3.13) represents a generic static and axisymmetric gravitational field.
Since we want to add black holes on the background represented by the gravitational field, we take the metric (3.13) as a seed for the inverse scattering procedure. Following the discussion in section 2, we need the generating matrix Ψ 0 , which serves as starting point to build the multi-black hole solution. The function which satisfies equations where (3.15) Now we can construct the BZ vectors (2.6): we parametrise m are constants that will be eventually related to the physical parameters of the solution. The BZ vectors are thus Depending on the value of C (k) 0 and C (k) 1 , the spacetime will be static or stationary, as we will see in the following.

Array of collinear static black holes in an external gravitational field
We now proceed to the generalisation of the Israel-Khan solution [26], which represents an array of collinear Schwarzschild black holes. The Israel-Khan metric is plagued by the presence of conical singularities which can not be removed by a fine tuning of the physical parameters 4 . On the converse, we will see that the external gravitational field will furnish the force necessary to achieve the complete equilibrium among the black holes. Given the seed metric (3.13) and the generating matrix (3.14), we construct a new solution by adding 2N solitons with constants This choice guarantees a diagonal, and hence static, metric. Each couple of solitons adds a black hole, then the addition of 2N solitons gives rise to a spacetime containing N black holes, whose metric is Metric (4.2) is, by construction, a solution of the vacuum Einstein equations (2.2), and it represents a collection of N Schwarzschild black holes, aligned along the z-axis, and immersed in the external gravitational field (3.13).
We limit ourselves to the case of real poles w k , since it represents the physically most relevant situation. These constants are chosen with ordering w 1 < w 2 < · · · < w 2N −1 < w 2N and with parametrisation The constants m k represent the black hole mass parameters, while z k are the black hole positions on the z-axis.
The black hole horizons correspond to the regions w 2k−1 < z < w 2k (k = 1, . . . , N ), while the complementary regions are affected by the presence of conical singularities, as happens for the Israel-Khan solution (cf. Fig. 1). Differently from that case, our solution can be regularised, as we will show in the next subsection.
It might be useful to write the metric (4.2) in the canonical Weyl form (3.1), in order to provide a comparison with the literature. By writing down explicitly the solitons (according to (2.5)) and the values for w k , one finds the Weyl functions where and we defined R ± k := ρ 2 + (z − z k ± m k ) 2 . Actually, this form of the metric corresponds to the choice for C f made below (cf. eq. (4.8)). From this form, it is evident that the Israel-Khan solution [26] is recovered simply by putting b n = 0.

Conical singularities and regularisation
The infinite multipole momenta b n allow to regularise the metric, i.e. to remove all the conical singularities, by tuning their values. More precisely, given N black holes, there will be exactly N + 1 conical singularities: two cosmic strings, one rear the first black hole and one ahead the last black hole, and N − 1 struts located between the N black holes. Hence, one will need at least N + 1 parameters to fix the singularities. We want to remove both the struts and the strings, since they are unphysical and lead to energy issues, as is shown in appendix A.
The manifold exhibits angular defects when the ratio between the length and the radius of small circles around the z-axis is different from 2π. Working in Weyl coordinates, a small circle around the z-axis has radius R = √ g zz ρ and length L = 2π √ g φφ [9]. The regularity condition corresponds then to L/(2πR) → 1 as ρ → 0. It is easy to prove that, for the static and axisymmetric metrics of the class (2.1), the above condition is equivalent to P ≡ f g tt → 1 as ρ → 0. We choose for convenience the gauge parameter C f as The quantity P = f g tt is equal to between the k-th and (k + 1)-th black holes (i.e. w 2k < z < w 2k+1 ), for 1 ≤ k < N . In the region z < w 1 we find while for z > w 2N we simply have P N = 1 , (4.11) thanks to our choice of C f . These expressions are the natural generalisations of the conical singularities for the Israel-Khan metric [28]. The above expressions provide a system of equations P k = 1 for 0 ≤ k < N , which can be solved, e.g., for the parameters b 1 , . . . , b N , with the result of fixing all the conical singularities. Hence the solution can be made completely regular outside the black hole horizons.
As a pictorial example we show in a figure the geometry of the black holes array in the triple hole configuration. As can be appreciated from Fig. 2, the surfaces of the event horizons are regular because, thanks to the procedure explained above, the spacetime is everywhere devoid from conical defects. At this purpose some integration constants, equal to the number of the spacelike rods of the rod-diagram, have to be fixed according to Eqs. (4.8)-(4.10).

Smarr law
We investigate the Smarr law for spacetime (4.2): to this end, we compute the total mass of the spacetime and the entropy and temperature of the black holes.
The mass is easily found by means of the Komar-Tomimatsu [29,30] integral. The result for the k-th black hole (i.e. the black hole in the interval w 2k−1 < z < w 2k ) is where α is a constant which takes into account the proper normalisation of the timelike Killing vector, generator of the horizon, ξ = α∂ t associated to (4.2). The entropy of a black hole is related to the area as S k = A k /4, hence where 14) The product f g φφ is independent of z in the limit ρ → 0, and that was crucial in the derivation of (4.13). Finally, the temperature is found via the Wick-rotated metric, and the result is It is easily shown, by using (4.12), (4.13) and (4.15), that the Smarr law is satisfied: We notice that the explicit value of α is not needed for (4.16) to work, while it is relevant in the study of the thermodynamics [9].

Distorted Schwarzschild black hole
The simplest non-trivial example we can consider for the complete external multipolar expansion from the general solution (4.2), is clearly the single black hole configuration for N = 1. In that case the functions that appear in the Weyl static metric (3.1) take the form This spacetime represents a static black hole embedded in an external gravitational field. The limit to the Schwarzschild metric is clear: it is obtained just by switching off all the multipoles b n = 0 ∀n. In order to recover the standard Schwarzschild metric in spherical coordinates the following transformation is needed see also [18]. However the form we are writing here is more general because, thanks to the extra parameter z 1 , allows to place the black hole in any point of the z-axis. In the absence of the external field the location of the black hole is irrelevant because the solution is symmetric under a finite shift of z 1 . But, when the external gravitational field is not null, a boost along the z-axis is significant, since the relative position of the black hole with respect to the external multipolar sources has some non-trivial effects on the geometry and on the physics of the black hole.
In particular the translation along the z-axis affects the event horizon shape. In fact, as can be evaluated by computing the equatorial and polar circles around the event horizon, it is possible to understand as the horizon surface stretches or contracts depending on the position of the black hole and the values of the external parameters. Some pictorial examples of the black hole horizon deformation for different external gravitational backgrounds are given in Fig. 3.
Moreover, the equilibrium constraint which removes the conical defects can be loosen with respect to the one usually found in the literature [18]. The imposition usually requested in the literature, i.e ∞ n=1 b 2n+1 = 0, is not fundamental when the black hole can be adjusted coherently with the external gravitational field. In fact when z 1 = 0 the two regularising constraints we have to impose to avoid conical singularities become, as seen in section 4.1, When b n = 0 for n > 2, we obtain a special subcase of the solution [9], in the limit where one of the black hole vanishes or where the two horizon rods merge, remaining only with a single black hole configuration 5 . Fig. 4 in this section refer, for simplicity, to this truncated expansion of the external field. A qualitatively analogous behaviour of the black hole horizon occurs in the full multipolar expansion.

Two Kerr-NUT black holes in an external gravitational field
We briefly discuss the rotating and NUTty generalisation of the binary black hole system immersed in an external gravitational field. The generalisation to an array of rotating and NUT black holes is straightforward, at least conceptually. We start again with the seed metric (3.13) and add four solitons (i.e. two black holes). However we make a different choice for the BZ constants: we choose [31] Here m i are the mass parameters, a i are the angular momenta and n i are the NUT parameters. We have also defined σ 2 i ≡ m 2 i − a 2 i + n 2 i . The poles w k are naturally defined as where z i are the positions of the black holes. The resulting metric is computed by following the inverse scattering method of section 2, but it is quite involved and we will not write it explicitly here. The simplest form of the metric is achieved by using the bipolar coordinates The metric regularisation from angular defects on the symmetry axis, can be pursuit as in the static case above, by tuning one physical parameter of the solution for each N + 1 spacelike rod.
In the simplest case of a dipole-quadrupole configuration, such a binary system is characterised by ten parameters: {m 1 , m 2 , a 1 , a 2 , n 1 , n 2 , z 1 , z 2 , b 1 , b 2 }. Contrary to the standard double-Kerr case, in which there is no external field, the relevant physical quantity is not the distance between the holes z 2 − z 1 . This happens because the external field does not act uniformly on the z-axis, hence the translation invariance along the axis is broken. This means that the positions of the black holes are two independent parameters.
However, when the external field parameters are set to zero, b 1 = b 2 = 0, one recovers the usual double-Kerr-NUT solution [31] (see [32] for a recent account on the equilibrium configurations). Moreover, one can see that the positions z 1 and z 2 can be reabsorbed into a single parameter l = z 2 − z 1 . Obviously, in such a limit the spacetime can not be regularised to give a physical solution, and the conical singularities can not be avoided. The spacetime is then affected by the presence of struts or cosmic strings, unless one admits "naked singularity-black hole" or "naked singularity-naked singularity" configurations.

Two Reissner-Nordström black holes in an external gravitational field
A natural generalisation of the solutions presented in the previous sections contemplates the addition of the electric charge, thus we will construct a multi-Reissner-Nordström solution immersed in an external gravitational field. Such a solution is the multi-black hole version of the metric presented in [18].

The charging transformation
There are several procedures to extend a stationary and axisymmetric solution of general relativity to support a Maxwell electromagnetic field. This methods allows, for instance, to generate the Reissner-Nordström spacetime from the Schwarzschild black hole solution, such as the Harrison [33] or the Kramer-Neugebauer [34] transformations. We want to present here perhaps the simplest version of this charging transformation, which maps a given static and axisymmetric vacuum solution to another static and axisymmetric electrovacuum spacetime, typically adding monopole electric charge 6 . Let us consider the most general static and axisymmetric metric for Einstein-Maxwell theory, the Weyl metric (3.1). Suppose the electromagnetic vector potential associated to this metric, compatible with the symmetries of this system, is null, i.e. A µ (ρ, z) = 0. Then an electrically charged solution can be generated by the following transformation on the ψ function of the metric which is be supported by an electric field given bŷ The continuous parameter ζ can be considered real and it is related to the electric charge of the spacetime.
In appendix B we derive this transformation from the Kramer-Neugebauer one and we show how the latter is contained in the Harrison transformation up to some gauge transformations. There we also provide a simple example of application of the charging transformation (6.1), (6.2), where the electrically charged black hole solution of Reissner-Nordström is generated from the Schwarzschild metric.

Generating two distorted Reissner-Nordström black holes
Now we want to charge a multi-Schwarzschild solution embedded in an external gravitational field. To keep the model as simple as possible, without constraining the physical parameters of the black hole, we consider, as a seed, the double-black hole spacetime immersed in an external gravitational field possessing dipole and quadrupole moments only. Note that this choice is done only for simplicity, but it could be chosen any external gravitational expansion endowed with multipoles of any order; likewise the charging method allows one to deal easily with an arbitrary number of sources. However we will act with the charging transformation (6.1), (6.2) on the solution (4.2) for N = 2 and with only the coefficients b 1 and b 2 different from zero. The resulting seed metric, which coincides with the one described in [9], is determined by the two functions of the static Weyl metric (3.1) where W ij = ρ 2 + µ i µ j and Y ij = (µ i − µ j ) 2 , while the solitons are defined as in (2.5). The new solution 7 will maintain the same γ as (6.3b) while transforming ψ according to (6.1). ζ is related to the electric charge of the black holes, therefore the transformed metric will represent a couple of charged black hole embedded in an external gravitational field. In the limit of null external field b n = 0 and vanishing one of the two black hole masses (m 2 = 0 for instance), we exactly recover the Reissner-Nordström solution. Otherwise when the external gravitational field is not present the solution approach the equal mass double charged masses of [35], [36].
We are interested in solutions regular outside the event horizons, therefore we have to consider the quantity P ≡ f g tt to avoid the possible conical singularities of the charged metric. In fact, when P differs from 1, it takes into account the deficit or excess angle along the three regions of the axial axis of symmetry outside the black holes event horizons, i.e. for ρ = 0 and z ∈ (−∞, w 1 ), z ∈ (w 2 , w 3 ) and z ∈ (w 4 , ∞). The solution is made regular from line singularities by imposing the following three constraints on the metric parameters: While this remains formally the same regularisation constraint of the uncharged case, the physical meaning of the parameters is different, as can be easily understood from the single source case treated in appendix B.
In fact the two black hole horizons are located in ρ = 0 and z ∈ (w 1 , w 2 ), z ∈ (w 3 , w 4 ), where Henceforward we consider C f and b n fixed, as in Eqs. (6.4)-(6.6), in order to assure the absence of conical singularities. Note that the proper distance between the two event horizon surfaces converges: It means that the balancing condition is non-trivial and can be realised for a finite separation between the sources.

Charges and Smarr law
The electric charge of each black hole can be computed thanks to the Gauss law [37] Note that non-null results occur only in the regions which define the event horizon of the black holes, as expected. Also note that, since the charging transformation is a one parameter transformation, it adds only an independent electric charge to the system. Thus the free physical parameters of the solution are five: z i , σ i and ζ. Hence the two black holes cannot vary independently their electric charge. More general solutions involving independent electric charge parameters can be built, but with more refined generating techniques such as [38] (or [39]). The mass of the charged black holes can be defined by evaluating, on their respective event horizon, the following integral 8 Considering the mass as a local quantity, i.e. defined close to the horizon, we can fix the gauge degree of freedom in the electric potential as A 0 = 0. Of course other gauge fixings can be pursued, for instance requiring that the electric potential vanish at large radial distances ρ 2 + z 2 → ∞. From the masses and electric charges of the black holes, eqs. (6.9) and (6.10), we can deduce, for any ζ = 1, the value of As expected the the masses and electric charges are not all independent, but they can be expressed just in terms of the parameters M i and ζ. In fact the electric charges can be written, thanks to eqs. (6.9) and (6.10), as The entropy for each black hole is taken as a quarter of the event horizon area which gives 14) The temperature of the event horizons, computed as in the previous section, can be written as From equation (6.12) it is easy to see that when the masses of the two black hole coincides, M 1 = M 2 , also Q 1 = Q 2 . Then it is possible to take ζ as in the single Reissner-Nordström case treated in appendix B In this symmetric case it can be straightforwardly checked that both the temperature and the surface area of the two black holes coincide. Anyway the thermal equilibrium can be reached also for more general sources configurations. Note that the event horizons become extremal in the limit case for the charging transformation, ζ → 1, as in the single black hole case. The Coulomb electric potential Θ evaluated on both the event horizons takes the same value We have now all the ingredients at our disposal to verify the Smarr law, both for the single element 19) and, thus, for the double black hole configuration where we defined Since distorted black holes have a preferred interpretation as local systems, we primarily focused on local quantities, basically defined on the horizon. Nevertheless the above results hold also in the case one considers the presence of the asymptotic Coulomb potential. The gauge freedom encoded in A 0 can be used to put to zero the value of the potential at large distance: when A 0 = ζ sign(b2) , then Θ ∞ = 0. The above results are valid for any α, the normalisation parameter of the Killing vector that generate the horizon ξ = α∂ t . Then, in this context, α can practically regarded as unitary. However, in discussing the first law of black hole thermodynamics, it is necessary to select a particular value for α, as described in [9], for a local point of view based on the assumption that the observer are located close the hole and they have no access to infinity.
This charged black binary configuration is one of the few multi-black hole examples where it is concretely possible to test the second law of black holes thermodynamics, as done for the uncharged case [9] or for the Majumdar-Papapetrou black holes [40]. For instance, when the system is isolated, it is easy to verify that for two configurations, with the same energy and background field, the disjoint state is always less entropic than a collapsed state, which can be thought as the final state: S < S . Anyway, for different boundary conditions this charged case has a rich phase transitions scenario, from adiabatic merging to black hole brimming. However it is outside the scope of this work and will be studied elsewhere.

Majumdar-Papapetrou limit
Inspecting the values of the regularising parameters b n in eqs. (6.5), (6.6) we notice that there is a special case for which they vanish: that happens for w 1 = w 2 and w 3 = w 4 , which, according to eqs. (6.7), (6.11) and (6.12) corresponds to extremality, M i = Q i , or ζ = 1 9 . In that case the standard Minkowskian asymptotics is retrieved.
That is not surprising because the charged double black hole configuration presented in section 6.2 naturally contains a particular subcase of the Majumdar-Papapetrou solution [5,6], the one which describes two identical black holes located along the z-axis. The black holes have to be identical because the charging transformation (6.1), (6.2) adds only one independent charge and the Majumbdar-Papapetrou solution possesses only extremal horizons. In fact the Majumdar-Papapetrou solution is the only configuration, within the double Reissner-Nordström system [35,36], that can reach the equilibrium by balancing the gravitational attraction thanks to the electric repulsion of the sources, without requiring any hyper-extremal event horizon, and thus preserving its black hole interpretation.

Summary and Conclusions
In this article several generalisation of the binary black hole system at equilibrium in an external gravitational field are built, thanks to various solution generating techniques. First of all the complete infinite multipolar expansion for the external field is considered as a background. These multipoles allow us to regularise an arbitrary number of collinear static Schwarzschild black holes at equilibrium. Therefore, thanks to the external gravitational field, we are able to remove all the conical singularities of the Israel-Kahn solution.
The physical quantities computed for this infinite array fulfil the Smarr relation. We notice that in the single deformed Schwarzschild case our solution is more general with respect to similar ones studied in the literature. This is because the one presented here has an extra physical parameter related with the position of the black hole with respect to the multipolar distribution. This novel feature is fundamental in order to have a more general balance constraint, hence an enriched physical description.
Then, with the aid of the inverse scattering method, we show how to add both the angular momentum and the NUT charge to each element of the black hole configuration. In this stationary case the metric describes a deformed multi-Kerr-NUT system.
Moreover we present, thanks to the Kramer-Neugebauer transformation, an electrically charged extension of the binary system at equilibrium. It represents two Reissner-Nordström black holes held in balance by the external gravitational field. Also in this case it is verified that the physical quantities of the double solution satisfy the Smarr relation. It is shown how in the charged case the constraint that assure the absence of conical singularities can be fulfilled also without the need of the external gravitational field because it is sufficient the balance between the gravitational attraction and the electric repulsion between the two black holes to maintain the equilibrium. Therefore this peculiar configuration has a standard flat Minkowskian asymptotics and can not be nothing else than the symmetric Majumdar-Papapetrou bi-hole, which is naturally contained into our charged binary system.
The presence of external gravitational fields such as the one produced by a distribution of matter around the black holes can be used to regularise the angular defects for a wide family of binary black hole systems. In our picture there is no need for struts or cosmic strings which are of uncertain experimental plausibility and of strident theoretical soundness. The multi-sources description proposed here is local, that is we believe have some phenomenological pertinence in the proximity of the black holes, at least till the galactic heavy matter that produce the necessary external gravitational multipolar expansion.
We believe that these new, well behaved, solutions of the Einstein equations can be of some astrophysical interest and also good reference for the Numerical Relativity community. In fact, until now, this latter was the only to posses tools to model and describe multi-black hole systems and to interpret observations. It would be interesting to study the stability of these metrics, even though, we expect that in the long term the final state might be a merging of the binary system. Therefore these models might be mainly useful in some metastable phase of the eventual collision between the black holes.
Several other generalisations could be performed relatively easily thanks to generating techniques both in the realm of General Relativity and in related gravitational theories for which the integrable methods are still applicable, from Brans-Dicke to other scalar tensor theories.

Aknowledgments
This work was supported in part by Conicyt-Beca Chile, in part by MIUR-PRIN contract 2017CC72MK003 and also by INFN.

A Conical singularities and energy conditions
Conical singularities, beyond making the spacetime manifold ill-defined from a mathematical point of view, give also rise to energy issues. In general, such singularities can be interpreted as strings or struts whose energy-momentum tensor has a δ-like nature. We show what are the physical issues that the conical singularities bring in when they are present in the spacetime.
Let us consider Minkowski spacetime with an wedge of angle 2πα artificially removed. By defining C = 1 − α , we can write the metric as where 0 ≤ ϕ < 2π. One can regard this spacetime as a field sourced by a cosmic string or a strut [41], whose non-vanishing energy-momentum tensor components are and δ(x, y) is the two-dimensional delta function depending on the coordinates orthogonal to the z-axis (x, y). µ is interpreted as the tension of the filament source. This result can be generalised to a generic four-dimensional spacetime [42], for which one finds that the Einstein equations G µν = 8πT µν give rise to an energy-momentum tensor C represents again the angular excess/deficit for the azimuthal angle. Let us consider, e.g., a two-black hole spacetime from (4.2) (N = 2): the result (A.4) clearly shows that for µ > 0 (positive tension), the source acts as a string that pull a black hole. This is the behaviour of the conical singularities that one finds at z < w 1 and z > w 4 There are no negative-energy issues in this case, but the string extends to infinity and the δ function gives rise to a divergent energy-momentum tensor.
In the case of µ < 0 (negative tension), the conical singularity in w 2 < z < w 3 acts as a strut which pushes apart the two black holes. The energy density associated to the energy-momentum tensor is negative (i.e. the strut is composed of anti-gravitational matter) and there is again a divergence due to the δ function.

B Harrison and Kramer-Neugebauer charging trasformations
Both the Harrison and the Kramer-Neugebauer [34] transformations are two symmetries of the Ernst equations for the Einstein-Maxwell theory. The Ernst equations are a couple of complex differential equations governing the axisymmetric and stationary electrovacuum of General Relativity [43] Re Greek letters (λ, β, α) represent continuous complex parameters, while latin letters, such as (b, c), label real ones. Some of these transformations, such as (I), (II) and (IV) are pure gauge transformations, hence they could be reabsorbed by a diffeomorphism. The Harrison transformation is (V), while the Kramer-Neugebauer, as defined in [44] to charge the Kerr metric embedded in an external gravitational field is The latter transformation reduces to the one in (6.1), (6.2) for static and uncharged seeds. Since both the Kramer-Neugebauer (B.7) and the Harrison transformation (V), have the same physical effects (they are know to add an electric monopole to an uncharged seed), we have the suspect that they are basically the same transformation, up to gauge transformations. In fact it can be shown that the subsequent composition of transformations (I), (IV) and (V) to a Ernst seed (E 0 , Φ 0 ) gives (V) • (IV) • (I) • E 0 Φ 0 = λλ * E0−β * (β+2λΦ0) 1+α 2 (−2β+α * ββ * −αλλ * E+2λ(αβ * −1)Φ0 β−αββ * +λλ * αE0+Φ0−2λαβ * Φ0 1+α 2 (−2β+α * ββ * −αλλ * E+2λ(αβ * −1)Φ0 . (B.8) Then considering a null electromagnetic Ernst potential, Φ 0 = 0, the imaginary part of the parameters α, β, λ zero and choosing the specific values we exactly recover the (KN) transformation (B.7). In case of static metrics the latter further simplifies to (6.1), (6.2). Therefore (KN ) and (V ) are basically equivalent, up to gauge transformations, so they might be called collectively Harrison-Kramer-Neugebauer transformation.
As an explicit example we show the efficacy of the charging transformation (6.1), (6.2) on an asymptotically flat, static and discharged metric. For instance acting on the Schwarzschild metric, we produce the Reissner-Nordström black hole.
For simplicity we take the seed in spherical symmetric coordinates Note that this procedure, even though it always produces a non-equivalent (to the seed) solution, it might act not so predictably on an non-asymptotically flat solution, even if static.