Electromagnetic Quasitopological Gravities

We identify a set of higher-derivative extensions of Einstein-Maxwell theory that allow for spherically symmetric charged solutions characterized by a single metric function $f(r)=-g_{tt}=1/g_{rr}$. These theories are a non-minimally coupled version of the recently constructed Generalized Quasitopological gravities and they satisfy a number of properties that we establish. We study magnetically-charged black hole solutions in these new theories and we find that for some of them the equations of motion can be fully integrated, enabling us to obtain analytic solutions. In those cases we show that, quite generally, the singularity at the core of the black hole is removed by the higher-derivative corrections and that the solution describes a globally regular geometry. In other cases, the equations are reduced to a second order equation for $f(r)$. Nevertheless, for all the theories it is possible to study the thermodynamic properties of charged black holes analytically. We show that the first law of thermodynamics holds exactly and that the Euclidean and Noether-charge methods provide equivalent results. We then study extremal black holes, focusing on the corrections to the extremal charge-to-mass ratio at a non-perturbative level. We observe that in some theories there are no extremal black holes below certain mass. We also show the existence of theories for which extremal black holes do not represent the minimal mass state for a given charge. The implications of these findings for the evaporation process of black holes are discussed.

for static and spherically symmetric metrics, are of second order. It was additionally checked that the linearized equations on constant curvature backgrounds are proportional to the linearized Einstein tensor [10]. Higher-order versions of Quasitopological gravity have been later constructed [39][40][41]. However, again these theories are constrained by the dimension of spacetime as they only exist in D ≥ 5.
More recently, a more general class of theories containing Quasitopological and Lovelock gravities has been identified. The theories of this class are known as Generalized Quasitopological gravities (GQGs) [42,43], and their main characteristic is that they allow for static and spherically symmetric (SSS) solutions of the form this is, with g tt g rr = −1. From their defining property, it can also be proven that the equation of motion of f (r) is partially integrable and that the linearized equations on constant curvature backgrounds are of second order. In addition, it is verified in all the known examples that black hole thermodynamics can be studied analytically in these theories. There also exist more refined subsets of GQGs that possess single-function Taub-NUT solutions [44] and/or second-order cosmological Friedmann equations [45][46][47]. Unlike Quasitopological gravities, their generalized counterparts do exist in D = 4, and in fact, they have been shown to exist in all dimensions and at all orders [41]. In addition, they are more general than one would expect, and it turns out that they provide a basis for at least all operators in the sets {Riem n |n ∈ N 0 } and {Riem n ∇Riem∇Riem |n ∈ N 0 } in the effective action of gravity once field redefinitions are taken into account [48]. The identification of GQGs was triggered by the earlier construction of Einsteinian cubic gravity [49] as a theory with second-order linearized equations on constant curvature backgrounds in arbitrary dimension. The special form of black hole solutions in this theory was soon noticed [50,51] and this motivated the definition of Generalized Quasitopological gravities established in Refs. [42,43]. Since then, many applications of these theories have been discussed in the literature in a wide range of topics including cosmology, black holes, holography and phenomenology .
The definition of GQGs does not involve matter, so as an important extension one may consider coupling other fields to these theories. This is a nice exercise if the matter fields also respect the property of producing SSS solutions satisfying g tt g rr = −1 [74]. A simple and relevant example of this is provided by a minimally coupled Maxwell field, a case which has been explored in Refs. [51,55,64,66] in the context of Einsteinian cubic gravity and higher-order GQGs. It is also a satisfactory strategy to couple these theories to non-linear electrodynamics, as recently shown in [75]. However, these examples only involve minimally coupled gauge fields, which is a very restricted way of coupling a vector to gravity. In general, higher-derivative effective actions may contain all types of couplings between the fields present in them. In addition, non-minimal couplings could have interesting effects that do not show up in the minimally coupled case, so they may be worth exploring. On the other hand, the study of charged black holes with higher-derivative corrections is a current matter of research in the context of string phenomenology due to its relation with the Weak Gravity Conjecture (WGC) [76] -see e.g. [77][78][79][80][81][82][83][84][85][86]. In this respect, it may be interesting to analyze the effect of vector-curvature couplings on charged black holes, since those couplings indeed appear in stringy effective actions [87][88][89].
Thus, one interesting question is whether it is possible to find a family of theories analogous to Generalized Quasitopological gravities in the case of non-minimal couplings between the curvature and a gauge field. The goal of this paper is to show that such generalization is indeed possible and to study charged black hole solutions in the new theories. In particular, we show that a satisfactory extension can be achieved by searching for theories whose magnetically-charged static spherically symmetric solutions satisfy the condition g tt g rr = −1. Note that, even though we focus on theories with "simple" magnetic solutions, one can generate theories with analogous electric solutions by dualizing the vector field. We obtain an infinite number of Lagrangians belonging to this new class of theories, that we denote Electromagnetic Generalized Quasitopological gravities (EGQG) 1 . As we show, for all of these theories it is possible to study the thermodynamic properties of black holes analytically, and for some of them we can even write exact black hole solutions. The paper is organized as follows.
• In Section 2 we review some basic aspects of general L(R µνρσ , F αβ ) theories, namely, equations of motion, duality transformations, conserved charges and first law of black hole mechanics.
• In Section 3 we analyze the equations of motion of particular L(R µνρσ , F αβ ) theories for charged static spherically symmetric configurations and we establish the definition as well as the properties of Electromagnetic Generalized Quasitopological gravities.
• In Section 4 we study the "quasitopological" subset of EGQGs, corresponding to those theories for which the equation of motion of the metric function f (r) ≡ −g tt is fully integrable. We find two infinite families of Lagrangians belonging to this class and we exactly solve the equations for magnetically charged spherically symmetric solutions. We show that in many cases the solutions are non-singular, corresponding to regular black holes or smooth horizonless geometries. We then study the thermodynamic properties of black holes in these theories, showing that the first law of black hole mechanics holds exactly. In addition, we analyze the properties of extremal black holes.
• In Section 5 we study an infinite family of proper Electromagnetic Generalized Quasitopological gravities, i.e., those for which the equation for f (r) is not algebraic. We analyze how this equation can be solved to search for black hole solutions and we man-age to determine analytically the thermodynamic properties of these black holes. The extremal limit is also discussed.
• We discuss our findings and point out new directions in Section 6.
In addition, we include two appendices with some technical results.
2 Some generalities on L(R µνρσ , F αβ ) theories As a preliminary step before focusing on the main topic of this paper, in this section we review some aspects of general L(R µνρσ , F αβ ) theories that will be useful for later purposes.

Equations of motion
Let us consider a general gauge-and diffeomorphism-invariant theory for the metric tensor g µν and a U(1) gauge field A µ . The Lagrangian of such theory must be constructed from contractions of the Riemann curvature tensor R µνρσ and the field strength 3 F = dA using the (inverse) metric g µν , and we denote it by L(R µνρσ , F αβ ). 4 Then, we define the action by where we consider an arbitrary spacetime dimension n. Since the results obtained in this subsection are valid for any dimension, we will keep n to be arbitrary for now, but we would like to emphasize that the rest of the manuscript will refer to 4-dimensional theories, as it will become apparent. Also, from now on we set G = 1. If we consider pure theories of gravity (that is, with no coupling to electromagnetism), a general formula is known for the associated equations of motion -see e.g. [35]. Our first task in this document shall be to derive an analogous formula in the case of L(R µνρσ , F αβ ) 2 For some people QG = Quantum Gravity. 3 The most natural way to achieve a gauge-invariant theory is to impose that the whole dependence of the Lagrangian on the gauge field A takes place through its curvature 2-form F . 4 For simplicity we assume no covariant derivatives acting on Riemann curvature tensors nor on field strengths.
theories. In order to obtain the generalized Einstein's equations, we consider the variation of the action (2.1) and evaluate it on a vector of the form (δg µν , 0): where the term δR αβργ is evidently not independent from δg µν . Defining, for convenience of notation, it is possible to find that, up to total derivatives, On the other hand, let us work out in two different manners the Lie derivative L ξ L with respect to an arbitrary vector field ξ ∈ X(M ): If we take into account that Consequently, equation (2.2) may be re-expressed as Hence the gravitational equations of motion of any theory given by the action (2.1), representing the most general theory of gravity coupled to electromagnetism which includes all 5 We recall that ∂L ∂g µν = −gµα ∂L ∂g αβ g νβ .
possible terms constructed out of curvature tensors, metrics and field strengths, are given by the formula 6 On the other hand, the derivation of the generalized Maxwell's equation is straightforward and it yields

Duality transformations
A duality transformation of a U(1) gauge field in four dimensions has the effect of exchanging the initial vector field by a dual vector field whose equations of motion correspond to the Bianchi identity of the former, and viceversa. Intuitively, this operation corresponds to the exchange of electric and magnetic fields. An interesting property of (Einstein)-Maxwell's theory -and also of many extensions of it, such as N = 2, d = 4 SUGRA -is the fact that it is invariant under duality transformations. However, in general the dual theory does not need to coincide with the original one, and thus the duality transformation establishes a map between two different theories. 7 This is usually the case when one considers non-minimal couplings between the curvature and the field strength.
Let us show how the transformation works on a general L(R µνρσ , F αβ ) theory. We consider the action (2.14) Then, we have to introduce an auxiliary field B µ whose equation of motion yields the Bianchi identity dF = 0, which is locally equivalent to F = dA. Thus, we add the term B µ ∇ ν F αβ µναβ in the action and now the variables are F (instead of A) and B. Integrating that term by parts we can write In the case of pure gravity we do not need to include an explicit symmetrization in the µν indices because the terms P ρσγ µ R νρσγ and ∇ σ ∇ ρ Pµσνρ are symmetric [35]. However, in the case at hands we were not able to prove that the different structures appearing in Eµν are automatically symmetric -although we highly suspect it -and hence we have to symmetrize explicitly. 7 Let us note that it is possible to define more general notions of electromagnetic duality than the one we are considering here. A generic duality transformation establishes an isomorphism between two different theories, that is, a bijection between their configuration and solution spaces [92][93][94], but a remarkable property of certain theories -such as ungauged four-dimensional supergravity -is the fact that they admit a non-trivial subgroup (the U-duality group) of the duality group that does preserve the theory [95]. We thank Carlos S. Shahbazi for clarifying these points to us.
is the dual field strength. Now, as required, the variation with respect to B yields the Bianchi identity of F , and the variation with respect to F gives a relation between F and G, This is a direct expression of the dual vector field G , but now we have to invert this expression to obtain F , which in general is some complicated function of G and of the curvature, Inserting this back in (2.15) we get the dual theory, which on general grounds is different to the original theory. In appendix A we compute explicitly the dual theory in the case of a Lagrangian quadratic in F .

Mass, charges and thermodynamics
For the purposes of this paper, let us very briefly review the issue of conserved charges and black hole thermodynamics. In the case of electric and magnetic charges, these are obtained directly from the equations of motion. In particular, the generalized Maxwell equation can be written as If a current three-form J is placed on the right-hand-side, the equation implies that the current is conserved, dJ = 0, and hence the natural definition of electric charge is where the integration is taken at spatial infinity. Let us note that in asymptotically flat spacetimes, with a Lagrangian L = R − F 2 +higher-order, we have M → F asymptotically, so in practice we can replace M by F as long as the integral is performed at infinity. In the asymptotically AdS case, there is a theory-dependent constant c q such that M → c q F at the boundary of AdS. On the other hand, the magnetic charge is defined in the standard way, Gravitational conserved charges in higher-order gravities were studied in Refs. [96][97][98], but here we are only interested in the total mass. In the case of asymptotically flat spacetimes, it turns out that the mass can be formally computed using the same prescriptions as for GR, for instance via the ADM [99,100] or the Abbott-Deser [101] formulas. Thus, the mass can be computed by examining the asymptotic behaviour of the metric in the usual way, i.e., identifying the term 2GM/r ∈ g tt . In the AdS case, a global theory-dependent factor should be added to those formulas, corresponding to replacement of Newton's constant by the socalled effective Newton's constant G eff [102]. For L(R µνρσ , F αβ ) theories we do not expect these results to be affected, since the higher-order operators formed from F αβ decay too fast at infinity to contribute to the mass.
Regarding black hole thermodynamics, it is known that higher-curvature gravities minimally coupled to a Maxwell Lagrangian F 2 satisfy the first law of black hole mechanics [103,104] Here M , Q and P are the mass and charges computed as we specified above, T is the Hawking temperature of the black hole [105] and S is Wald's entropy [106,107], given by 8 where the integral is carried out on the bifurcation surface of the event horizon and µν is the binormal to this surface. In addition, Φ h is the electrostatic potential at the horizon, while Ψ h is the electrostatic potential of the dual vector field whose field strength is given by (2.17). These can be computed according to the following relations where ξ ν is the Killing vector that generates the horizon, with the condition that Φ and Ψ vanish asymptotically. It was shown in Ref. [110] that the form of the first law is unchanged when instead of the Maxwell Lagrangian one considers non-linear electrodynamics minimally coupled to Einstein gravity. However, such analysis does not include the case of non-minimally coupled terms, and one may wonder if the form of the first law could be altered in that case. This is, just like the entropy is no longer proportional to the area when there are higher-curvature terms, S = A/4, the question is whether the quantities Φ or Q in (2.22) could have to be replaced by different ones in order for the first law to hold. Another non-trivial question is whether the Noether's charge and the Euclidean path integral approaches yield equivalent results for black hole thermodynamics [111,112].

Static and spherically symmetric solutions
In this section we address the problem of finding static, spherically symmetric (SSS) solutions of L(R µνρσ , F αβ ) theories. Obviously, the equations of motion are far too general to be solved without specifying a Lagrangian, so instead our aim is to understand the structure of those equations and to describe a class of theories for which the problem can be simplified.
For a SSS configuration, we can write a general ansatz for the metric in the usual form which depends on two functions N and f , while in the case of the vector field we will assume either an electric or magnetic ansatz, as given by One may also consider dyonic vectors, but this increases the difficulty of the problem taking into account the non-linearity of Maxwell equations and that duality invariance is generically lost, so we will discuss purely electric or purely magnetic configurations only.

The reduced Lagrangian
The easiest way to find the equations of motion consists in evaluating the Lagrangian on this configuration so that we obtain a reduced Lagrangian, defined as for the electric and magnetic configurations respectively. Then, the equations of motion for the variables N , f and Φ in the electric case are given by 5) and similarly in the magnetic case, where δ/δN , etc, represent the Euler-Lagrange variation.
Using the chain rule, it can be shown that the equations E N = E f = E Φ = 0 are equivalent to some components of the equations of motion obtained from direct evaluation of (2.12) and (2.13). The rest of the equations are then satisfied on account of the Bianchi identities [43]. So far we have made no assumptions on the form of the Lagrangian L(R µνρσ , F αβ ), besides it being an invariant formed from the curvature, the metric, and the field strength of the vector field. In order to make further progress, we are going to assume that the Lagrangian is a polynomial in F αβ and R µνρσ , i.e, it is composed of terms of the form F 2m R n , with the indices contracted appropriately. The discussion proceeds differently for electric or magnetic configurations, so let us study both cases separately.

Electric solutions
Let us first consider the case of an electric vector field. By looking at the structure of the curvature tensor on a SSS metric [113], it is not difficult to show that a monomial composed out of 2m field strengths and n curvatures has the following structure when evaluated on such configuration, In addition, F has the property of being homogeneous of degree 0 in N , and hence it can be expressed as where the sum is always finite, and the functions F ij are polynomial in f , f and f . Now, schematically the reduced Lagrangian is L N,f,Φ = N r 2 sin θ n,m F 2m R n . Therefore, the Maxwell equation, obtained from variation with respect to Φ, reads where Q is an integration constant. Since the left-hand-side of the last equation is a (polynomial) function of Φ , we can in principle invert it so that we get 9 In this way we have eliminated one of the variables in the system of equations. Now we have to plug the value of Φ in the equations for N and f and we get two differential equations for these functions, Since Φ sol is generically a highly nonlinear (not even polynomial) function of the variables f and N and their derivatives, solving these equations is in general an inaccessible problem.

Magnetic solutions
Let us turn now to the case of magnetic configurations. It can be seen that the monomials F 2m R n have the form where F has the same structure as in the electric case. Then, the reduced Lagrangian has the following form,  We see that this equation is always solved for χ (θ) ∝ sin θ, and hence we have where P is the magnetic charge. Therefore, magnetic vectors are not affected by the corrections and they have the usual form. Now we can just plug this value of χ back in the equations of N and f . However, let us note that, since χ does not depend on these functions, in this case we can insert the on-shell value of χ on the reduced Lagrangian, from where we get a Lagrangian for N and f only. Thus, for magnetic configurations we can consider from the beginning F m = P dθ sin θ ∧ dφ , (3.15) and we can derive the equations for N and f from the effective Lagrangian where we are dropping the factor of sin θ since it is irrelevant. In this case, the equations are much simpler than for the electric vector field.
3.2 The condition g tt g rr = −1: Generalized Quasitopological theories Following the previous discussion, we have come to the conclusion that magnetic solutions are much simpler to study than electric ones. However, even if we consider restricting ourselves to these magnetic solutions, the equations for N and f are typically too complicated to obtain relevant information in general, so further simplification is desirable.
In the case of pure gravity, an interesting class of theories has been identified in recent years. These theories are known as Generalized Quasitopological gravities 10 (GQGs) [42,43] and they are characterized by possessing SSS solutions of the form i.e., with N = 1, and where in addition the equation of motion for f can be partially integrated. Other remarkable aspects of these theories are that the thermodynamic properties of black holes can be studied fully analytically and that they only propagate a massless graviton on constant curvature backgrounds. Interestingly, GQGs can be nicely combined with minimally coupled vector fields, since they respect the property of having solutions with g tt g rr = −1. However, one may wonder if it is possible to generalize these theories in the case of non-minimally coupled vector fields. The defining property of GQGs, from where all the rest can be derived 11 [42,43], is that their reduced Lagrangian becomes a total derivative when evaluated on the single-function 10 The works [42,43] can be considered as the ones establishing the general properties of the new family of theories, but these are motivated by earlier works on Einsteinian cubic gravity [49][50][51]. 11 See also chapter 3 of [114] for a refinement on some of the results in [43].
ansatz (3.17), so we may try to extend this definition when a non-minimally coupled vector is present. The main problem in that case is that the reduced Lagrangian L N,f,A will not in general enjoy the same structure as for pure gravities, since it will strongly depend on the electromagnetic potential A µ . Thus, for instance, it does not seem possible to impose that L 1,f,A be a total derivative without specifying the value of A µ . A more reasonable property to ask would be that the Euler-Lagrange equation of f vanishes identically when it is evaluated on N = 1 and on a gauge field that solves the Maxwell equation, Still, in the electric case we have seen that A sol typically has a non-polynomial dependence on f and its derivatives, so it seems difficult to find by brute force theories satisfying this property. Fortunately, the situation is different for magnetic configurations. In fact, we have seen that in the magnetic case we can work with a reduced Lagrangian that depends only on N and f , since we can set F = F m as in (3.15) from the start. Then, we can see that the reduced Lagrangian L N,f for non-minimally coupled theories F 2m R n with a magnetic vector field has the same structure as for pure gravity theories, and therefore one can extend straightforwardly the definition of Generalized Quasitopological gravities to these non-minimally coupled terms. In particular, we can present the following Theorem 1 Let us consider a theory with a Lagrangian L(R µνρσ , F αβ ) of the form L(R µνρσ , F αβ ) = R − F 2 + higher-derivative terms, (3.19) where the higher-derivative terms are formed from monomials of the Riemann tensor and the field strength, 12 schematically R n F 2m . Let us consider a SSS configuration given by the metric (3.17) and by a magnetic vector field with field strength (3.15) and let us define the reduced Lagrangian of the system as (3.20) If the Euler-Lagrange equation for the reduced Lagrangian L f vanishes identically, i.e., The points 1 and 2 follow immediately from the results in [43] taking into account that the reduced Lagrangian has the same structure as in the case of pure gravity. Point 3 also follows from the results there -see also [114] -, but it is somewhat trivial, since the terms of the form R n F 2m with m > 0 do not contribute to the linearized equations of the metric, while the pure curvature terms satisfying (3.21) are known to yield Einstein-like linearized equations. In addition, let us note that the degrees of freedom associated to the gauge field are the same as in Maxwell theory, since the generalized Maxwell equation for any theory of the form (3.19) is of second-order in A µ . As for point 4, it technically stands as a conjecture since no formal proof has been offered so far. However, we highly suspect a proof must exist due to the large evidence collected in the case of Generalized Quasitopological gravities [42,[50][51][52][53], and we also show in the next sections that it holds for all the non-minimally coupled theories that we construct.
So far, these results involve magnetically-charged black hole solutions. However, following the procedure explained in Section 2.2, one can dualize any theory satisfying (3.21) and obtain a new theory with electrically-charged solutions of the form (3.17). Hence all the items in the Theorem 1 hold for the dual theories after replacing "magnetically-charged" by "electrically-charged". This motivates the following definition of Electromagnetic Generalized Quasitopological Gravities (EGQGs): Definition 1 We say that a theory L(R µνρσ , F αβ ) belongs to the family of Electromagnetic Generalized Quasitopological Gravities (EGQG) if and only if its Lagrangian or the Lagrangian of its dual theory satisfies the condition (3.21).
One could refer to these theories as "Magnetic" and "Electric" GQGs, respectively, but we call them collectively Electromagnetic GQGs for two reasons. Firstly, because it makes sense to use the adjective electromagnetic to express that these theories are non-minimally coupled to an electromagnetic field. Secondly, because theories of one and another class are simply related by duality transformations, and hence they are equivalent. Let us also note that, even though the condition (3.21) can be satisfied by simple polynomial Lagrangians (see next sections), the dual (electric) theory is typically much more involved and will contain an infinite number of terms, justifying thus why Definition 1 was made in terms of the magnetic ansatz (3.15).
EGQGs, just like their pure gravity counterparts 13 , come in two different classes. In the general case, the equations of motion of these theories for charged SSS configurations are 13 There is an important qualitative difference between EGQGs and purely gravitational GQGs though.
While Lovelock gravities belong naturally to the GQG family, EGQGs do not include Lovelock-like theories in which a gauge field is non-minimally coupled to gravity, like the ones defined at Refs. [112,115,116].
reduced to a second-order equation for the function f . However, there is a special subset of these theories for which the order of this equation is reduced twice again, and we are left with an algebraic equation. In this case we say that the theories belong to the Quasitopological class (without the "generalized" part). For pure metric theories, Quasitopological gravities only exist D ≥ 5, but as we show below, infinitely many Electromagnetic Quasitopological gravities exist in D = 4.

Electromagnetic Quasitopological gravities
As we stated above, the family of theories allowing for single-function SSS solutions come in two classes: those for which the equation for f is algebraic are called "Quasitopological", while those for which f satisfies a 2nd order equation are called "Generalized Quasitopological". In this section we focus on the former case. In order to find theories of the EGQG class, one first writes down a general Lagrangian (including for instance all the densities of the form F 2m R n up to a given order), then evaluates the Lagrangian on the configuration given by (3.17), (3.15), and finally demands that the condition (3.21) is satisfied. At the end, that condition gives us a number of constraints on the couplings of the higher-derivative terms. In addition, if we want to restrict to Quasitopological theories, we must only keep the subset of those theories that yield an algebraic equation for f .
At lower orders in the derivative expansion, one can easily find all the theories of this type, but the process becomes more and more involved at higher orders, as the number of independent densities one can include grows very fast. Thus, a general analysis does not seem a priori accessible. Nevertheless, from the analysis of the lower-order densities we can probably extract a general conclusion on the structure of EQGs. Indeed, in the appendix B, we observe the following two facts. First, that there are only two Lagrangians of the form F 2 R belonging to the Electromagnetic Quasitopological class 14 . Second, that at order F 2 R 2 , despite the larger number of linearly independent invariants one can construct, there are only two different ways in which these densities modify the equation of f . Thus, if we are only interested in studying SSS solutions, it suffices to keep two representative EQ densities at a given order. By repeating an analogous analysis at higher orders, we observe that the situation seems to be general, in the sense that there are many independent EQGs but their equations on SSS metric are degenerate so that there are only two different contributions at every order. So, instead of studying the whole set of EQGs, we will provide a set of representative theories which -we conjecture -span all the possible modifications to SSS solutions. This will be enough for our goal, which is to study black holes in these theories.
It turns out that it is not difficult to provide a set of two representative Lagrangians of the EQ type at every order. In order to write them, let us introduce the following notation: 14 In this very particular case it even happens that the EGQ family coincides with the quasitopological one. This is of course not a general feature.
with the convention that R 0 µν σ] . Then, we have the following Lagrangians of order R n F 2m : so that these Lagrangians are well-defined for all integers n ≥ 0 and m ≥ 1. Let us also mention at this point that, in the case n = m = 1, the existence of Lagrangians with simple magnetic spherically symmetric solutions has been previously noticed in the literature [117][118][119], although, to the best of our knowledge, generalizations have not been constructed.
In order to show that the Lagrangians above belong to the EGQ family, we have to check that they become a total derivative when evaluated on the single-function ansatz (3.17) with a vector field strength given by (3.15). For that, we evaluate the Lagrangians on the general SSS metric ansatz given by (3.1). We get: where, following the notation of [113], we have introduced which represent several components of the Riemann and Ricci tensors. We can immediately check that both Lagrangians L n,m , defined back at (4.2) and (4.3), belong to the EGQ class. For that, we use the expressions (4.4) and (4.5) above and evaluate them on N = 1. A direct computation shows that n,m , (4.9) where I (a) n,m = 2 m+n P 2m d dr r 3−4m ψ n , (4.11) Since these Lagrangians are total derivatives, the corresponding Euler-Lagrange equations for the single-function SSS ansatz vanish identically, showing that they are truly EGQ theories. Actually, we can refine a bit more this statement by noticing that the equation of motion for the metric function f (r) is algebraic, so that both L n,m are of the quasitopological type. We check this behaviour in the following subsection.

Spherically symmetric solutions with magnetic charge
Let us consider the most general extension of the Einstein-Maxwell theory that can be constructed out of the Lagrangians L where (4.14) As explained above, L 0,1 = 0, so the usual Maxwell term is included in the action (4.13) as long as we set λ 0,1 = 1. Let us also clarify that we formally include an infinite number of terms in the action out of generality, but at any moment we may consider that only a finite number of couplings are non-vanishing.
Our goal in this section is to obtain magnetically-charged SSS solutions in this set of theories. We have already determined that the Maxwell equation is always solved by magnetic vectors with field strength (3.15) and that, due to the properties of the theory, the Einstein's equations allow for solutions with N = 1. Thus, we only have to determine the equation of the function f in the SSS metric (3.1), which can be obtained by evaluating the action on the SSS ansatz (3.1) (with the help of (4.5)), varying with respect to N , and then evaluating at N = 1. The resulting equation takes the form of a total derivative, dÊ(f, r)/dr = 0, and upon integration we obtain the algebraic equation where M is an integration constant that can be straightforwardly related to the mass of the solution. In this expression we have introduced two functions α n and β n that are given by where the coefficients depend on the magnetic charge and on the coupling constants of the theory as β n,m = 2 m+n−2 P 2m 2(n+m−1) 2(n − 1)λ n,m + (n 2 − 4n + 2 + m(−2 + 4n))γ n,m , (4.18) As we remarked earlier, the equation of motion for f (r) (4.15) seems to be the most general equation one can get for the set of all Electromagnetic Quasitopological gravities. This is, we suspect that any other EQG will only have the effect of changing the value of the coefficients α n,m and β n,m above. We explicitly show in appendix B that any Quasitopological theory built with two curvature tensors and two gauge field strengths indeed satisfies this property. Along with the metric and the (magnetic) field strength there is an additional physical magnitude of interest that we can find: the electric potential associated to the dual field strength. Indeed, if (g µν , F µν ) is a magnetic solution of a given theory, then (g µν , G µν ), where G µν is the dual field strength as defined in Equation (2.17), is a solution of the associated dual theory. In this latter theory, the potential of G µν will be electric. This electric potential will make its appearance in the subsequent first law of black hole thermodynamics (see Subsection 4.3), so it will be useful to compute it. For that, using the notation employed in Eq. (2.3), first we notice that where we have implicitly defined Imposing the field strength to be magnetic (3.15) and evaluating on the general SSS ansatz (3.1), we find that the dual field strength G, given by Equation (2.17), takes the form where we defined Amusingly, the latter can be explicitly rephrased as G = −Ψ (r)dt ∧ dr, Ψ being the electric potential. It takes the form

Explicit non-singular solutions in quadratic-curvature theories
In general, an explicit solution of Eq. (4.15) in a theory involving an arbitrary number n of Riemann curvature tensors is not available, since we would need to obtain the roots of an nth-degree polynomial. Therefore, for the sake of simplicity, let us analyze the solutions of theories that are only of second order in the curvature (n ≤ 2) but including an arbitrary number of gauge field strengths (m ≥ 1). In this case, the equation of motion for the metric function f (r) takes the form where we have taken into account that α 0 (r) = −β 0 (r). Since this is a quadratic polynomial in f one may solve the equation right away to obtain two possible solutions (4.25) Now, one can check that the solution f + is asymptotically flat and that it satisfies f + (r) = 1 − 2M/r + . . . when r → ∞. In addition, this solution reduces to the Reissner-Nordström one in the limit in which the corrections vanish → 0, where one has to take into account that α 0 → P 2 /r 2 while α 1 , α 2 , β 1 and β 2 vanish in that limit. On the other hand, f − has an exotic asymptotic behaviour, f − (r) ∼ −1/β 2 (r), and it does not have an Einstein gravity limit, so we will consider only f + as the physically relevant solution.
A very remarkable property of these solutions is that, in many cases, the curvature singularity at r = 0 is regularized by the higher-derivative corrections. In order to analyze the behaviour of the solution near r = 0, let us assume that we only include terms containing up to 2m c field strengths (so that 1 ≤ m ≤ m c ). Then, the functions α i and β i read and we can see that, except for certain fine-tuned values for the couplings, α 2 (r) and β 2 (r) are the dominant terms in the limit r → 0. Hence, from (4.25) we get Thus, whenever α 2,mc + β 2,mc > 0 and β 2,mc = 0 we have f + (r) ∼ 1 + O(r 2 ) near r = 0, implying that the geometry is regular there. Notice that this is quite remarkable, since we do not need to fine-tune the couplings, only demand that they satisfy a bound. In order for the solution to be globally regular we also have to make sure that there are no other singularities, e.g., the term inside the square root in (4.25) should not become negative, but again this is easily achievable (for instance, if all the couplings are positive). Let us note that the regularization is also possible if we only include linear curvature terms (n = 1) -in fact this has been previously observed in the literature in the case of F 2 R theories [118], although in that case the couplings must be related in a specific way. By analyzing the solutions of Eq. (4.15) near r = 0, one can see that the regularization at r = 0 can also be achieved in the general case in which we include terms F 2m R n of arbitrary order. The form of the equation (4.15) forces f (r) ∼ 1 + O(r 2 ) near r = 0 in most of the cases in a quite natural way. Therefore we conclude that, without the need of much tuning on the coupling constants, the magnetically-charged SSS solutions of the EQG theories (4.13) are singularity-free. 15 Notice that we have made not reference yet to black hole solutions, because, as in Einstein-Maxwell theory, not all the charged solutions are black holes. If the charge is too large compared to the mass, the solution does not have a horizon and in GR this means that we have a naked singularity. However, in our theories the gravitational field does not diverge, so horizonless regular solutions exist. In Fig 1a we show the profile of f (r) for a black hole solution while in Fig. 1b we show the one corresponding to a gravitating point charge. In both cases, the gravitational field is regular everywhere.
Another physical quantity of interest is the dual electric potential Ψ, which was calculated back in Equation (4.23). In the particular Electromagnetic Quasitopological theories we are  Figure 1: The metric function f (r) and the potential Ψ (in appropriate units) for two given particular sets of couplings, magnetic charge P and mass M . Note that Fig. 1a represents a black hole solution with an inner and outer horizon while Fig. 1b is an instance of a non-black hole solution. In both cases, couplings have been chosen so that Ψ is regular at r = 0, and we see that it vanishes in this limit.
considering, it turns out that the electric potential takes the form

(4.29)
While the geometry is generally regular, this is not the always case for the electric potential. If we want to have a regular potential at r = 0, not any set of couplings is allowed. Indeed, since around r = 0 the metric function f (r) can be approximated by f (r) ∼ Therefore regularity requires that α 0,m − A(α 1,m + β 1,m ) − A 2 (α 2,m + 1 2m (β 2,m − α 2,m )) = 0 for all m. This will happen for a certain subset of the whole moduli space of couplings, but it is a realizable feature, as it is shown in Fig. 1a and Fig. 1b. Thus, in these cases black holes and horizonless solutions have regular gravitational and electromagnetic fields everywhere. These solutions are generalizations of the example recently reported in Ref. [128].

Black hole thermodynamics
After describing the static spherically symmetric solutions of Electromagnetic Quasitopological gravities (4.13), in this subsection we focus on black holes and their thermodynamic description. One of our goals is to check that the first law of black hole mechanics holds in these theories and to identify the relevant thermodynamic potentials.
Let us begin with the solution given by (3.17) and assume the metric function f (r) has some zero for r ∈ R + . The black hole horizon would be consequently located at r h = max{r ∈ R + |f (r) = 0}. Using the equation of motion for f (r) (4.15), after evaluation on r = r h we find that From here we can solve for the mass M of the black hole and get We can also obtain the temperature, for which we first work out the derivative of (4.15) at Substituting the expression for the mass found at (4.32) and taking into account that the temperature T of the black hole is given by 4πT = f (r h ), we are left with .

(4.34)
Our next objective is the computation of the black hole entropy S, which is given by Wald's formula (4.36) Plugging this result into Eq. (4.35) and performing the integration on the angular variables, we find the black hole entropy: (4. 37) We see that the entropy is no longer just the black hole area divided by 4 and we have corrections. In particular, we check that the type (a) theory described by (4.2) does not introduce any corrections to the entropy, being just the type (b) theory (and more concretely, the term involving a Ricci scalar) (4.3) the one which causes deviations from the Bekenstein-Hawking result.
In order to check the first law of black hole mechanics we need to bear in mind that we have a magnetically-charged solution. Consequently, if we consider the dual theory, this magnetic solution will become an electric one after dualization. However, we know that electric solutions satisfy a first law which includes the associated electric potential. Therefore, taking into account that black hole thermodynamics remain unchanged under duality transformations, we conclude that the dual of any electric solution, which will be magnetic, will satisfy a first law including the aforementioned electric potential. Hence we can consider this argument backwards to justify that it is the dual electric potential the one entering in the first law of thermodynamics of these magnetic solutions.
From (4.23) we get the following value of the dual electrostatic potential evaluated at the horizon where at the same time we may use the expression for T in (4.34). At this point we have the quantities M , T , S and Ψ h expressed explicitly as functions of r h and P , and it is not hard to check using (4.32), (4.34), (4.37) and (4.38) that the following relations hold, where ∂ r h and ∂ P denote partial differentiation with respect to r h and P , respectively. When expressed this way, one may directly check that the following first law is satisfied. Hence we have shown that there exists a first law of thermodynamics for the Electromagnetic Quasitopological gravities described by the action (4.13) which holds exactly.
Despite the presence of non-minimally coupled terms, this result shows that the first law is formally unchanged, i.e., the effect of the charge appears through the standard term Ψ h dP , where Ψ h is the electrostatic potential at the horizon.
In order to complete our study of the thermodynamic properties of these black holes, let us compute the free energy. On the one hand, this can be defined from the rest of thermodynamic potentials as F = M − T S . On the other hand, F should be obtained from the on-shell evaluation of the Euclidean action according to F = T I E . Regarding the computation through the Euclidean action, we have to include an appropriate boundary term and suitable counterterms. Finding these boundary terms for higher-curvature theories of gravity is a highly non-trivial issue, e.g. [129][130][131][132][133], but nevertheless one can see that, whatever these terms are, they do not contribute to the on-shell evaluation of the action in the case at hands. On general grounds, we expect that the boundary terms will be proportional to the first derivative of the Lagrangian with respect to the curvature [58], and since we are considering asymptotically flat situations, all such terms decay too fast at infinity to make a finite contribution. 16 Thus, we may use as a boundary term the standard Gibbons-Hawking-York term [134,135] minus its background contribution. Therefore we propose the following Euclidean action where we have already Wick-rotated the time coordinate. In order to evaluate this action on our black hole solutions, let us note that the on-shell Lagrangian takes the form of an explicit total derivative when evaluated on the single-function metric (3.17). This follows from (4.9), (4.10) and from a similar property satisfied by the Ricci scalar. Thus we have (4.43) and the Euclidean action reads Then, one can check that the evaluation of I(r) at r → ∞ is exactly cancelled by the boundary terms, so we are left with the evaluation at the horizon, I E = βI(r h )/4. This yields the following value (4.45) By comparison with (4.41) we check that F = I E /β, and consequently, we show the inner consistency of these computations and that the Noether-charge and that the Euclidean path integral approaches to black hole thermodynamics give equivalent results. 16 The situation is different for asymptotically AdS solutions, but in that case one may introduce an effective boundary term which is proportional to the Gibbons-Hawking-York term [58]. This procedure is known to work at least for theories of the GQG class and has been tested in several occasions [44,59,67].
Finally, let us work out the specific heat C P at constant magnetic charge. A direct application of the inverse function theorem shows that We check that C P generally vanishes in the extremal limit T → 0, which we study in more detail next.

Extremal black holes
We finish the study of black holes in Electromagnetic Quasitopological theories by pinpointing some characteristic features of their extremal limit. For simplicity, we restrict ourselves to the particular class of theories which are quadratic in the vector field strength (m = 1). If we define the dimensionless parameter ρ = r h / , then we can express the black hole mass M for this class of theories as where we have introduced the function 2 n ρ 2n+1 λ n,1 .  From here it is trivial to solve for P 2 and obtain the following extremal charge-to-mass ratio: . (4.50) Thus, we have an explicit a non-perturbative expression for the extremal charge-to-mass ratio in terms of the radius. Note that if only a finite number of terms is included in the action, the function U is a polynomial in 1/ρ. However, if an infinite number of them is added, U can actually be any function of the form U (ρ) = u(ρ −2 )/ρ, where u(x) is an arbitrary analytic function (to recover Einstein-Maxwell theory at low energies it must satisfy u(x) → 1 when x → 0, though). In Fig. 2 we show P/M | ext as a function of the mass for several choices of this function. The effect of higher-derivative corrections on extremal black holes has a particular interest in the context of the weak gravity conjecture (WGC) [76]. In fact, a mild form of the WGC states that the extremal charge-to-mass ratio in a consistent theory of quantum gravity must not decrease as the mass decreases. Thus, P/M ext must be a growing (or constant) function when when we move from larger to smaller masses. This condition ensures that the decay of an extremal black hole into a set of smaller black holes is possible, at least from the point  Fig. 2a we notice the existence of two branches (just the upper one would be connected to the Reissner-Nordström black hole), and the charge-to-mass ratio is monotonically decreasing with the mass for both branches. However, there are no extremal black holes below the mass at which both branches merge. On the other hand, in Fig. 2b we consider two different functions U (ρ) and we realize that in both cases the extremal charge-to-mass ratio is monotonically growing.
Our study, on the other hand, is fully non-perturbative, and we can analyze what happens when the corrections become important.
According to the WGC, just the theory depicted at Fig. 2a would be admissible, since it satisfies (for the two branches) that the charge-to-mass ratio decreases when the mass grows. However, note that extremal black hole solutions cease to exist below a minimal mass (when the two branches merge). Although this might seem to be a peculiar feature of the particular model considered, this behaviour turns out to appear quite generally. Additional examples of this situation are shown in Section 5.3 for a different family of theories. One should wonder what happens with the evaporation process of black holes at this point. For that, let us consider an initially large (M >> ) non-extremal black hole. Due to Hawking radiation, it loses mass until it approaches extremality. At that moment, it also needs to lose charge in order to continue evaporating, and this is achieved if the WGC holds by emitting a particle with charge-to-mass ratio p/m ≥ 1. Note that since our black holes satisfy P/M ext > 1 and this quantity becomes larger for smaller black holes, evaporation is not obstructed. In addition, a process by which an extremal black hole decays into a set of smaller, non-extremal black holes would be in principle allowed in terms of energy and charge conservation. Through this process, the black hole evaporates down to arbitrarily small masses, following approximately the line of extremal black holes in Fig. 2a (we may assume the black hole remains near-extremal during the evaporation process). Then, it will reach the minimal mass in order for extremal black holes to exist, and this can have several meanings. One possibility is that below that line there are no black holes at all, e.g., all the solutions are naked singularities or horizonless smooth configurations. The black holes could then transition to one of these objects, but this is highly speculative. A more interesting possibility is that below that mass any black hole is non-extremal, and this is precisely the case with the model depicted in Fig. 2a, corresponding to λ 1,1 = −2, λ 2,1 = 1/4. The black hole solutions of that theory with M = 5 (below the minimal extremal mass) are shown in Fig. 3. We see that, no matter how large the charge is, the black hole is non-extremal. Thus, in this case, the black hole can always lose mass by means of Hawking radiation without imposing any conditions on kind of particles emitted. The fact that there is no obstruction to achieve the black hole's evaporation is in fully agreement with the spirit of the WGC. Finally, let us comment on another interesting property that we can study in the extremal limit, namely, the value of the electrostatic potential Ψ at the event horizon. Indeed, in the limit T → 0 the quantity Ψ h takes the following surprisingly simple expression: Let us remark that Ψ h = 1 for extremal Reissner-Nordström black holes, but this is no longer the case for our black holes. In Fig. 4 we represent the electrostatic potential for the same three theories we considered in Fig. 2. We observe that in general such electric potential does not necessarily monotonically increase or decrease with the mass. We discover a rather counter-intuitive fact: rapid decreases in the charge-to-mass ratio plots seem to correspond to increases of the electric potential. Indeed, one would expect that as the charge diminishes, the potential to decrease as well, but we are finding precisely the opposite behaviour. This phenomenon not only happens for the theories forbidden by the WGC, but it also takes place for theories which, a priori, would be allowed. This is explicitly seen for one branch (the one disconnected from Reissner-Nordström) of the theory with U (ρ) = −1/ρ + 4/ρ 3 − 1/ρ 5 .

Electromagnetic Generalized Quasitopological gravities
We have just studied two families of EGQGs belonging to the quasitopological subset, i.e., they yield an algebraic equation for the metric function f allowing for the analytic study of the black hole solutions. However, these are not the only type of EGQGs that exist. Analogously to the case for pure gravity, there are some theories for which f does not satisfy an algebraic equation, but a 2nd order differential equation -these are the proper Generalized Quasitopological theories. One may wonder why study these theories since we have already described two infinite classes of theories with simpler black hole solutions. There is a good reason, though. It turns out that, even though for Generalized Quasitopological theories we usually are not able to provide an explicit black hole solution, we can, nevertheless, obtain the thermodynamic properties of black holes exactly. As we show below, the thermodynamic relations can have a quite different form with respect to the quasitopological case, so that these new theories provide us with qualitatively different modifications of the Reissner-Nordström solution.
Based on our previous experience with (purely gravitational) GQGs, we expect that there are many of these theories at each order, so we will not attempt to provide a complete classification of this family of theories. Instead, our goal is to show that these theories indeed exist and to study some of their properties. A general characterization of EGQ theories may be addressed elsewhere.
We have found a simple family of EGQ Lagrangians, which read where (R n ) µν is the n-power of the Ricci tensor, with the convention that R 0 µ ν = δ µ ν . Let us first of all show that these Lagrangians indeed satisfy the GQG condition, as given by Eq. (3.21). Evaluation on the the general magnetic SSS ansatz, yields the following reduced Lagrangian, This expression will be useful to compute the equations of motion later. Further evaluation on N = 1 shows that the Lagrangian becomes a total derivative, and therefore it belongs to the EGQG class. Let us now study the black hole solutions of these theories.

Black holes
Let us consider an extension of Einstein-Maxwell theory with the terms of the EGQ class (5.1) above: Here is an overall length scale while µ n,m are dimensionless couplings. Note that the usual Maxwell term −F 2 is included in the sum, since we have L EGQG 0,m = −(4m−3)(F 2 ) m . Thus, by convention we set µ 0,1 = 1. By construction, this theory has magnetically-charged solutions of the form (3.17), and hence we only have to determine the equation of motion for f . This can be done most easily by taking the variation of the reduced Lagrangian L N,f with respect to the function N and evaluating at N = 1. The result is a third-order equation for f which takes the form of a total derivative: Note that all the terms in the sum are of second order in derivatives except for those with n = 0 and n = 1, since in those cases the Lagrangians reduce to minimally coupled terms (F 2 ) m and to Electromagnetic Quasitopological theories as the ones studied in the Section 4, respectively. Now, integrating the equation above we get where M is an integration constant that, as we show next, turns out to be the mass. Since in general the equation above is a differential equation of second order, there are two more integration constants which need to be fixed by the boundary conditions. This is analogous to the case of purely gravitational GQ theories which has been studied in various papers [42,[51][52][53], so let us just comment briefly on it. First, we impose that the solution is asymptotically flat (we do not have a cosmological constant), which implies that f (r) → 1 when r → ∞. In the asymptotic region, we may expand the general solution in the following form: where f p is a particular solution while f h represents a deviation with respect to that solution (and will satisfy a homogeneous equation). We can obtain a particular solution by assuming a 1/r expansion, which yields the following result On the other hand, the boundary conditions imply that f h → 0 asymptotically, and hence we can assume that it is arbitrarily small. Thus, plugging (5.13) into (5.12) and expanding linearly in f h we get where the asymptotic expansion of the coefficient reads a = 2µ 2,1 4 P 2 The equation above can be solved in terms of Bessel functions, but for our purposes it suffices to note that the asymptotic solution behaves as where A and B are integration constants. Thus, when µ 2,1 > 0 one of the modes is exponentially growing and the other one is exponentially decaying. By setting the appropriate constant to 0 we achieve an asymptotically flat solution with a free integration constant. When µ 2,1 < 0, the solutions become highly oscillating at infinity and the only way to obtain a regular solution is to set A = B = 0, thus there are no further boundary conditions that one can fix. This is problematic because we cannot impose regularity at the horizon (see below), and therefore there are no regular black hole solutions in this case. Thus, we only consider µ 2,1 > 0. If this coefficient is 0, the constraint will appear in the next coefficient in the expansion.
On the other hand, we impose the existence of a regular horizon, i.e., a point r h at which f (r h ) = 0 and around which f is analytic. In particular, we assume that f has a Taylor expansion of the form where we are making explicit that f (r h ) = 4πT , where T is Hawking's temperature. When we insert this expansion into the equation (5.12), we get a system of equations that relate the coefficients a n . Nonetheless, the first two equations are special, since they only involve r h and T . These read These two equations allow one to get (implicitly) the temperature T and the radius r h once M and P are given. The rest of the equations provide relations for the coefficients a n . A simple inspection reveals the only free parameter in the expansion is a 2 , and the rest of the a n are fixed in terms of it. Finally, a 2 is fixed by demanding that the solution be asymptotically flat. The full solution f (r) can be obtained by a numeric integration of (5.12) using (5.18) as initial condition, and implementing a shooting algorithm to search for the value of a 2 that yields the correct asymptotic behaviour. Such numeric resolution will be carried out elsewhere, but comparing with previous works on neutral black holes in GQ theories, we expect that the solution exists providing the condition on the couplings discussed above is satisfied. Fortunately, a great deal of information about these black holes can be obtained without resorting to the numeric solution.

Black hole thermodynamics
Even though the profile of the solutions has to be determined numerically, one remarkable property of the theories in Eq. (5.9) -which is shared by all the theories of the GQ class -is that the thermodynamic properties of black holes can be found analytically. First, note that the two relations (5.19), (5.20) above give us the relation between M , P and T . Unfortunately, such relation cannot be written explicitly due to the complicated form of the equations, but it is nevertheless possible to solve the system of equations parametrically. Let us first introduce two dimensionless parameters p and x defined as Then, we can define the 2-variable function In terms of these quantities we can rewrite (5.19) and (5.20) as follows Whenever ∂ x W = 0, we can obtain explicitly r h (x, p) from the second equation. On the other hand, if W does not depend on x, the same equation determines the relation x(p), while r h is free. Note that this only happens in the trivial case in which the higher-order Lagrangians do not depend on the curvature, and hence it is not relevant for our purposes. Then, inserting r h (x, p) in (5.23) we obtain the explicit relation M (x, p), and analogously, we get T (x, p) and P (x, p) from (5.21), namely, (5.25) Thus, we have been able to write all these thermodynamic quantities, as well as the radius, in terms of two independent parameters x and p. This is a useful way to study the thermodynamic phase space of these theories. Let us now compute the rest of thermodynamic properties of these black holes.
The entropy is computed by Wald's formula, as introduced previously in Section 4 When computing the derivative with respect to the curvature of the Lagrangians (5.1), each time we derive one of the Ricci tensors R αβ we end up generating a contraction between the binormal αβ and a field strength. Note that such contractions are always 0 for magnetic configurations, since and F are orthogonal in that case. Thus, only the action of the derivative on the Ricci scalar appearing in (5.1) yields a non-vanishing contribution. The result reads where we have already performed the integration on the horizon. Evaluating this expression at r = r h and using the parameters (5.21) and the function (5.22), we get the simple result Again, using (5.24) we obtain the explicit relation S(x, p).
Let us now compute the electrostatic potential at the horizon. As we saw in Section 2.2, the dual field strength is given by (2.17). The derivative of our Lagrangians (5.1) with respect to the field strength yields where Z ρσαβ = R n−1 ρσ (nRg αβ − (4n + 4m − 3)R αβ ) . Additionally, we can obtain the free energy from the on-shell Euclidean action. The computation can be done using the same prescription for the boundary terms as in Eq. (4.42). The bulk action can be evaluated right away thanks to our on-shell Lagrangians being total derivatives (5.8). Then, the evaluation at infinity gets canceled with the contribution from the boundary terms and we are left with the evaluation of the quantities I n,m in (5.8) at the horizon. This yields the following result for the free energy, F = I E /β: Summarizing, the equations (5.23), (5.25), (5.28), (5.33) and (5.34), together with the relation (5.24), give us explicit expressions for all the thermodynamic quantities M , P , T , S, Ψ h , F and the radius r h in terms of two independent parameters x and p. The theorydependence of all these formulas is encoded in the function W defined in (5.22). Let us now check that these quantities satisfy consistent thermodynamic relations. In particular, they should satisfy the 1st law of black hole mechanics, This relation is in fact verified. The easiest way to see this consists in assuming first that r h is an independent variable in the expressions of M , S and P (Eqs. (5.23), (5.28) and (5.25), respectively). Then, the variations of those quantities with respect to just x and p automatically satisfy the first law above. Afterwards, we may take the variation only with respect to r h (assuming that now x and p are independent variables) and we check that it also satisfies the first law once we notice the constraint (5.24). Hence when the dependence of r h on x and p is taken into account, the first law holds too for arbitrary variations of the free parameters.
On the other hand one can also check that F = M −T S, which is a non-trivial consistency test of our results, indicating that the Wald's entropy (Noether charge) and the Euclidean action approaches are equivalent.

Extremal and near-extremal black holes
Let us study how the corrections affect extremal black holes. In terms of the variable x, the extremality condition T = 0 implies that x and r h are related according to (5.36) Due to this, (5.24) becomes a complicated equation that determines the relation between p and x at extremality. Namely, we have To simplify the discussion, let us consider the subset of theories that are only quadratic in the Maxwell field strength (but which have an arbitrary number of higher-curvature terms). In such case, the function W has the form For this function, it is possible to solve (5.37) explicitly to obtain p(x) at extremality: Then, from (5.23) and (5.25) we obtain the mass and the charge 17 40) and the extremal charge-to-mass ratio The entropy in turn reads Then, as we did in Section 4.4 we can check some particular cases to see if it is possible to satisfy the mild form of the Weak Gravity Conjecture at a non-perturbative level. Since x = 2 /r 2 h , we must demand that P/M ext is monotonically growing with x, although in that case we also have to make sure that M is a decreasing function of x. As an example, let us consider the case in which there is a single higher-derivative term in the action so that U = 1 + µ n x n . Then we have M ext = √ x 1 + (n + 1)µ n x n 1 + (2n + 1)µ n x n , P M ext = 1 + (2n + 1)µ n x n 1 + (n + 1)µ n x n , (5.43) For µ n > 0, the extremal charge-to-mass ratio and the mass are actually monotonically decreasing with x, so this case should be discarded according to the WGC. On the other hand, if we take µ n < 0 we observe that for small x (large M ) the charge-to-mass ratio is in fact growing with x. However, it soon reaches a maximum value and then decreases again. Moreover, M has a minimum value, so there are no extremal black holes below certain 17 We stumble upon the following fact: the results for the extremal mass and charge given by Eq. (5.40) coincide exactly with those for EQs -that one may obtain from Eqs.   Fig. 5. We find the same qualitative behaviour in all of these cases, namely, P/M has a maximum value which happens for the minimum mass. Thus, it seems quite difficult for P/M ext to be a growing function all the way down to M = 0, at least within this family of theories. Nevertheless, this can be interesting from the point of view of the WGC, since, as we saw in Section 4.4, it may imply that below the minimal mass all the solutions are non-extremal black holes, and hence there is no obstacle to prevent the evaporation of these black holes.
Another intriguing fact about these examples is that the corrections to the extremal entropy are negative. For instance, in the case of U = 1 + µ n x n we have S ext = π 2 x 1 + (4n + 1)µ n x n 1 + (2n + 1)µ n x n < πP 2 if µ n < 0. (5.44) This seems in contradiction with some claims and results in the literature [77,82] which relate positive corrections to the extremal charge-to-mass ratio with positive corrections to the entropy. However, the contradiction is not such, since, as noted in Ref. [84], the comparison must be done with the corrections to the near extremal entropy, while the corrections to the extremal entropy are independent. This is an example of that situation.  Figure 6: Temperature vs mass diagram at fixed charge for near-extremal black holes with ∂ 2 M/∂T 2 P ext < 0. The extremal black hole is not the state with the minimal mass. We consider the model U (x) = 1 + x 2 , but nevertheless the profile of the curve will be similar for any other case in which ∂ 2 M/∂T 2 P ext < 0.
Finally, let us also briefly comment on near-extremal black holes. A characteristic property of extremal black holes in Einstein gravity is that the specific heat at constant charge goes to zero, while its first derivative is positive. This means that near-extremal black holes satisfy M − M ext = cT 2 with c > 0, and therefore the mass of the black hole grows as we increase the temperature. Interestingly enough, this is not always the case for our black holes with higher-derivative corrections. The specific heat, defined as C P = ∂M ∂T P , vanishes at extremality, but its first derivative reads instead One can see that this quantity can have either sign, depending on the model and on the value of x. If it is positive, then near-extremal black holes behave as in Einstein-Maxwell theory and they are stable, in the sense that when we increase the temperature (hence we depart from extremality) the mass also increases. The case in which this quantity is negative is quite intriguing. It implies that in order to get away from extremality, the black hole must lose mass. Therefore, extremal black holes are thermodynamically unstable and they do not represent the minimal mass state for a given charge. Instead, the minimal mass state will take place at a different point in which C P = 0, and this is the solution to which the black hole tends when it evaporates. An example of this situation is represented in Fig. 6, where we show T vs M at fixed charge for a particular set of higher-derivative terms. Another consequence of this effect is that, in a certain region of the parameter space, there exists more than one black hole solution with the same mass and charge. This non-uniqueness of solutions can be thought as a discrete violation of the no-hair conjecture and is analogous to the situation with charged black holes in Einsteinian cubic gravity that was recently reported in Ref. [66].

Conclusions
In this paper, we have introduced a new class of non-minimally coupled higher-derivative extensions of Einstein-Maxwell theory. These theories are characterized by possessing magnetic SSS solutions characterized by a single metric function f (see (3.17)) whose equation of motion is (at least partially) integrable. In addition, within this set of theories, the thermodynamic properties of black holes can be computed exactly. Such theories are analogous to the Generalized Quasitopological gravities and thus we refer to them as Electromagnetic Generalized Quasitopological gravities. As in the case of pure gravity, we have seen EGQGs come in two main classes: those for which the SSS equations of motion can be reduced to an algebraic equation for f belong to the "quasitopological" class, while if the equation is of second order we say that the theory is properly of the "generalized quasitopological" class. We have constructed an infinite number of densities of both types, although we suspect that there are many others, especially in the case of Generalized Quasitopological theories. Determining the most general structure of these Lagrangians would be an interesting problem.
In the case of Quasitopological theories, we have shown some explicit examples of black hole and non-black hole solutions -see Section 4.1. We observed that, in a quite remarkable and general way, these solutions possess globally regular geometries, i.e., the timelike singularity at r = 0 characteristic of charged black holes or point charges is smoothed away by the higher-derivative corrections. In slightly more restrictive cases, we showed that the electrostatic potential of the dual theory also remains finite everywhere, thus making these solutions particularly appealing. In particular, in the horizonless case, one may regard these objects as solitons or even as four-dimensional fuzzballs.
For both, Quasitopological and Generalized Quasitopological theories, we have performed a detailed study of black hole thermodynamics -see Sections 4.3 and 5.2. We have been able to provide explicit expressions for all the relevant thermodynamic potentials and we have shown that the first law of black hole mechanics, holds exactly. Here S is Wald's entropy and Ψ h is the electrostatic potential of the dual theory evaluated on the event horizon. Thus, the first law is formally unchanged with respect to the case of a minimally coupled gauge field. In addition, we have checked that the Euclidean methods provide the same answer for black hole thermodynamics than the Noether's charge approach. In particular, we have seen that the on-shell Euclidean action yields indeed the free energy, T I E = F = M − T S. Due to the large space of theories that we consider, we have not made a general analysis of the features of the new thermodynamic relations, so a more detailed study is left for future work.
Motivated by the weak gravity conjecture, we did study the properties of extremal and near-extremal black holes in these theories. A mild form of the WGC states that, in a consistent quantum theory of gravity, the charge-to-mass ratio of extremal black holes should grow monotonically as the mass decreases. This would allow for the decay of extremal black holes in terms of energy and charge conservation. Previous literature had studied perturbative corrections to the extremality bound in a variety of theories, ranging from general EFTs to Stringy effective actions [77][78][79][80][81][82][83][84][85][86]. Although our theories do not belong (a priori) to those categories, they have the advantage of allowing us to perform exact, non-perturbative computations. Thus, they may be used to learn about the corrections to extremality at large coupling. As we observed in Sections 4.4 and 5.3, it is always easy (e.g., by choosing the signs of the couplings appropriately) to get P/M ext to satisfy the WGC when the mass is large (i.e., in the perturbative regime). However, when the curve P/M ext vs M ext is continued to lower masses, one often finds that it stops at a minimal mass, meaning that there are no extremal black holes below that mass. There can be different reasons for this behaviour, but we have shown with an example (see Fig. 3) that a possibility is that below that mass all solutions are non-extremal black holes, regardless the value of the charge. This is actually appealing from the point of view of the WGC, since it implies that, below the minimal mass, charged black holes find no obstruction to evaporate.
Higher-derivative corrections can introduce new effects into the game, and in particular we observed another situation that has implications for black hole evaporation. In some instances -as shown in Fig. 6 -it may occur that the extremal black hole is not the one with a minimum mass for a given charge. In those cases, (near-)extremal black holes are unstable, and tend to decay to this minimum mass black hole, which has a non-vanishing temperature. Thus, in that situation one does not have to worry about the charge-to-mass ratio of the extremal black hole, but about the one of the minimum mass black hole. An analogous example has been recently reported in Ref. [66] in the context of Einsteinian cubic gravity with a (minimally-coupled) Maxwell field.
The new theories offer various possibilities since they allow us to perform many explicit computations that are inaccessible in general higher-derivative theories. Thus, let us close the paper by commenting on future directions. As we have already mentioned, it would be interesting to complete the characterization of EGQ Lagrangians to find the most general action of this type. On the other hand, here we have focused on asymptotically flat solutions, so one could extend this work by including a non-vanishing cosmological constant. The asymptotically anti-de Sitter case is particularly relevant due to its connection to holography. In fact, it is known that higher-derivative gravities with a negative cosmological constant are very useful holographic toy models that can be used to learn non-trivial information about Conformal Field Theories -see Refs. [59,67] for recent results involving Generalized Quasitopological gravities. Since EGQGs contain higher-derivatives not only of the metric but also of a vector field, these may be used to probe additional aspects of a CFT.
One could also characterize subsets of these theories satisfying additional properties. For instance, we expect that some of the EGQGs allow for single-function Taub-NUT solutions whose thermodynamic properties can be studied exactly, as the ones in [44]. In particular, some theories of the quasitopological subclass might allow for explicit Taub-NUT solutions. In addition, non-minimally coupled electrodynamics may be of interest in cosmology -see e.g. [136] -so we may wonder if some of these theories could be useful in that context, as the ones in [45][46][47].
Regarding higher-dimensional generalizations, we find several possibilities. Since the defining property of EGQGs is related to the structure of their magnetic SSS solutions, a straightforward generalization can be achieved in the case of gravity coupled to a (D − 3)form in D dimensions. In that case we ask for the same property to hold, namely, that there are magnetic SSS solutions with g tt g rr = −1. Note that the dual of these theories corresponds to gravity coupled to a vector field, and the magnetic solutions become electric ones, so they are more natural in this dual frame. On the other hand, in higher-dimensions there are solutions in the form of extended objects, such as black strings. One may search for higher-derivative theories in which the structure of these solutions remains simple, in a way analogously to the condition g tt g rr = −1. As an example, one may consider the case of a sixdimensional theory with a non-minimally coupled 2-form field, and search for magnetic black string solutions satisfying such type of condition. These higher-dimensional generalizations may be addressed elsewhere.
Taking into account Eq. (2.17), we have that Therefore the action I dual dual to (A.1) turns out to be where we have included the contribution from the Lagrange multiplier term which imposes the Bianchi of F µν . Define now a tensor Q −1 µνρσ with the same symmetries as Q µνρσ which satisfies Q −1 µνρσ Q ρσαβ = δ αβ µν . (A.6) Using this inverse tensor of Q, it is clear that Therefore, on taking into account Eq. (A.4), we find that the dual theory may be expressed in the following compact form: where we have defined the tensor χ µνρσ as This procedure can be used to obtain electrically-charged solutions from E(G)Qs with magnetic ones, as it was explicitly done in Ref. [128]. In general, the dual theory (A.8) will not be polynomial. However, it is possible to write it as a formal power series if we decompose Q as Q µνρσ = δ µνρσ +Q µνρσ . (A.10) This decomposition is rather natural since it corresponds to an expansion around the pure Maxwell term given by δ µνρσ . The inverse Q −1 can in turn be written as where we have defined (Q 0 ) µνρσ = δ µνρσ and (Q n ) µνρσ =Q µ 2 ν 2 µνQ µ 3 ν 3 µ 2 ν 2 · · ·Q ρσ µnνn . Substituting back in Eq. (A.8), we obtain the dual action in terms ofQ: In this appendix we are going to build all Electromagnetic (generalized) quasi-topological gravities constructed by linear combinations of terms up to quadratic order both in the Riemann curvature tensor R µνρσ and in the gauge field strength F µν . For that, we first classify all E(G)Qs RF 2 , made up of scalar terms with one Riemann and two field strength and, afterwards, we proceed analogously for E(G)Q theories R 2 F 2 , whose constituent terms contain exactly two Riemann tensors and two field strengths.

B.1 RF 2 theories
Our first task is to find a set of diffeomorphism-invariant terms which span all possible scalars of the form RF 2 . This can be done straightforwardly and we find the following basis of invariants containing one Riemann tensor and two field strengths: Now we build the Lagrangian density For the theory to be of the (Generalized) Quasitopological type, we must ensure that the latter expression vanishes. This is accomplished by where M is an integration constant appropriately chosen to be identified with the mass, as done in the main text. There are two important conclusions to extract from Eq. (B.6). Firstly, we recognize precisely the same structure as in Eq. (4.15), if we limit ourselves to the Einstein-Hilbert term and the terms with n = 1, m = 2. Secondly, we check that the set of EGQs and EQs coincide for theories of the form RF 2 , since Eq. (B.6) is algebraic. This property does not hold generally of course and is very particular of RF 2 theories. As a matter of fact, the special properties of these Lagrangians had been previously noticed in the literature [117][118][119].
The theory is a (Generalized) Quasitopological one if all A i vanish simultaneously. Such a system of linear equations is solved by: b 3 = − 2b 1 − 2b 10 , 7b 7 = −40b 1 + 40b 10 − 4b 11 − 2b 15 − 18b 2 − 4b 4 − 4b 5 + 8b 6 , 7b 9 = −80b 1 − 8b 12 − 4b 13 − 2b 14 + 2b 15 − 36b 2 − 8b 4 − 8b 5 − 96b 6 + 14b 8 . (B.11) Imposing these constraints, one obtains the generic expression for any EGQ constructed out of terms with two Riemanns and two field strengths. However, we still need to figure out which of these EGQs are actually quasitopological. For that, we must learn when the equation for f (r) is algebraic. As aforementioned, this equation is derived after evaluating the Lagrangian (B.8) on our magnetic SSS ansatz, varying the subsequent action respect to N and afterwards setting N = 1. One obtains: and we recognize the same structure as in Eq. (4.15) after restricting ourselves to those terms with n = 2, m = 1. Hence we have proven that the equation for f (r) of the most general EQ constructed from terms with at most two Riemann tensors and two gauge field strengths is indeed represented by Eq. (4.15), after an appropriate choice of couplings.