Dyonic black holes at arbitrary locations

We construct and study stationary, asymptotically flat multicenter solutions describing regular black holes with non-Abelian hair (colored magnetic-monopole and dyon fields) in two models of $\mathcal{N}=2,d=4$ Super-Einstein-Yang-Mills theories: the quadratic model $\overline{\mathbb{CP}}^{3}$ and the cubic model ST$[2,6]$, which can be embedded in 10-dimensional Heterotic Supergravity. These solutions are based on the multicenter dyon recently discovered by one of us, which solves the SU$(2)$ Bogomol'nyi and dyon equations on $\mathbb{E}^{3}$. In contrast to the well-known Abelian multicenter solutions, the relative positions of the non-Abelian black hole centers are unconstrained. We study necessary conditions on the parameters of the solutions that ensure the regularity of the metric. In the case of the $\overline{\mathbb{CP}}^{3}$ model we show that it is enough to require the positivity of the"masses"of the individual black holes, the finiteness of each of their entropies and their superadditivity. In the case of the $ST[2,6]$ model we have not been able to show that analogous conditions are sufficient, but we give an explicit example of a regular solution describing thousands of non-Abelian dyonic black holes in equilibrium at arbitrary relative positions. We also construct non-Abelian solutions that interpolate smoothly between just two aDS$_{2}\times$S$^{2}$ vacua with different radii (dumbbell solutions).

10-dimensional Heterotic Supergravity. These solutions are based on the multicenter dyon recently discovered by one of us, which solves the SU(2) Bogomol'nyi and dyon equations on E 3 . In contrast to the well-known Abelian multicenter solutions, the relative positions of the non-Abelian black-hole centers are unconstrained.
We study necessary conditions on the parameters of the solutions that ensure the regularity of the metric. In the case of the CP 3 model we show that it is enough to require the positivity of the "masses" of the individual black holes, the finiteness of each of their entropies and their superadditivity. In the case of the ST[2, 6] model we have not been able to show that analogous conditions are sufficient, but we give an explicit example of a regular solution describing thousands of non-Abelian dyonic black holes in equilibrium at arbitrary relative positions. We also construct non-Abelian solutions that interpolate smoothly between just two aDS 2 ×S 2 vacua with different radii (dumbbell solutions).

Introduction
One of the most fascinating features of extremal black-hole solutions is that they can be superposed or combined, following certain rules, into solutions that describe several of these objects in equilibrium. Nowadays, these solutions are referred to as multicenter solutions, to encompass more general cases in which some of the objects associated to the "centers" are not black holes.
These solutions exhibit very interesting properties which we are going to review later on, but the most striking of them is that they exist at all. The existence of stationary solutions describing several gravitating objects in equilibrium is commonly (and correctly) attributed to cancellation between attractive gravitational forces and repulsive electric or magnetic forces. However, apparently, there are no self-interaction terms for the electromagnetic fields in the actions of the theories in which these solutions exist (e.g. in the Einstein-Maxwell theory, which admits the Majumdar-Papapetrou (MP) solutions [1,2] describing extremal Reissner-Nordström black holes in static equilibrium). It is, therefore, bewildering that the electromagnetic fields know that two centers with fields that correspond to charges of the same kind must repel each other.
It is useful to remember what the situation in absence of gravity is like. In that case, we are used to place point-like charges of arbitrary values at arbitrary points in space and then find the corresponding electrostatic field which solves all the Maxwell equations with those sources. This is possible because the Maxwell field is Abelian and it does not know what is the interaction between those centers nor whether they can be in static equilibrium or should be hold by other forces in the chosen positions unless interaction terms such as the worldline actions for charged particles, embodying the Lorentz force, are added to the theory.
In contrast, in a non-Abelian theory such as General Relativity, the interaction between two mass centers and their motion is completely determined by the field equations, as shown by Einstein, Grommer, Infeld, Hoffmann and others in Refs. [3,4,5,6,7,8]. General Relativity knows that it is not possible to have two Schwarzschild black holes in static equilibrium because self-interaction is built-in and regular static multicenter solutions simply do not exist.
Regularity is, evidently, a very important condition in this discussion because there are indeed solutions describing an arbitrary number of Schwarzschild black holes placed at arbitrary points in a straight line: the Israel-Khan solutions [9]. However, these solutions have conical singularities in the lines that join every two contiguous black-hole centers, the deficit angle being related to the Newtonian force acting between them. These singularities can be interpreted as struts exerting an additional force to compensate the gravitational attraction and hold the black holes in their positions. 1 Many of the singularities that occur in multicenter solutions can be interpreted along the same lines: they show that external forces are needed to hold the configuration in equilibrium. Therefore, we will be interested in the conditions required to make the singularities disappear and, ultimately, we will only consider regular solutions.
As we have stressed, the Maxwell equations in curved backgrounds do not contain any electromagnetic self-interaction terms. The reason why the MP solutions are possible must, therefore, lie entirely in the gravitational interaction and, more specifically, in the electromagnetic interaction energy which is implicitly contained in the electromagnetic energy-momentum tensor. Gravity may not know directly about electromagnetic interactions between charged particles but it does know about all the interaction energies. In the end, this is equivalent to knowing the interactions themselves well enough as to determine the equations of motion of the mass centers, as shown by Einstein et alia, and also of charge centers, as shown by Wallace and Infeld in the interesting but less well known Refs. [11,12,13]. 2 1 In an infinite periodic array of Schwarzschild black holes the total gravitation force over each of them vanishes and the conical singularities disappear [10]. 2 In Ref. [14], Brill and Lindquist studied the time-symmetric initial-date problem for several nonextremal Reissner-Nordström black holes and considered the contribution to the total energy in the common asymptotically-flat region of the gravitational and electrostatic interaction energies, but no connection between their values and the possibility of evolving the initial data into a completely regular static solution (a MP solution) was made. Similar solutions for the time-symmetric initial-data problem in Einstein-Maxwell-dilaton gravity and in models of N = 2, d = 4 supergravity are known [15,16] and This mechanism is, obviously, much more general and explains, for instance, the existence of static multi-D-brane solutions in superstring theory effective field theories (supergravities) in spite of the fact that, in the underlying fundamental theory, Dbranes have very complex interactions with a very delicate cancellation associated to supersymmetry [17].
The existence of static multicenter solutions, that we have almost grown used to take for granted, is, therefore, a small wonder and a signal that supergravities and other effective field theories including gravity encode a large amount of information of the original (superstring or other) theory.
The equilibrium of forces (or, actually, of interaction energies) required by these solutions can be achieved in more complicated situations, giving rise to stationary multicenter solutions. The first family of solutions of this kind was found by Perjés, Israel and Wilson (PIW) [18,19] in the Einstein-Maxwell theory. They are an extension of the MP family in which each center can have higher momenta of the gravitational and electromagnetic fields and the whole spacetime can also global momenta. If one wants the centers to be regular black holes 3 the dipole and higher momenta of each center must vanish, in agreement with the no-hair theorem. 4 Still, these solutions can present other pathologies such as Misner strings or closed timelike curves (CTCs) and Hartle and Hawking proved in Ref. [21] that the only regular solutions in the PIW family are those of the MP subfamily. 5 Generalizations of the PIW family were found in pure N = 4, d = 4 supergravity ("SWIP" [22,23,24]) and later in N = 2, d = 4 supergravity coupled to vector multiplets, [25,26,27] typically as timelike supersymmetric solutions. However, these families were only used to construct static (MP-like) multicenter solutions and the full potential of these solutions remained hidden until Denef and Bates showed in Refs. [28,29] how to construct completely regular multicenter solutions describing many static black holes with electric and magnetic charges (dyons) with global angular momentum in N = 2, d = 4 supergravity coupled to vector multiplets.
The source for the angular momentum of these solutions is the angular momentum of the electromagnetic fields due to the presence of Dirac monopoles and electric charges in different places. As a matter of fact, the global angular momentum is proportional to the symplectic-invariant Dirac-Schwinger-Zwanziger quantization condition for dyons and it should, therefore, be quantized. Generically, these solutions have Misner strings (the gravitational analog of Dirac strings) which can only be avoided at the price of introducing CTCs [30]. The cancellation of the sources of Misner strings imposes constraints on the charges and location of the centers 6 . These constraints are extremely hard to solve for 3 or more centers.
The fact that all the multicenter solutions mentioned so far carry Abelian dyonic could also be studied from the same point of view. 3 In the Einstein-Maxwell theory this seems to be the only way to avoid having naked singularities. 4 This is also the only way to have globally defined unbroken supersymmetry [20]. 5 We review this result in Appendix B from our own point of view. 6 We will review these constraints in Section 2. charges only suggests that a possible reason for the typical presence of Misner strings is, precisely, the Abelian nature of the fields and it also suggests that they could be avoided by the use of non-Abelian fields. However, no non-Abelian multicenter families of solutions have been constructed so far and it is the purpose of this paper to do it for the first time. 7 Here we are going to focus on the 4-dimensional case and in a forthcoming paper we will consider the 5-dimensional one [34], although some of the considerations made here will also apply to that case. Finding gravitating solutions with genuinely non-Abelian fields is a very complicated problem due to the non-linearities of the equations and, therefore, so far there have been no attempts to construct multicenter solutions beyond those mentioned in footnote 7. Actually, in the context of the Einstein-Yang-Mills (EYM) and Einstein-Yang-Mills-Higgs (EYMH) theories, even the single-center solutions are only known numerically [35,36] 8 . This makes them very difficult to study, interpret and generalize. As argued in Ref. [39], embedding a field theory in a supersymmetric one provides new tools to tackle the problem and construct new solutions, specially if one assumes that they preserve some supersymmetry, but it is necessary to use N > 1 supersymmetry 9 for the kind of solutions we are after. In general, the EYM and EYMH theories are not consistent truncations of any N > 1 supergravity, and, therefore, if we want to use the solution-generating techniques provided by supersymmetry, we must consider the simplest extended supergravities that include non-Abelian Yang-Mills fields, which we have called Super-Einstein-Yang-Mills (SEYM) theories. Typically, they come equipped with scalar fields that play the rôle of (usually adjoint) Higgs fields and, typically, low values of N give more freedom to choose the gauge group. Thus, N = 2 is the optimal value in 4 and 5 dimensions. 10 N = 2, d = 4 SEYM theories are theories of N = 2, d = 4 supergravity coupled to vector multiplets in which some subgroup of the isometry group of the Special Kähler scalar manifold has been gauged. These theories are the simplest which include YM fields and have a positive semidefinite scalar potential. This forces the timelike supersymmetric solutions to be asymptotically flat because the asymptotically-DS 4 ones cannot be supersymmetric. If the gauge group has an SU(2) factor, it is also possible to use it to gauge simultaneously the SU(2) factor of the U(2) R-symmetry group. The resulting theory has a potential that allows for asymptotically-aDS 4 solutions, but it is a much more complicated theory and only a few solutions (none of them describing 7 A 2-center solution that describes two SU(2) gravitating BPS magnetic monopoles in equilibrium (an 't Hooft-Polyakov and a Wu-Yang monopole) in N = 2, d = 4 Super-Einstein-Yang-Mills (SEYM) theories was constructed in Ref. [31], using the solutions of the Bogomol'nyi equations found in Refs. [32,33]. However, these are just particular solutions which are very hard to generalize to a higher number of centers. It is also possible to construct solutions with many Wu-Yang magnetic monopoles, but these are equivalent, up to a singular SU(2) gauge transformation to solutions with as many Dirac monopoles embedded in SU(2) and should not be considered as genuinely non-Abelian. 8 See also See Refs. [37,38]. 9 That is, more than 4 supercharges. 10 The 5-dimensional supergravity theories with 8 supercharges will be referred to as N = 1 theories because it is the minimal value in d = 5. black holes) are known [40], even though the most general timelike supersymmetric solutions have been characterized in Ref. [41]. Further generalizations are possible in presence of hypermultiplets, but here we are going to stick to the simplest possibility. 11 The solution-generating methods needed to construct non-Abelian solutions of N = 2, d = 4 SEYM theories were found in Ref. [45] and they have been successfully applied to construct, in fully analytical form, several interesting supersymmetric single-center solutions with genuine non-Abelian hair such as global monopoles and extremal static black holes [46,47,31,48] and the non-Abelian 2-center solutions mentioned in footnote 7.
As we will see, the supersymmetric solution-generating technique employed requires solving the non-Abelian Bogomol'nyi equations in E 3 [49]. The solutions to these equations are BPS magnetic monopoles such as the SU(2) 't Hooft-Polyakov monopole in the Bogomol'nyi-Prasad-Sommerfield limit [50,51,52], the SU(2) Wu-Yang monopole [53] or the colored monopoles found by Protogenov as part of its full classification of the possible spherically-symmetric solutions [54] and which can be extended to other gauged groups [48]. All of them have been used to construct regular black holes or gravitating (global) monopoles in N = 2, d = 4 SEYM theories. The multicenter solutions of the Bogomol'nyi equations are expected to describe several of these magnetic monopoles in equilibrium but only the very restricted or trivial examples discussed in footnote 7 were known until very recently.
In Ref. [55], based on the results of Ref. [56], one of us found a multicenter solution of the SU(2) Bogomol'nyi equations on E 3 describing an arbitrary number of colored monopoles in equilibrium. Furthermore, this configuration was generalized to describe colored dyons through the inclusion of electric non-Abelian sources. This multi-colored dyon solution will provide the basis to construct non-Abelian multicenter solutions in N = 2, d = 4 SEYM theories.
Colored magnetic monopoles are very interesting solutions that behave as Wu-Yang monopoles near the origin but have asymptotically vanishing magnetic monopole charge. 12 In Ref. [59] we showed that they are related via dimensional oxidation à la Kronheimer [60] to the BPST instanton [61] and the multi-colored monopole solution corresponds to a multi-instanton solution in a non-trivial hyper-Kähler space [34,56].
As in the Abelian case, multi-colored dyon solutions are stationary, rather than static. However, as we are going to see, these never gives rise to Misner strings and the positions of the dyons can be chosen completely at will. This is one of the main properties of the non-Abelian multicenter solutions that we are going to construct here. Another important property is that, due to the rapid fall-off of the non-Abelian fields at spatial infinity, the non-Abelian field do not give a net contribution to the global angular momentum.
The regularity of multicenter solutions is not guaranteed by the absence of Misner strings alone. It is necessary to study the complete metric and, in particular, the so-called "metric function" e −2U defined in Eq. (3.3), whose behavior determines the regularity of the black-hole horizon at each center and which must not vanish anywhere else. We are going to look for general conditions guaranteeing that this is the case and, at least for some models, we are going to see that they have very reasonable physical interpretations.
This paper is organized as follows: in Section 2 we set up the problem of finding non-Abelian, timelike supersymmetric, multicenter solutions of N = 2, d = 4, 5 SEYM theories, introducing the multi-colored dyon solution. In Section 3 we focus on the 4-dimensional case and apply the technique to two models of SU(2) N = 2, d = 4 SEYM (the CP 3 model in Section 3.1 and the ST [2,6] model, which can be embedded in Heterotic Supergravity, in Section 3.2), finding solutions whose regularity conditions we will study in full detail in terms of masses and entropies. Section 4 contains our conclusions. In Appendix A we briefly review N = 2, d = 4 SEYM theories. In Appendix B we revise Hartle and Hawking's result on the non-existence of stationary multi-black-hole solutions in the Einstein-Maxwell theory (which is just the bosonic sector of pure N = 2, d = 4 supergravity).

Setting up the problem
The problem of finding timelike supersymmetric solutions of N = 2, d = 4 SEYM theories and timelike or null supersymmetric solutions with an additional isometry of N = 1, d = 5 SEYM theories 13 boils down to the far simpler problem of finding functions Φ Λ , Φ Λ and vector fieldsȂ Λ r 14 in Euclidean 3-dimensional space E 3 solving the following three sets of equations: whereD r is the gauge covariant derivative in E 3 with respect to the connectionȂ Λ r . 13 These theories are briefly reviewed in Appendix A. The N = 1, d = 5 SEYM case will be dealt with in a forthcoming paper [34].
The first set of equations (2.1) are just the Bogomol'nyi equations [49] for a set of real, adjoint, Higgs fields Φ Λ and gauge vector fieldsȂ Λ r on E 3 . Due to their nonlinear structure (when the gauge group is non-Abelian) one has to solve simultaneously for Φ Λ andȂ Λ r . In the Abelian case, the integrability condition for these equations is the Laplace equation in E 3 , i.e. ∂ r ∂ r Φ Λ = 0; the Abelian vector fields are completely determined by the choice of harmonic functions Φ Λ and usually they are not written down explicitly.
For the SU(2) gauge group, which will be our main interest, all the sphericallysymmetric solutions were found by Protogenov in Ref. [54]. The BPS limit of the 't Hooft-Polyakov monopole [50,51,52], the SU(2) Wu-Yang monopole [53] and the so-called colored monopoles considered in Refs. [47,48] are, perhaps, the most interesting solutions. Only the 't Hooft-Polyakov monopole is regular, but, just as in the Abelian case, the singularity of the solution in E 3 needs not imply the existence of a spacetime singularity in the complete supergravity solutions. Actually, the singularities are typically associated to extremal black hole horizons.
Multicenter solutions of these equations, specially with the right properties necessary to construct multi-black-hole solutions, are extremely hard to find. In this paper we will use the multicenter solutions of the Bogomol'nyi equations found by one of us in Ref. [55] to construct multi-center black-hole solutions in 4 and 5 dimensions. This solution, which will be reviewed in the next section, is based on the multi-instanton solutions of Etesi and Hausel Ref. [56] and on the general relation between instantons in hyperKähler spaces and BPS monopoles on E 3 found by Kronheimer in Ref. [60].
The second set of equations (2.2) is a set of linear equations for the scalar fields Φ Λ . For SU (2) and, more generally, for compact groups, one can always use the trivial solution Φ Λ ∝ Φ Λ , which also satisfy Eqs. (2.3). However, a more interesting set of solutions has been found in Ref. [55] and we will make use of them. In the Abelian case, again, the Φ Λ are harmonic functions in E 3 : The third equation, (2.3) is the integrability condition of the equations that defines the 1-form ω r that appears in the 4-and 5-dimensional metrics. If we use the other two sets of equations, it seems to be automatically satisfied. However, since, typically, the fields Φ Λ , Φ Λ have singularities, the first two sets of equations may not be identically satisfied at the locus of the singularities. When this happens, the 1-form ω r still exists, but it can only be defined locally: it will exhibit Dirac-Misner string singularities [30] that can only be cured by defining different ω r which are regular in different patches and identifying these solutions in the overlaps up to "gauge transformations" that can be identified as coordinate transformations in the time direction. The consistency of these construction requires a periodic identification of the time coordinate with the consequent loss of asymptotic flatness. For this reason, Eq. (2.3) is required to hold everywhere and, at the loci of the singularities, this condition leads to nontrivial equations in the Abelian case which generically (for non-trivial Φ Λ ) constrain the relative distances of the pairs of black holes in terms of their charges and the moduli [28,29]. We will see that the solutions found in Ref. [55] do not imply any such constraints because they solve identically Eqs. (2.1-2.3) at the would-be singularities.
Given a solution Φ Λ , Φ Λ ,Ȃ Λ r of the above equations there are three sets of rules that allow us to construct timelike supersymmetric solution of N = 2, d = 4 SEYM theories and timelike or null solutions with an additional isometry of N = 1, d = 5 SEYM theories respectively. The functions and 1-forms Φ Λ , Φ Λ ,Ȃ Λ r will be the building blocks of the physical fields of the solutions. We will review the rules for the 4-dimensional case in Section 3 where we will construct and study explicit solutions of several supergravity models with a single non-Abelian SU(2) sector. Now we are going to set up the general problem of solving those equations and we are going to review the solutions found in Ref. [55] to which we will henceforth refer to as the multi-colored dyon.

The multi-colored dyon solution
The indices Λ, Σ, . . . that label the vector fields can be split into those corresponding to the Abelian and non-Abelian (SU(2)) sectors. Labeling the former with λ, σ, . . . and the latter with A, B, . . ., which will only take three values 15 (2.11) 15 We can always call these values 1, 2, 3 for convenience. Then, we can use the same labels for the Cartesian coordinates in E 3 , which simplifies considerably the notation. 16 Our conventions for the SU(2) objects are as follows: the structure constants are f AB C = +ε ABC = +ε AB C (the upper or lower position of the indices, which we will choose for essentially esthetic reasons, is irrelevant) and the covariant derivative and gauge field strength arȇ In some cases we use the following vector notation We will also use the notation The integrability conditions of Eqs. (2.7) are ∂ r ∂ r Φ λ = 0, which are solved by harmonic functions in E 3 , as mentioned above. The explicit form of the corresponding Abelian vector fieldsȂ λ r will not be required in what follows. It will be sufficient to know that they exist. Eqs. (2.9), which are also solved by harmonic functions in E 3 , can be interpreted as the integrability conditions of Abelian Bogomol'nyi equations for dual vector fieldsȂ λ r , but we will not need to know their explicit forms, either.
In order to obtain multi-center black-hole solutions, the harmonic functions Φ λ , Φ λ must be of the form for some points x α whose positions may be constrained by the integrability equations (2.11). 17 As shown in Ref. [55], Eqs. (2.8) are solved by 18 for real functions P satisfying Harmonic functions P of the form satisfy the above equation everywhere in E 3 , including at the locus of the singularities x = x α . For just one singularity ( x 1 = 0) and positive coefficients P 0 , P 1 , the corresponding solution of the Bogomol'nyi equations corresponds to a colored monopole [47,48]. The behavior of the gauge fields at infinity is such that using the standard definition of magnetic charge one gets zero. The non-Abelian fields, in fact, do not seem to contribute to any of the conserved charges defined at spatial infinity (mass or angular momentum, as we are going to see). The behavior of the gauge fields near the singularity r = 0, though, is the same as in the SU (2) Wu-Yang monopole case and they seem to contribute to the quantities that can be defined in the near-horizon limit, such as the Bekenstein-Hawking entropy, in exactly the same way as the Abelian fields corresponding to electric or magnetic charges. Then, one would naively conclude that the addition of a colored monopole to an Abelian black hole does not modify the asymptotic behavior (a clear violation of the no-hair and uniqueness "theorems") but it does modify the entropy, diminishing it both in 4 and in 5 dimensions [31,62,63]. However we have recently shown that, at least in the simpler 5-dimensional cases studied in [39,57], this is just an illusion caused by an inadequate identification of the charges of the solution in terms of fundamental objects in string theory; actually the non-Abelian sources modify the asymptotic charges but not the entropy. We expect this to be the appropriate interpretation in more complex configurations as well [58].
Let us now consider Eqs. (2.10). Apart from the trivial possibility Φ a = KΦ a , the following solutions were found in Ref. [55]: The simplest way to satisfy this equation is to choose Q as a harmonic function on E 3 with the same poles as P: With this choice, Eqs. (2.10) are satisfied everywhere in E 3 , including at the singularities of Q and P. Since Eqs. (2.8), whose integrability conditions arȇ are also satisfied everywhere for the chosen P, it is to be expected that Eq. (2.11) do not get any contribution from the non-Abelian sector. As a matter of fact, The Abelian sector of Eq. (2.11) contains terms proportional to δ (3) ( x − x α ) for all the poles x α and we need them to vanish identically for the reasons explained above. Requiring the coefficient of each delta function to vanish leads to Summing these equations over the index N we get a constraint relating the coefficients of the poles Φ λ α (which are proportional to each center's charges) to the constant terms Φ λ 0 which are related to the values of the scalars at infinity (moduli): This condition can be interpreted as requiring the vanishing of the global NUT charge of the spacetime [20] to avoid global Dirac-Misner strings or global periodic time. The conditions derived above for each center have the same meaning and, if the charges have been chosen, they constrain the relative positions of the centers. These constraints must be compatible with the triangle inequalities | any triplet of poles β, α, γ and this may not always be possible. Since our main interest lies in the non-Abelian sector, we will not discuss these equations in more detail, as they have already been thoroughly studied in the literature. It suffices to stress that the non-Abelian solution of Ref. [55] does not lead to any restrictions on the relative positions of the centers whatever the choices of coefficients P 0 , P α , Q 0 , Q α . Since Eqs. (2.11) are the integrability conditions of another set of equations, it is worth taking a look at the solutions of the latter associated to the choices made here. The equations we are talking about are those determining the components of the 1form ω r defined on E 3 : where ω (N)A stands for the (non)-Abelian contribution: If the integrability equations are satisfied, ω A can be defined in a single patch. The construction of the exact solutions is reviewed, for instance, in Ref. [44]. ω N A was found in Ref. [55] to be given by ; this is too fast to contribute to the asymptotic charges. Near the center x * In order to determine if ω N A contributes to the near-horizon limit we consider ω N A r dx r in spherical coordinates centered at x = x * to find that it ω N A r dx r ∼ rdϕ where r is the local radial coordinate. Then, if in this limit the tt component of the 4-dimensional metric e 2U ∼ r 2 , the contributions of ω N A will be subleading and the solutions will have the usual aDS 2 ×S 2 near-horizon limit.
This concludes the general discussion. We are now ready to construct 4-dimensional solutions from the building blocks we have introduced and studied here.

The metric takes the form
where ω = ω r dx r is the above spatial 1-form and the metric function e −2U is given by the Hesse potential e −2U = W(I ) .

5.
The scalar fields are given by 6. The components of the vector fields are given by

6)
3.1 Solutions of the CP 3 model

The model
The CP 3 model is characterized by the quadratic prepotential All we need to construct supersymmetric solutions is the CP 3 Hesse potential More details on these models can be found in Refs. [31,44].

The solutions
The Abelian sector of the model is determined by the complex harmonic function H ≡ Φ 0 + 2iΦ 0 and the non-Abelian one by the two triplets of real functions Φ A and J A ≡ −2Φ A . According to the general discussion, if we use the multi-colored dyonic solution we only need to solve the Abelian part of the integrability equations (2.11). For just one Abelian vector the only possibility is e H ∝ m H or, equivalently, H = e iγ H for some real harmonic function H and a constant phase γ. Then, according to the discussion in Section 2, the solution is given in terms of three harmonic functions H, P, Q with singularities at the same N isolated points x = x α by The metric function, the 1-form ω = (ω r ), the scalar fields and the scalar potential can be written as The vector fields of the solution can be constructed using the general recipe, Eqs. (3.6) and (3.7), but we will not do it explicitly here as we are more concerned with the regularity of the metric and scalar fields.

Spherically-symmetric and dumbbell solutions
As a warm-up exercise, it is convenient to start by the construction of a single-center solution of this model, which is static because with a single center necessarily must have Φ ∝ J . This was already done in Ref. [47] (with less independent parameters), but here we will show that there is also a Robinson-Bertotti dumbbell solution similar to the one recently discovered in a 6-dimensional context in Ref. [68]. These singlecenter dumbbells are obtained by setting to zero the constant term in the harmonic functions of the Abelian sector. Without the non-Abelian colored monopole, we would simply obtain the standard Robinson-Bertotti aDS 2 ×S 2 solution, which is sometimes called a double extreme black hole. When the colored dyon is included the geometry gets modified. However, the non-Abelian field decays very fast with the distance and the original aDS 2 ×S 2 asymptotic is recovered. On the other hand, near the origin, the colored dyon contributes as just another "Abelian" charge and one also gets an aDS 2 ×S 2 spacetime, albeit with different radius (smaller than the original). Thus, the Robinson-Bertotti dumbbell solution interpolates between two aDS 2 ×S 2 spacetimes of different radii. 20 Taking N = 1 (and suppressing the indices that label the centers), we get Let us analyze the asymptotically-flat (h 2 = 1) case first. It is convenient to define 18) in terms of which the metric function takes the form where we have defined the manifestly positive function which varies smoothly from 0 at r = 0 to 1 at r = ∞.
In the above form the metric function is, therefore, manifestly positive if M (which is the mass) and E (which will be seen to be the entropy times π) are both positive. In the asymptotic and near-horizon limits we find respectively showing that the colored dyon field cannot be seen asymptotically but does contribute to the near-horizon geometry: i.e. it appears in the entropy E and in the covariant attractor value of the scalars [46]. Setting h = 0 (with λ = 0) we get the dumbbell solution The metric function interpolates smoothly between E/r 2 at r ∼ 0 and p 2 /r 2 at r ∼ ∞ while staying always positive. The scalars interpolate between two covariantly-constant attractors which have different r-dependence because the gauge connection behaves differently in both limits.

Multicenter solutions
For more than one center the metric function is given by Inspired by the single-center case, we now define as then the term in line (3.25) can be expressed as whereas the term in line (3.26) can be written as The last terms in Eqs. (3.30) and (3.31) are easily seen to be positive. First, we define the positive functions K α and then, we write 1 − s α s β r α r β P 2 n α · n β = (3.34) from which the positivity is paramount because the functions K γ are positive and n α · n β ∈ [−1, 1]. Since there is a term (1 + η α η β ) multiplying the whole second term we need to impose the condition that The function R α is a generalization of the function R defined in Eq. (3.20) for the singlecenter case and varies from 0 at r α = 0 to 1 at infinity or at any other point r β = 0 β = α.
The functions R αβ are also bound by 0 and 1 and are equal to 1 at all the points r γ = 0 and at infinity. The metric function takes the final form and its positivity can be guaranteed by imposing the conditions for all α, β and the sign condition (3.35). The only poles in the metrical factor (the zeroes of g tt = e 2U , and, hence, the horizons) are the ones at the points r γ = 0. As can be seen in the asymptotic expansion r → ∞, the physical meaning of the first set of conditions is that the mass that each individual black hole would have if it were isolated must be positive. The meaning of the other two sets of conditions comes from the study of the near-horizon limits r α → 0. In that limit the dominant term is the coefficient of 1/r 2 α the value of α we are dealing with. Since R α vanishes precisely at r α = 0 only the constant part of the coefficient, E α , survives and we get an aDS 2 × S 2 geometry with metric so E α , as the notation suggests, is the entropy of the α'th black hole up to a factor of π. Thus, we are asking for all the individual extremal black holes to have a regular horizon.
The third set of conditions amounts, then, to the requirement that the entropy of a black hole whose charges are those of the pair αβ combined should be larger than the sum of the individual entropies, i.e. we are assuming the superadditivity of the entropy.
In some special cases, though, the third set of conditions is more restrictive than necessary to ensure the regularity of the metric. Notice that the metric function can be positive everywhere even if the second line in Eq. (3.36) has negative constant coefficient. For instance, in the two centers case, the constant coefficients of the 1/r 2 α and 1/(r α r β ) terms are and can be rewritten in this form: This combination is non-negative everywhere if which is a weaker condition for which we have, however, no clear physical interpretation.
The physical scalars are regular everywhere and can be written as and vanish as O(r −2 ) at infinity while. At the αth center they take the covariantlyconstant attractor value e −iγ (1+iη α ) gp α n α . In the previous discussion we have ignored the presence of a non-trivial 1-form ω r dx r in the metric given by Eqs. (2.26) or (3.13) because, asymptotically, it vanishes faster than any other function in the metric and, in the near-horizon limits, they are also subleading. However, we must see if its presence gives rise to pathologies such as closed timelike curves. For the harmonic functions P and Q in Eq. (3.10) it takes the explicit form Far away from the centers, and for λ = 0, and choosing coordinates such that v is parallel to the z axis and g ϕϕ = e 2U ω 2 ϕ − e −2U r 2 sin 2 θ is clearly negative in that limit. When λ = 0, then ω ϕ ∼ r −2 asymptotically, decaying still too fast to contribute to the angular momentum or to modify the sign of g ϕϕ .
In the near-horizon limit x → x * , where x * denotes the coordinates of the center we are zooming on We can choose adapted coordinates such that now u is parallel to the z axis, so we can write to leading order ω ∼ u(ydx − xdy) r = u * r sin 2 θdϕ , (3.47) and In that expression the second term is always smaller than the first term in this limit. Beyond this analysis, we have explored numerically the value of g ϕϕ for several, simple, multicenter configurations and have found that it can vanish (for instance,in a 2-center example, all along the the axis that contains both centers) but it never changes sign. See Fig. 1 for a simple two-center example.  (1, 2). The blue surface always lays below the yellow surface, whence g ϕϕ remains finite and positive.

The model
The ST [2,6] model is the cubic model with prepotential (3.49) where i = 1, 2 · · · , 6 labels the vector multiplets and where the fully symmetric tensor d ijk has as only non-vanishing components d 1αβ = η αβ , where (η αβ ) = diag(+ − · · · −) , and α, β = 2, · · · , 6 . (3.50) The 6 complex scalars parametrize the coset space (5) , (3.51) and the group SO(3) acts in the adjoint on the coordinates α = 4, 5, 6 that we are going to denote with A, B, . . . indices. These are the directions to be gauged. In order to construct solutions we only need the Hesse potential of this theory, which is given by We could have gauged any three of the directions 3, 4, 5, 6, and, therefore, the ungauged one could have been truncated from our model. However, as shown in Ref. [68], it is necessary to have one additional Abelian vector field to be able to uplift the solution to 6 dimensions and then to Heterotic supergravity [57]. Furthermore, with the extra Abelian vector multiplet, the model can be seen as the STU model coupled to an SU(2) triplet. This can be made manifest by combining the Abelian directions 2 and 3 as follows whose positivity leads to a constraint on the possible values of the imaginary parts of the scalar fields, a constraint that we will use later. More details on this theory and, in particular, on its relation with the toroidal compactification of the Heterotic string can be found in Refs. [31,62,39,57].

The solutions
For the sake of simplicity, we are going to consider solutions with non-vanishing functions I 0 , I A , I 1 , I + , I − , I A only. 21 It is convenient to redefine these functions, bringing the metric function e −2U , the 1-form ω and the scalar fields to the form while the vector fields are given by

63)
The explicit magnetic part of the SU(2) vector field, Ȃ , is determined by Φ, J , which we will choose as in the CP 3 model Eqs. (3.10) and (3.11). We rewrite them here for convenience: The Abelian functions Φ 0 , J 1 , J + , J − will be given by The above form of the metric function (3.57) has the interesting feature that the 1-form ω appears in it (the last term). If we switch off all the functions but Φ and J , e −2U = |ω| and we get a metric which is completely determined by ω, but which is not asymptotically flat neither free of singularities since |ω| can vanish.
We are going to work with the following alternative form of the metric function If we plug into Eq. (3.55) the values of the scalars, we find the condition (3.70) and, using this condition in the above form of the metric function we find a second regularity condition These conditions are necessary but not sufficient to ensure the regularity of the solution, which also requires (3.72)

Spherically-symmetric and dumbbell solutions
Again, we start by studying solutions with a single center that we conveniently place at x = 0, suppressing all indices α, β, . . . Since Φ ∝ J the 1-form ω vanishes and the solutions are necessarily static.
Imposing the standard normalization of the metric at spatial infinity and studying the asymptotic behavior of the scalar fields we identify the integration constants and we will take q 1 , p 0 > 0 and q + q − > 0 with sign(q ± ) = sign(h ± ) = −sign( mZ ∓ ). Let us consider the first regularity condition Eq. (3.70). Expanding the functions in the left-hand side we find where, given the values of the h constants, is the non-negative function given in Eq. (3.20) so the combination Σ must be positive Doing the same with the second regularity condition Eq. (3.71) we get where which will be manifestly positive if the combination Finally, let us consider the third regularity condition Eq. (3.72). All the terms that originate in the product of the first two terms are manifestly positive if Ω and Σ are positive. The only negative terms come from the last term and are of O(r −4 ) We just need to compare them with the positive O(r −4 ) terms coming from the first two terms, i.e. we have to consider Thus, the third regularity condition is fulfilled if we require that Observe that, if this condition is satisfied, the entropy is given by The conditions that we have imposed on the charges and the asymptotic values of the scalars automatically ensure the positivity of the mass, which is given by Setting all the h constants to zero, we get a dumbbell solution with metric function and the square root function interpolates smoothly between E at r = 0 and p 0 q 1 q + q − at r → ∞, which is the value one would get in the purely Abelian solution. The scalars also interpolate between a covariant attractor and an Abelian attractor.
Let us now move to the multicenter case.

Multicenter solutions
The presence of more centers does not change the asymptotic values of the scalars and, therefore, the values of the constants h are unchanged and given by Eqs. (3.73). We impose on the charges of each center the same conditions as in the single-center case, that is: Moreover, the four harmonic functions cannot change sign anywhere, as if any of them becomes zero then metric function is imaginary, among other pathologies. Then we can include the conditions The first regularity condition Eq. (3.70) can be rewritten in the form where R α and R αβ are the functions defined in Eqs. (3.33) and (3.34), respectively, and and its positivity is manifest by requiring Only the first of these conditions (Σ α > 0 , ∀α) is independent, though. It implies that q + α > η 2 α g 2 q − α and, substituting in In a similar way, we rewrite the second condition Eq. (3.70) in the form where now The positivity bound is obviously satisfied by requiring and one can show, as before, that the first condition implies the second. Finally, let us study the third condition Eq. (3.72). Again, all the terms that come from the first two factors (corresponding to the first two conditions) are positive if the conditions that we have derived above are met. The negative contributions come from 22 and they have to be compared with other (positive) terms of the same order, O(r −4 ) and with the same structure. Let us first consider terms of the form r −4 α , which are dominant in the α th near-horizon region, The positivity of these terms is guaranteed by the positivity of Σ α and Ω α , which we have required before, and the reality of the entropy of each black hole: This implies that the metric function is well-defined in the neighbourhood of each black hole, provided the corresponding entropy is real. On the other hand, asymptotic flatness and the sign conditions on the parameters described above, which in turn imply positivity of the "masses", guarantee that the metric is also regular far away from any center. However, contrary to our experience with the CP 3 model, we have not been able to find a general analytical proof of the regularity of the metric due to the complexity of the ST[2, 6] model, as we will shortly see. Nevertheless, we expect most multicenter solutions asymptotically flat and with well-defined individual entropies to be regular everywhere, as the non-Abelian terms generally decay faster with distance than the Abelian harmonic functions. Because of their simplicity, let us consider the terms of the form r −2 α r −2 β : The coefficient of r −2 α r −2 β has constant terms and other terms which are linear and quadratic in R α and R αβ . The linear ones are manifestly positive. The quadratic terms add up to which is clearly negative when n α · n β = 0. However, it is bounded from above and below and its negative contribution can still be cancelled by the other terms. The constant terms are where we have defined ∆ αβ is positive under the assumptions we have made, because, for instance The second term is also positive, but the third is negative. Based on our previous experience with the CP 3 model, we can try to relate this coefficient to the superadditivity of the entropy, rewriting it as follows: This expression is not very enlightening as there is no simple way to show that the would-be positive terms in the first line are actually larger than the negative ones in the second. 23 Summarizing, not all the terms that appear below the square root sign in the metric function are positive definite and we have not been able to determine a set of conditions ensuring the positivity of the whole expression and the regularity of the metric 23 The superadditivity condition E αβ ≥ E α + E β or E 2 αβ − (E α + E β ) 2 ≥ 0 does not seem to lead to any identity that can be used directly in the terms at hands. From the conditions Σ αβ ≥ Σ α + Σ β and Ω αβ ≥ Ω α + Ω β we find Center 1 2 1 1 3 1 1 2 1 Center 2 1 2 1 1 1 1/2 3/4 1/2 function, which still might possible, in accordance with our experience with the CP 3 model.
To conclude this subsection we are going to give an explicit example of a completely regular two-center dyonic solution of this model. Our choice of charges for the two centers is given in Table 1. In this case, the coefficient of the r −2 1 r −2 2 term is the only one which is not manifestly positive and is given by However, since R 2 12 ≤ 1, this term is positive everywhere. Observe that the 1-form ω has exactly the same form as in the CP 3 model case and, as the analysis made in that case showed, it will have no effect on the regularity of the metric.
Adding and substracting 2E α E β we arrive to which cannot be used for our purposes.

Thousands of dyonic black holes
While we have not been able to prove the reality of the metric function for completely general configurations, we have argued that most solutions (if not all) described by our construction are well-behaved, provided the "masses" and entropies of the individual black holes are positive and real. To add further support to this thesis, we now describe a very general solution composed of 6060 black holes whose regularity we have checked by numerical analysis.
The system is composed of two well differentiated clusters. The first cluster describes a set of 1480 pairs of black holes with the same charges as the two-center system presented at the end of previous section. The second cluster contains 3100 black holes whose charges have been chosen with a random generator, provided the conditions (3.86), (3.87) and (3.103) are met. Since the position of each black hole is free, those have been placed as depicted in Figure 2 for aesthetic reasons 24 .

Conclusions
In this paper we have constructed and studied the very first multicenter black-hole solutions with non-trivial non-Abelian fields corresponding to colored monopoles and dyons. These solutions describe regular black holes in equilibrium when certain conditions (which we discuss below) are met. In general, they are stationary, although they have vanishing angular momentum unless the Abelian fields contribute to it. If these Abelian contributions are absent, the black holes can be have arbitrary positions. The main ingredients in the construction of these solutions are 1. Unbroken supersymmetry, which provides us with a very powerful solutiongenerating technique [46]. The use of this technique is only possible if one considers (as we have done here) the simplest N > 1 supersymmetric generalizations of the Einstein-Yang-Mills system. As a reward for considering this generalization the solutions are obtained in a completely analytical form. This, in its turn, allows for a deeper understating of the solutions.
2. The multi-colored dyon solution of Ref. [55], which is the main building block of the physical fields of the 4-dimensional spacetime solution. This solution solves the integrability equations (2.3) everywhere independently of the positions of the centers.
Profiting from the analytical form of the metrics obtained, we have tried to determine general conditions on the charges and moduli guaranteeing regularity. In the CP 3 model with any number of dyonic centers at arbitrary positions, we have shown that the positivity of each of the "masses" and entropies and the superadditivity condition for every pair of black holes are sufficient to guarantee regularity. Actually, as we have seen for just two centers, a condition weaker than superadditivity can also be sufficient. In the ST [2,6] model with only some Abelian vectors active, we have not been able to prove that similar conditions for an arbitrary number of dyonic centers are sufficient, although we have explicitly constructed and checked numerically highly non-trivial regular solutions with thousands of black holes. Also, we have shown that very simple conditions are sufficient when there are only magnetic monopoles at the centers.
We have also found that, removing the constant part of the harmonic functions in the spherically symmetric (single-center) solutions one can obtain solutions that interpolate between two aDS 2 ×S 2 vacua with different radii that we have called dumbbell solutions. They are the 4-dimensional version of similar 6-dimensional solutions found in Ref. [68] interpolating between two aDS 3 ×S 3 vacua with different radii, also in a non-Abelian context. The existence of these solutions in the non-Abelian case 25 sug-gests the possible existence of an Euclidean instanton describing the decay of one vacuum into the other one. The aDS/CFT interpretation of the corresponding transition (if found) should provide interesting insights into this correspondence.
As we have discussed in Section 2 solutions to the same three sets of equations (2.1)-(2.3) can be used to construct timelike supersymmetric solutions of N = 1, d = 5 SEYM theories using different rules to relate the building blocks that occur in those equations and the physical 5-dimensional fields. Typically, building blocks that lead to regular 4-dimensional solutions produce singular 5-dimensional solutions and vice versa. This means that the construction of 5-dimensional solutions will have to be studied independently. Work in this direction is well under way [34].
We have deliberately set aside for future work (already in progress [58]) the paradoxes created by the strange properties of the colored dyons which do not seem to contribute to the mass or any other asymptotic charge (so they behave as non-Abelian hair) but, nevertheless, do seem to contribute to the entropy. In the 5-dimensional case an analogous paradox was completely solved in Ref. [57] by the correct, string theoryinspired, reinterpretation of the Abelian charges and the identification of a globally regular solution supported by the non-Abelian field (a BPST instanton) [39]. Although we have not yet found globally regular solutions associated to the 4-dimensional colored dyons, we expect a similar resolution for this paradox, at least in the case of the ST [2,6] model, because the string theory embedding of the CP 3 model is unknown (or inexistent). 26 As mentioned in the introduction, the non-Abelian asymptotically-aDS case is much harder to deal with in SEYM theories. We are currently working on the generalization of the methods and solutions used here and we expect to report on our results soon [70].
To conclude, SEYM theories provide new tools to study the interplay between non-Abelian Yang-Mills and gravitational fields through the construction of a wealth of new, fully analytical solutions, some of which can be reinterpreted in the framework of string theory. As we have discussed, there are many directions to be explored and it is our purpose to follow some of them in the near future.
John Templeton Foundation grant 48222. TO wishes to thank M.M. Fernández for her permanent support.
A N = 2, d = 4 SEYM theories N = 2, d = 4 Super-Einstein-Yang-Mills (SEYM) theories can be seen as the simplest N = 2 supersymmetrization of the Einstein-Yang-Mills (EYM) theories. They are nothing but theories of N = 2, d = 4 supergravity coupled to n vector multiplets in which a (necessarily non-Abelian) subgroup of the isometry group of the (Special Kähler) scalar manifold has been gauged using some of the vector fields of the theory as gauge fields. 27 The necessary and sufficient conditions for the gauging of a non-Abelian subgroup of the global symmetry group to be possible are: 1. It must act on the vector fields in the adjoint representation.

2.
It must be a symmetry of the prepotential; see e.g. Ref. [46] for more details.
We will only be concerned with the bosonic sector of the theory, which consists on the metric g µν , the vector fields A Λ µ (Λ = 0, 1, · · · , n) and the complex scalars Z i (i = 1, · · · , n). The action of the bosonic sector reads In this expression, G ij * is the Kähler metric, D µ Z i is the gauge-covariant derivative N ΛΣ is the period matrix and, finally, V(Z, Z * ) is the scalar potential Since the imaginary part of the period matrix is negative definite, the scalar potential is positive semidefinite, which leads to asymptotically-flat or -De Sitter solutions.
In the above equations, k Λ i (Z) are the holomorphic Killing vectors of the isometries that have been gauged 28 and P Λ (Z, Z * ) the corresponding momentum maps, which are related to the Killing vectors and to the Kähler potential K by for some holomorphic functions λ Λ (Z). Furthermore, the holomorphic Killing vectors and the generators T Λ of the gauge group satisfy the Lie algebras For the gauge group SU (2), which is the only one we are going to consider here, we use lowercase indices 29 x, y, z = 1, 2, 3 and the structure constants are f xy z = ε xyz , so [k x , k y ] = −ε xyz k z , [T x , T y ] = +ε xyz T z . (A.8) The equations of motion of the theory can be written in the following form: +4M MN F M µ ρ F N νρ + 1 2 g µν V(Z, Z * ) = 0, (A.9) where G Λ µν is the dual vector field strength F M µν is the symplectic vector of vector field strengths F M ≡ F Λ G Λ , (A.13) 28 The employed notation associates a Killing vector to each value of the index Λ in order to avoid the introduction of yet another class of indices and the embedding tensor (See e.g. the reviews [71]); it is understood that not all the k Λ need to be non-vanishing. 29 These will be a certain subset of those represented by Λ, Σ, . . .. M MN is the symmetric 2(n + 1) × 2(n + 1) matrix defined by and

B Supersymmetric multi-BH's in pure EM theory
Einstein-Maxwell gravity is equivalent to minimal N = 2, d = 4 supergravity (in fact it could be called the CP 0 model). The timelike supersymmetric solutions of this supergravity theory are nothing but the Perjés-Israel-Wilson family of solutions [18,19] which we can, then, study using the language and methods we use in other models of N = 2, d = 4 supergravity in the main text, recovering Hartle and Hawking's result [21] that the only regular solutions in this family of solutions are those of the Majumdar-Papapetrou subfamily [1,2]. Our starting point will be that one can only use in the construction of regular solutions harmonic functions with point-like singularities corresponding to electric or magnetic monopoles, but no higher multiplets of the electromagnetic field [20]. The metric function of pure supergravity is given by By assumption, the complex function H has the form and the metric function, conveniently normalized at infinity (h = e iγ ) can be written in the form where M α ≡ e(e iγ Γ * α ) , (B.5) is the mass of the α th black hole, is (up to a factor) the entropy of the α th black hole, and is (up to a factor) the entropy of a black hole with the charges of the α th and β th black holes combined. It is evident that the metric function will be regular if the masses are non-negative M α ≥ 0, the entropies corresponding to centers with non-vanishing mass are strictly positive E α > 0 and the entropy of the combination of two black holes is not smaller than the sum of the entropies of the individual black holes E αβ ≥ E α + E β . Given the expressions for the masses and entropies, it is also evident that the condition E α > 0 for M α > 0 is, actually, redundant.
We also have to examine Eq. Furthermore, the condition E αβ ≥ E α + E β is equivalent to This condition is also trivially valid for α = β, which corresponds to the condition E α ≥ 0. For two black holes N 2 = −N 1 and the α = 1 equation takes the form which, if the masses are positive, as required by the regularity of e −2U , is only be solved by N 1 = N 2 = 0 so the phases of Γ 1 and Γ 2 are both equal to e iγ . Then, and the 1-form ω vanishes identically.
For three black holes, if one of the N α vanishes, we recover the equations of the two-black-hole case, and the same conclusion. Let us, then, consider the case in which the three N α are different from zero. Eqs. (B.11) imply that the three masses are also different form zero. Due to Eq. (B.12), two of the N α will have the same sign and the third will have the opposite sign. With no loss of generality we can consider N 1 > 0 and N 2,3 < 0 (the other case differs only in a global sign). This means that whose l.h.s. and r.h.s. are, respectively, positive and negative definite by assumption.