{'fraction_non_alphanumeric': 0.03677001534434516, 'fraction_numerical': 0.01276062821554292, 'mean_word_length': 4.455681398498092, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermodynamic analysis of black hole solutions in gravitating nonlinear electrodynamics

We perform a general study of the thermodynamic properties of static electrically charged black hole solutions of nonlinear electrodynamics minimally coupled to gravitation in three space dimensions. The Lagrangian densities governing the dynamics of these models in flat space are defined as arbitrary functions of the gauge field invariants, constrained by some requirements for physical admissibility. The exhaustive classification of these theories in flat space, in terms of the behaviour of the Lagrangian densities in vacuum and on the boundary of their domain of definition, defines twelve families of admissible models. When these models are coupled to gravity, the flat space classification leads to a complete characterization of the associated sets of gravitating electrostatic spherically symmetric solutions by their central and asymptotic behaviours. We focus on nine of these families, which support asymptotically Schwarzschild-like black hole configurations, for which the thermodynamic analysis is possible and pertinent. In this way, the thermodynamic laws are extended to the sets of black hole solutions of these families, for which the generic behaviours of the relevant state variables are classified and thoroughly analyzed in terms of the aforementioned boundary properties of the Lagrangians. Moreover, we find universal scaling laws (which hold and are the same for all the black hole solutions of models belonging to any of the nine families) running the thermodynamic variables with the electric charge and the horizon radius. These scale transformations form a one-parameter multiplicative group, leading to universal"renormalization group"-like first-order differential equations. The beams of characteristics of these equations generate the full set of black hole states associated to any of these gravitating nonlinear electrodynamics...

and the horizon radius. These scale transformations form a one-parameter multiplicative group, leading to universal "renormalization group"-like firstorder differential equations. The beams of characteristics of these equations generate the full set of black hole states associated to any of these gravitating nonlinear electrodynamics. Moreover the application of the scaling laws allows to find a universal finite relation between the thermodynamic variables, which is seen as a generalized Smarr law. Some particular well known (and also other new) models are analyzed as illustrative examples of these procedures.

Introduction
In General Relativity the theorems on singularities, on the uniqueness of the black hole solutions of the Einstein vacuum field equations, and the Cosmic Censorship conjecture [1] lead to the strong belief that the Kerr spacetime describes accurately the outcome of the gravitational collapse of rotating electrically neutral matter distributions. When such distributions support a nonvanishing total electric charge, similar arguments, now involving the minimally coupled Einstein-Maxwell theory, lead to Kerr-Newman black hole configurations as the final state of the collapse [2].
However, there are several arguments justifying the study of other BH configurations 1 , which are solutions of more complex gravitating NEDs, as possible final states of the gravitational collapse of electrically charged matter. A first argument comes from the low energy regime of string theory, where fundamental excitation is gravity. If one accepts the fundamental character of this theory, then the dynamics of the emerging gauge fields in this regime is governed by non-minimal effective Lagrangians [3]. In this context some nonlinear electrodynamics, as the Born-Infeld model [4], are to be seen as more "fundamental", in describing the dynamics of the electromagnetic fields, than the usual Maxwell theory, which is regarded, within this picture, as a weak field approximation. Besides that, Born-Infeld generalizations of SU (2) non-Abelian gauge field theory, coupled to tensor-scalar gravity, have been also used for the description of dark energy in a cosmological context [5].
A second argument concerns the quantum vacuum effects on the electromagnetic field, which are neglected in the classical Einstein-Maxwell theory. Now accepting the fundamental character of the Maxwell theory and its quantum formulation, the corrections of the Dirac vacuum on the dynamics of electromagnetic fields in flat space can be described, at a phenomenological classical level, through nonlinear effective Lagrangians that define NEDs, the first historical example being the EH Lagrangian [6]. It accounts for the nonlinear effects of the Dirac vacuum on the low energy electromagnetic wave propagation, calculated to lowest order in the perturbative expansion. When higher order operators are included in the perturbative calculation we are led to a sequence of effective Lagrangians which take the form of polynomials in the field invariants, arranged as an expansion in operators of increasing dimensions [7]. These effective Lagrangians exhibit finite-energy elementary solutions in flat space, even though the bare (Maxwell) Lagrangian does not [8]. When minimally coupled to gravity such Lagrangians describe satisfactorily the vacuum effects in the slowly varying curvature regions [9]. Disregarding the rotational degrees of freedom, which are not essential for the present considerations, the elementary solutions of any physically admissible NED (a concept to be defined below) minimally coupled to gravity, are asymptotically flat, exhibit a curvature singularity at the center and their electrostatic field strengths are monotonically decreasing functions of the radius [10,11]. Near the singularity this effective approach fails, and the description of the vacuum effects would require a consistent quantum treatment of both the electromagnetic and gravitational fields, which is not available at present. But in moving away from the singularity the curvature stabilizes very quickly and becomes slowly varying in the characteristic distances of the vacuum screening. Consequently, the effective approach remains accurate everywhere in the manifold, excepted in a small region close to the singularity. Therefore, in describing the collapsed state of non-rotating electrically charged matter, the BH configurations associated to some NEDs (as EH) could be more "realistic" than the Reissner-Nordström one.
A general analysis of the gravitational configurations and geometric structures of the elementary solutions associated to the whole class of admissible NEDs, minimally coupled to gravity (including some well known models such as Maxwell, Born-Infeld or Euler-Heisenberg, as particular representative examples) has been performed in Refs. [10,11]. Nevertheless, the issue of the thermodynamic properties and laws of the corresponding BH solutions was not considered in those works. Such a study is the main object of this paper.
The thermodynamic analogies for asymptotically flat BH solutions of the Einstein equations were established in the early seventies for the cases of the Kerr-Newman BH and its non-rotating and uncharged limits [12]. The subsequent discovery of Hawking's radiation [13] turned these analogies into actual laws of the thermodynamics of these configurations. Since then, an important amount of literature has dealt with this issue in depth. In this way many papers have analyzed the thermodynamics of the Reissner-Nordström black hole (see e.g. [14] and references therein, also [15]). Let us mention some extensions and developments going beyond the analysis of this solution: i) the validity of the thermodynamic laws has been established for asymptotically flat BH solutions of several classes of field theories minimally coupled to gravity, including NEDs supporting asymptotically coulombian elementary solutions [16]. In this way the thermodynamic analysis has been performed for asymptotically flat BH solutions associated to some particular NEDs, as the gravitating BI model [17] or, more recently, the family of NED models defined through Lagrangian densities which are constructed as powers of the Maxwell one [18]; ii) this analysis has been also performed for BH solutions of the gravitating Born-Infeld NED in asymptotically cosmological backgrounds, as (anti-)de Sitter spaces in (3 + 1) dimensions [19,20], in higher dimensions [21,22] and in (2 + 1) dimensions [23], extending the Hawking-Page results of the uncharged case [24] and motivated by the AdS/CFT correspondence [25]; iii) several extensions of these solutions, resulting from the minimal coupling of NEDs to gravitational dynamics defined beyond the Einstein-Hilbert action, have been considered, as in the cases of Gauss-Bonnet and Lovelock [26], f (R) [27] or scalar-tensor theories [28]; iv) statistical computations of the entropy for the extreme and near extreme five dimensional RN black holes have been performed in the framework of string theory [29].
However, most of the aforementioned works deal with some particular models, while a comprehensive study of the thermodynamic properties of BH configurations associated to general gravitating NEDs is possible, but still lacking. The aim of this paper is to provide such a study, by analyzing the full set of NEDs minimally coupled to gravitation in (3 + 1) dimensions and defined by Lagrangian densities which are arbitrary functions of the two field invariants, but restricted by some consistency requirements (such as the positivity of the energy, the regularity of the Lagrangian and the parity invariance) leading to what we shall call admissible theories. In Ref. [11] the gravitating admissible NEDs were classified into two sets, according to the form of the asymptotic behaviours of their elementary solutions. These solutions, which are always asymptotically flat, can approach flatness as the Schwarzschild field (normal case) or slower than it (anomalous case). In the anomalous case some thermodynamic variables (such as the ADM mass) cannot be defined and the thermodynamic analysis is not possible, at least in the usual way. Consequently, in the present work we shall deal only with the analysis of the thermodynamics of asymptotically normal BH solutions supported by the corresponding families of G-NEDs.
In section 2 we summarize some results of Ref. [8] regarding the elementary solutions of NEDs in flat space, which are necessary for the present study. In particular, we introduce the classification of these theories that, when generalized to the gravitating versions, defines the different families of possible BH solutions.
In section 3 we briefly describe the geometrical structures characterizing the BHs associated to each family. This will be necessary in the next section for a clear understanding of the thermodynamical behaviours. A detailed report of the contents of this section has been published in Ref. [10]. Section 4 is devoted to the thermodynamical analysis itself. First, we prove the validity of the expression of the first law for any admissible gravitating NED supporting asymptotically normal BH solutions. In particular this validity is established for those models with non Maxwellian weak field limit, for which the theorems of Ref. [16] do not immediately apply. Next we analyze the behaviour of the main thermodynamical variables and specific heats for the BH configurations associated to the different families. In particular the study of the thermal properties of the extreme BH configurations suggests the possible presence of van der Waals-like first-order phase transitions between these states for models belonging to several families. The results of this analysis are illustrated with some well known examples (as the RN, EH and BI) representative of different families, while other examples are specifically built to exhibit some particular behaviours.
In section 5 we find some universal relations between the thermodynamic variables of the asymptotically-normal BH states associated to admissible G-NEDs. First, we obtain a universal expression for the generalized Smarr law, valid for all such BH solutions. This law reduces to the usual Smarr formula in the RN case and to the corresponding expressions obtained for the few particular examples of G-NEDs studied in the literature. Next, we obtain scaling laws satisfied by the functions relating the thermodynamic variables, which are universal (i.e. valid for any G-NED). We make a preliminary discussion on some implications of these laws. In particular, we show the existence of a one-parameter multiplicative group structure underlying the scaling transformations and obtain the associated first-order differential equations ("renormalization group"-like equations) relating the derivatives of several state variables with respect to the charge (Q) and the horizon radius (r h ) of the BH solutions. These equations are the same for any G-NED and the associated characteristic curves in the corresponding three-dimensional space define the running of the thermodynamic variables and generate the "equations of state" for the sets of BH solutions associated to any gravitating NED (defined as surfaces supported by the bean of characteristics).
Finally, we conclude in section 6 with some comments and perspectives.
2 General results on NEDs in flat space-time Let us summarize the results concerning the complete characterization and classification of the whole set of admissible NEDs through the properties of their ESS solutions [8], which will be necessary for the analysis of the gravitating problem. The dynamics of these fields in flat four-dimensional space-time are governed by Lagrangian densities which are single-branched functions of the two gauge field invariants is the field strength tensor of the vector potential A µ , F * µν = 1 2 ǫ µναβ F αβ its dual and E and H the electric and magnetic fields, respectively. For physical consistence of the corresponding theories these functions must be defined in a open and connected domain of the X − Y plane, including the vacuum, and be restricted by the conditions of parity invariance (ϕ(X, Y ) = ϕ(X, −Y )) and positivity of the energy. The latter condition reads [8] ρ ≥ with vanishing vacuum energy (ϕ(0, 0) = 0). Moreover, the ESS solutions E(r) must be single-branched, defined and regular for r > 0. This requires ϕ(X, 0) to be a strictly monotonically increasing function of X in the range of values X = E 2 (r > 0) (∂ϕ/∂X | Y =0 is required to be continuous and strictly positive there). The above set of requirements define what we have called admissible NEDs, the only ones that we shall consider here. In addition one can assume ϕ(X, Y ) to satisfy the condition (for X = E 2 (r > 0), Y = 0) which is necessary and sufficient for the linear stability of the ESS solutions, as shown in Ref. [8].
The Euler field equations associated to the Lagrangians (1) read and the corresponding energy-momentum tensor takes the form If we look for ESS solutions of the field equations (E(r) = E(r)r/r, H = 0) we find the first integral where the integration constant Q is identified with the point-like charge source of the solution. Once the form of ϕ(X, Y ) is specified Eq.(6) gives in implicit form the expression of the ESS field E(r, Q) as a function of R = r/ |Q|, in such a way that it scales as |E(r, Q)| = |E(R, Q = 1)|.
As we shall see below, this scaling law induces similar laws in the different thermodynamic functions describing the systems of BH solutions of any gravitating NED. We shall call "characteristic" the Q = 1 configurations, from which the whole thermodynamic structure of the set of BH solutions associated to a given G-NED is determined through these scaling transformations (see sections 4 and 5 below). The positivity of the derivative of the Lagrangian function in Eq.(6) allows to restrict the analysis to the case (E > 0, Q > 0) without loss of generality. Assuming behaviours for the ESS solutions around the center and asymptotically of the form the admissibility conditions restrict the values of the exponents to p ≤ 0 and q < 0. At the center the fields diverge for p < 0 and behave as for p = 0. In the latter case the parameter a (the maximum field strength) and the exponent σ > 0 are universal constants for a given model, whereas the coefficient b(Q) is related to the charge of each particular solution as b 0 = b(Q = 1) being also a universal constant of the model. Asymptotically E(r → ∞) vanishes for the ESS solutions of any admissible model. The corresponding behaviours of the Lagrangian densities around these values of the fields (X = E 2 (r)) are given by the expression where α i , (i = 1, 2) are positive constants which are related with the coefficients in Eq. (8). For X → ∞ (and p < 0) we have whereas for X → 0 we have In both cases γ i > 1/2. If p = 0 the Lagrangian densities behave around the center (X = E 2 (r = 0) = a 2 , Y = 0) as if σ = 2. If σ = 2 this behaviour becomes In these equations ∆ = ϕ(a 2 , Y = 0) < ∞ if σ > 2, and it is a universal constant of the model. If σ ≤ 2 the value of ∆ can be calculated from the explicit expression of the Lagrangian density (see Eqs.(81) and (82) below). Consequently, for σ ≤ 2 the Lagrangians exhibit a vertical asymptote at X = a 2 . For σ > 2 they take finite values there, with divergent slope (the Born-Infeld model [4] belongs to this case with σ = 4). Equation (5) leads to the expression for the energy density distribution of the ESS solutions. As easily seen, the spatial integrals of these functions around the center of the ESS solutions converge if −1 < p ≤ 0 and diverge if p ≤ −1. This allows a classification of the admissible NEDs in terms of the central behaviour of their ESS solutions, leading to three families. We have called A1 the first family, corresponding to the range of values −1 < p < 0, for which the ESS solutions diverge as r → 0, but the integrals of energy converge there. For these models the Lagrangian functions behave, when X → ∞, as in Eq.(11) with γ > 3/2. We have called A2 the second family, for which p = 0, corresponding to the NEDs supporting ESS solutions which are finite at the center. Obviously, the integrals of energy converge there for these solutions (see Eqs. (14)- (16)). The third family corresponds to p ≤ −1 and supports ESS solutions diverging at the center and whose integrals of energy diverge also there (we have called these models UVD, acronym for ultraviolet divergent). When the integral of energy is convergent at the center (cases A1 and A2) the "internal energy function" giving the energy of the ESS field of a point-like charge Q inside a sphere of radius r, is well defined. The asymptotic behaviour of the ESS solutions allows for a second classification of the admissible NEDs. For large r the integrals of energy converge asymptotically if q < −1 and diverge if −1 ≤ q < 0. Thus we have defined four classes of models differing by this behaviour: B1 models (when −2 < q < −1) whose ESS solutions are asymptotically damped slower than the Coulomb field and the integrals of energy converge asymptotically; B2 models (when q = −2) whose ESS solutions behave asymptotically as the Coulomb field; B3 models (when q < −2) whose ESS solutions vanish asymptotically faster than the Coulomb field; IRD models (acronym for infrared divergent) when −1 ≤ q < 0, whose ESS solutions are asymptotically damped slower than the Coulomb field and their integrals of energy diverge asymptotically. When the integral of energy converges asymptotically (cases B1, B2 and B3) the "external energy function" giving the energy of the field outside the sphere of radius r, is well defined. The combination of these central and asymptotic behaviours allows for the exhaustive classification of the admissible NEDs in twelve classes. Six of them (combinations of A1 and A2 with B1, B2 and B3 behaviours) support finite-energy ESS solutions, which are to be interpreted as non-topological solitons if the necessary and sufficient stability condition (3) is fulfilled by the corresponding Lagrangian densities in the domain of definition of these solutions [8]. For these six families the convergence of the integral of energy in the whole space leads to the relations for the energy functions associated to the ESS solutions of charge Q. Using Eqs. (6) and (16) it can be easily seen that these energies, when they are well defined, scale as where the total energy of the characteristic configuration ε(Q = 1) is a universal constant for a given model. By using the last expression of ̺ in Eq. (16) for the calculation of (18) as r → 0, after some partial integrations and taking account the first integral (6), as well as the central and asymptotic behaviours of the fields in these finite-energy cases, we are led to the useful expression relating the total energy of the soliton with the integral of the field. The ESS solutions of the remaining six cases (combinations including UVD, IRD or both behaviours) are energy divergent. When coupled to gravity all the models of the twelve families lead to consistent gravitating field theories. However, we shall see that the families involving the IRD behaviour fail to allow for a thermodynamic analysis in the usual way.
In figure 1 we have plotted the profiles of the Lagrangian density functions of the admissible NEDs (for Y = 0) around the vacuum and on the boundary of their domain of definition.

General results on gravitating NEDs
The full set of admissible NEDs minimally coupled to gravity, and the structure of their G-ESS solutions have been classified and extensively analyzed in Refs. [10,11]. However, as mentioned in section 1, this analysis concerned mainly the geometric structure of these solutions and excluded the study of their thermodynamic properties. Before tackling this thermodynamic problem let us summarize some pertinent results of these references.
As usual, the action corresponding to the minimal coupling of gravity and any NED reads Fig. 1 Qualitative behaviour of the admissible Lagrangian densities ϕ(X, Y = 0) i) for small ESS fields (E 2 = X ∼ 0, corresponding to the B and IRD asymptotic behaviours of the ESS solutions), ii) for large ESS fields (E 2 = X → ∞, corresponding to the A1 and UVD centralfield behaviours) and iii) for finite maximum field-strength models (X = E 2 max = a 2 ), corresponding to the A2 central-field behaviour. The Lagrangian densities behave as X γ 2 around X = 0 (in the B and IRD cases) and as X γ 1 around X → ∞ (in the A1 and UVD cases) where the γ i constants are related to the central and asymptotic behaviours of the ESS solutions through Eqs. (12) and (13) respectively. In the A2 cases the Lagrangian density exhibits a vertical asymptote at X = a 2 (if σ ≤ 2) or takes a finite value with divergent slope there (if σ > 2). In the intermediate range of X values, matching the central and asymptotic regions, ϕ(X, Y = 0) must be strictly monotonically increasing, for admissibility.
g and R being the metric determinant and the curvature scalar, respectively. The associated Einstein equations for the G-ESS solutions can be written in an adapted Schwarzschild-like coordinate system, with an interval of the form where dΩ 2 = dθ 2 + sin 2 (θ)dψ 2 is the angular metric and g 00 = g(r) is the only metric function to be determined. The simplicity of the form of this interval is a consequence of the Einstein equations together with the equality of the mixed components of the energy-momentum tensor, T 0 0 = T 1 1 , which remains the same as in flat space in spherical coordinates. These equations reduce to whose compatibility can be easily established. The general solution of this system reads where C is an arbitrary integration constant. The characterization and classification of the admissible NEDs in flat space, summarized in section 2, allows for the complete characterization and classification of the corresponding G-NEDs. The asymptotic behaviour of the metric function is dominated by the behaviour at large r of the ESS solutions in flat space and is independent of their properties at finite r. Indeed, when the asymptotic behaviour of the ESS solutions in flat space belongs to the cases B1, B2 or B3 the metric function can be written as where ε ex (r, Q) is the external energy function, defined by (18), and the constant M is identified as the ADM mass. Owing to the asymptotic behaviour of the external energy function in these cases (∼ r q+1 , q < −1), the last term in (26) can be neglected at large r and the metric function behaves asymptotically as the Schwarzschild one [10]. The combinations of these asymptotic B-cases with the central-field A1 and A2 cases, which support finite-energy ESS solutions in flat space, admit an alternative expression for the metric function in terms of the internal energy function, which reads where the constant U is now related to the gravitational mass and the electrostatic soliton energy in flat space through This constant may be roughly interpreted as the gravitational binding energy of the configuration plus the proper energy of the point-like source (see, however, the discussion about this kind of interpretations and the associated difficulties in Ref. [30] and the pertinent references therein). In reference [10] we have called "critical" the configurations with U = 0. As we shall see in the next section, in the finite field amplitude cases (A2) with the particular value of the charge the critical configurations correspond to extreme black points [31], which exhibit a rich and peculiar thermodynamical behaviour. With Q < Q c they correspond to non-extreme BPs. In A2 cases with Q > Q c (or for A1 cases with any Q), the U = 0 configurations correspond to BH with an external event horizon and a vanishing-radius Cauchy inner horizon. For models belonging asymptotically to B-cases and with UVD centralfield behaviour (e.g. Maxwell Lagrangian) the electrostatic energy diverges and only the expression (26) for the metric function remains valid. For NEDs with IRD solutions in flat space the G-ESS solutions are asymptotically flat, but not Schwarzschild-like (asymptotically anomalous behaviour) and the ADM mass cannot be defined [11]. In these cases the metric function takes the form (27) if the central-field behaviour belongs to cases A1 or A2 (but now the interpretation (28) of the constant U makes no sense). For the families with UVD-IRD behaviour, only the expression (25) describes properly the form of the metrics, which are thermodynamically meaningless.
At finite r the geometric structure of the G-ESS solutions is characterized by the central behaviour of the electric field in flat space. The event horizons (if any) are located on the largest solution of the equation g(r) = 0. For models with UVD central-behaviour solutions the metrics exhibit a structure similar to that of the Reissner-Nordström solution (naked singularities, extreme black holes or two-horizons black holes). If the central behaviour belongs to cases A1 or A2 we have, in addition, new kinds of solutions: non-extreme singlehorizon black holes (for families A1 and A2) and extreme and non-extreme black points (for families A2) [10].
All the G-ESS solutions of the Einstein equations minimally coupled to admissible NEDs must belong to one of the twelve classes obtained as combinations of these asymptotic and central-field behaviours. Consequently, these results, which are summarized in table 1, characterize qualitatively the possible geometric structures of the G-ESS solutions of all admissible G-NEDs.

Thermodynamic analysis
Let us consider now the thermodynamic properties of the G-ESS black hole solutions of the admissible NEDs. In Ref. [16] D. A. Rasheed established the validity of the zeroth and first laws of black hole thermodynamics for BH solutions of large classes of field theories minimally coupled to gravity, including the stationary solutions of G-NEDs with Maxwellian weak field limit. In the present context the rotational degrees of freedom are absent and the Rasheed proof is restricted to the non-rotating ESS black hole solutions, corresponding to our asymptotically coulombian cases B2. As we shall see, it is easy to extend the validity of this first law for non-rotating solutions including the cases B1 and B3 (the zeroth law is trivially satisfied for the ESS black hole solutions). However, as already mentioned, in the asymptotically IRD cases, for which the ADM mass is not defined and (as we shall see) the electrostatic potential diverges at infinity, the thermodynamic relations cannot be established, at least in the usual way. Consequently, the following considerations concern the B cases only.
The set of G-ESS (non-rotating) black hole solutions associated to a given G-NED can be assimilated to "states" of a system characterized by two independent parameters, the most immediate ones being the integration constants: the mass M , well defined for asymptotically-B NEDs, and the charge Q. The location of the black hole event horizons is the largest solution r h (M, Q) of the equation obtained from the condition g(r h , Q, M ) = 0 on Eq. (26). Owing to the admissibility requirements, ε ex (r, Q) is a monotonically decreasing and concave function of r at constant Q, vanishing asymptotically as r q+1 with q < −1 in the B cases [8], [10]. It exhibits a vertical asymptote at r = 0 in case UVD or takes, in cases A1 and A2, a finite value there (= Q 3/2 ε(Q = 1); see Eq. (20)) with divergent slope in case A1. In case A2 we distinguish three qualitatively different subcases, according to the values of the charge: case A2a if Q > Q c , case A2b if Q < Q c and case A2c when Q = Q c . The geometrical structures associated to these three subcases have been extensively discussed in Ref. [10]. The behaviour of M as a function of r h for several representative values of Q, obtained from Eq. (30), is plotted in the figures 2 and 3 for the A2 and A1-UVD cases, respectively. The locations of the minima of these curves (when present) correspond to the horizon radii (r he ) and masses (M e ) of extreme black holes, determined by equation (30) and the extremum condition The cut points between the constant-mass horizontal lines and the positiveslope branches of the curves of constant charge, define the event horizon radii and the masses of the different black hole configurations. For a fixed Q, if M lies under the corresponding curve, there are no horizons and the solutions are naked singularities. Otherwise, besides extreme BHs, we have two-horizon (Cauchy and event) BHs, single horizon BHs (in cases A1 and A2) and extreme (in case A2c) and non-extreme (in case A2b) black points. Moreover, since the radii of the external (event) horizons of the possible black hole configurations lie always on the positive-slope region of the constant-Q curves, the mass-radius relation at fixed Q is a monotonically increasing function for all the black hole solutions of the admissible G-NEDs. It approaches the Schwarzschild relation M = r h /2 at large r h , as well as at vanishing Q. Let us define other important "state functions" and establish the validity of the first law for the ESS black hole solutions associated to the set of all admissible NEDs of asymptotic types B.
From its general definition [32] the surface gravity (κ) for these charged BHs takes the form where Eqs. (16), (18) and (26) have been used. If we look now for the total differential of the "energy" M as a function of r h and Q in Eq. (30), and using again Eqs. (16) and (18) we obtain for the partial derivatives the expressions and where A 0 (r h , Q) is the covariant time-like component (in the coordinate system (23)) of the four-vector potential in the appropriate gauge, which is fixed through the conditions Thus, the total differential of M reads With the usual definitions of the "temperature" T and the "entropy" S in terms of the surface gravity κ and the horizon area A, given by together with the usual notation for the electric potential, the differential expression (36) takes the more familiar form This is the expression of the first law, which holds for the black hole solutions of any G-NED of the nine families defined by the combinations of the three Basymptotic behaviours and the three central-field behaviours UVD, A1 and A2.

The extreme black holes
Let us now analyze the important case of the extreme BHs. From Eq.(31) and the first of (37) it is obvious that the temperature of these objects vanish.
The condition (31) does not involve the mass and then, defines the radius of the extreme BH configurations as a function of the charge (r he (Q)), which is monotonically increasing. Indeed, the general expression of the derivative of this function can be obtained from (31) and the first integral (6) in terms of the values of the field on the extreme-BH horizon. It reads and is always positive for admissible models (for Q > 0 solutions). The explicit form of r he (Q) which, as we shall see at once, determines the whole thermodynamic behaviour of the subclass of extreme BHs, depends on the particular model and can be directly obtained from Eqs. (6), (16) and (31), once the form of the Lagrangian is specified. Alternatively, it can be obtained by integrating Eq.(40) with the appropriate boundary conditions. In fact this integration can be formally performed using the scale properties (7) leading, after some manipulations, to the general expression where the function Ω is defined by the integral Q 0 and x 0 being related integration constants. For the Einstein-Maxwell theory this relation reduces to the simple and well known expression for Reissner-Nordström extreme BHs. But for general NEDs it can become largely involved.
Nevertheless, the r he − Q relation for small and large radii is typical of each family. In both limits it can be determined by replacing the expression (11) (with i = 1, 2) in (31) and using the central or large-r behaviours of the field (E(r → 0) ∼ r p ; E(r → ∞) ∼ r q ). This calculation leads to relations of the generic form where γ i , (i = 1, 2) are the exponents defined in Eqs. (12) or (13) and the λ i are given by the α i being the coefficients in Eq. (11).
In the A2 cases the small-radius extreme BHs are present for charges slightly over the critical one (Q Q c , see Fig.2) and the charge-radius relation in this region can be obtained by a similar procedure from Eqs. (9) and (31), leading to if σ < 2, if σ > 2 and when σ = 2. Once this r he − Q relation is obtained, Eq. (30) gives the mass of the extreme BHs as a function of their radius or their charge. The same will arise for the other state functions defined so far, fully determining the thermodynamic properties of the set of extreme BH configurations associated to any admissible NED. The lines of extreme BHs in the corresponding state diagrams are the T = 0 isotherms which, as we shall see at once, define basic boundaries separating different phases in these diagrams and allow for the search of first-order phase transitions.

The state variables and phase diagrams
As already mentioned, the elementary solutions of a given gravitating NED are characterized by two "state" parameters and it is possible to divide the plane of these parameters into separate regions corresponding to different phases of the set of solutions. The most immediate choice for this characterization is the pair of integration constants M − Q. It is possible to determine the qualitative form of the phase diagram for the different families of G-NEDs in terms of these parameters. To this end, let us analyze the behaviour of M as a function of the charge at fixed r h , which can be obtained from the scaling law 2 (30) and the last of Eqs. (20). This function is plotted in figure 4 for the cases A1-A2 (main frame) and UVD (small frame). The gravitating BI and RN models, respectively have been taken as representative examples for this figure, but the phase diagrams for the models belonging to the corresponding families are qualitatively similar. The curves representing the extreme BH states are also displayed on both frames of this figure. They are the envelopes of the beams of constant r h curves (see Eq. (31)) and begins tangent to the r h = 0 curve at Q = Q c in case A2. In case A1 the corresponding curve is obtained by moving this tangency point to the origin (Q c → 0, in the limit of large maximum field strength, a → ∞). In cases A1 and A2, the region over the r h = 0 curve contains the single-horizon BH states, the region between this curve and the extreme BH curve contains the two-horizon BH states and the region under these curves and over the Q axis corresponds to naked singularity configurations. Similar comments can be done for the UVD case (small frame), where the beam of constant-radius curves approach the M axis as r h → 0.
As a function of r h at fixed Q the behaviour of the mass has been already discussed at the beginning of this section and plotted in Figs. 2 and 3.
Phase diagrams involving other state variables can be obtained by similar methods. Let us consider the case of the field on the horizon E(r h , Q). As a function of r h , at constant Q, its form is obtained from the first integral (6). The behaviour of this function around r h = 0 and as r h → ∞ has been completely determined in section 2 for the elementary solutions of all classes of admissible NEDs and is plotted qualitatively in figure 5. It should be stressed that the only requirement to be satisfied by the field as function of r at constant Q, in all admissible cases, is to be monotonically decreasing, but its second derivative can change sign many times at finite r, leading to less smooth lines than the dashed ones plotted in this figure (for an example, see the upper small frame in figure 12 below). In determining the behaviour of E(r h , Q) as a function of Q at fixed r h we can use the scaling law (7), which leads to the condition For Q > Qc this curve splits the regions associated to the single horizon and two-horizons BH configurations. The diagram corresponding to the A1 cases is similar but with vanishing critical charge (Qc = 0). The analysis of the diagram in the small frame (UVD cases) is similar, but now the r h → 0 curves approach the Q = 0 axis (which corresponds to the Schwarzschild BH limit) and the charged single-horizon BH configurations are absent.
and we see that E, at constant r h , is an increasing function of Q. Equations (7) and (49) allow to perform a complete analysis of the behaviour of E around the limit Q = 0 (R h → ∞) and as Q → ∞ (R h → 0) for the different families of NEDs considered here. This behaviour is plotted in figure 6.
For the extreme BHs we can obtain easily the behaviours of the field on the horizon at small and large Q by replacing the r he − Q dependencies obtained for the different families. The large Q behaviours and the small Q behaviours in the cases UVD and A1 are obtained from Eq.(43) and have the same form, which reads as Q → 0 (i = 1) and as Q → ∞ (i = 2). In the same way, for the A2 cases the behaviour of E e (Q) for small-radius extreme BHs (Q → Q c ) is straightforwardly obtained from equations (44)-(47) in models with different values of σ. For example, for σ > 2 we have Qualitative behaviour of E(r h , Q) as a function of r h at fixed Q (main frame). At the center, in cases UVD and A1 the field diverges, while it takes a finite value in case A2. In this latter case the slope of E(r, Q) at vanishing r h depends on the value of σ (see Eq. (9)). The small frame exhibits the corresponding behaviours of the potential Φ(r h , Q). The dashed pieces of the lines, matching the central and asymptotic field behaviours, are always monotonically decreasing but can exhibit more involved forms than those showed here for many models (see the upper small frame in Fig.12 below). The thermodynamically pertinent parts of these curves start at the EBH radii r he (Q), obtained by solving Eq. (31) for the corresponding value of Q.
Another important phase diagram can be obtained from the external energy function ε ex (r h , Q). As a function of r h , for fixed Q, it has been extensively discussed and employed in Ref. [10] for the characterization of the different families. Its typical qualitative behaviour is plotted in figure 7 for all cases. The dependence of ε ex on Q at constant r h can be determined from the scaling law (20) and Eq. (34). This leads to the curves plotted in figure 8, corresponding to the cases A1-A2 and UVD (small frame). The representative region for the event horizons of the BH configurations in the A2 cases is bounded by above by the line of extreme BH configurations and the line r h = 0, which intersect in a point with Q = Q c , corresponding to the extreme black point configuration A2c. In case A1 the line of extreme BHs lies always under the line r h = 0 and defines the upper boundary of the BH region. This A1 case is obtained from the main frame of Fig.8 in the limit of divergent maximum-field strength (a → ∞ or Q c = 0). For the UVD case (in the small frame) this boundary is also given by the line of extreme BHs. The Q axis, corresponding to r h → ∞, is the lower boundary of this thermodynamically pertinent region.

4.3
The potential Φ(r h , Q) and van der Waals-like transitions.
As seen from the first law (39) the electrostatic potential is the conjugate state variable of the charge. Consequently the search for van der Waals-like phase  as a function of Q at constant r h we first obtain its scaling law from Eqs. (7) and (35), which reads (R h = r h / √ Q) and shows that Φ is a homogeneous function of r h and √ Q of degree one. From this law we obtain the Q-derivative of Φ as which shows that Φ is a monotonically increasing function of Q at constant r h . The form of this derivative and of the scaling law allows straightforwardly for a complete discussion of the behaviour of the potential at small and large Q for constant r h . This behaviour is displayed in figure 9 for the different families. In the A2 case, for Q < Q c , the curve r h = 0 corresponds to the non-extreme BPs and defines a part of the upper boundary of the thermodynamically pertinent region. The Q = Q c point of this curve defines the extreme BP. For large values of Q the potential at constant r h behaves as in the A2 cases. In the A1 cases we have whereas in the UVD cases we have (p ≤ −1). At small Q the behaviour of the constant r h curves is dominated by the large-r behaviour of the electrostatic field, and reads where −1 > q > −2 (B1), q = −2 (B2) and q < −2 (B3).
For Q > Q c the thermodynamically pertinent region lies under the curve of extreme BHs, which is defined by eliminating r he between equation (31) and the expression of Φ(r h , Q). The shape of this curve, which begins on the extreme BP (Q = Q c ) in the cases A2, and at Q = 0 in cases A1 and UVD, depends of the particular model, but its behaviour around these points and also for large values of Q is universal for every family. It can be determined explicitly from equations (43) to (47) and the corresponding expansions of Φ(r, Q) around r → 0 or r → ∞, obtained from the integration of Eqs. (8) and (9). For large values of Q these equations lead to where the positive coefficient is (q < −1). We see that the potentials on the horizon for large-Q extreme BH configurations diverge slower than √ Q in the B1 cases, approach a constant value 8π/α 2 (the same in all the asymptotically coulombian (B2) cases if the normalizing α 2 constants of the Lagrangians in Eq.(11) are fixed to the same value) and vanish in the B3 cases. For small Q values in the A1 cases the electrostatic potential on the horizon of extreme BHs behaves as In the UVD cases the same expressions (58) and (59) hold, with the replacements q → p, α 2 → α 1 , γ 2 → γ 1 and µ 2 → µ 1 (excepted in cases with p = −1, where logarithmic dependence arises, but does no modify the behaviours described here). Now, Φ e (Q) vanishes at Q = 0 with divergent slope for −2 < p ≤ −1 cases, takes a finite value there (= 8π/α 1 ) if p = −2 and exhibits a vertical asymptote if p < −2 (see figure 10).
In the A2 cases the integration of Eq.(9), together with (45)-(46) lead, for σ < 2, to  54)). The asymptotic shapes of the A1 and UVD curves depend on the degree of divergence of the field r = 0 (see Eqs. (55) and (56)). The form of the dashed pieces of constant r h curves, matching the small and large-Q regimes, depends on the particular model but must be always monotonically increasing. In fact, the thermodynamically meaningful region, where the external horizons are present, lies under the r h = 0 curve for Q < Qc and the line of EBH solutions for Q > Qc in the A2 cases. For the A1 and UVD cases this region lies under the line of EBHs, which is calculated in the text for the different cases. It is qualitatively represented in the main frame for Q Qc in a particular A2 model with σ > 1 (dashed-dotted line). The small-Q, Q Qc and large-Q behaviours of these EBH curves defining those upper boundaries are plotted in Fig.10) for all possible cases. where whereas for σ > 2 we obtain The derivative of Φ e (Q) can be directly obtained from Eq.(53) and (40) in the general case and reads This equation allows the calculation of the slope of the extreme BH curves around Q c , which is positive and finite in the cases σ < 1, finite positive, null or finite negative if σ = 1 and negative divergent if σ > 1. All these features are displayed in figure 10. As already mentioned, the function Φ e (Q) involves the two thermodynamically conjugate state variables Φ and Q, and defines a vanishing temperature "isotherm" of the set of BH states associated to a given model (see equation (37) and the extremality condition (31)). Thus the analysis of these curves can reveal the presence of van der Waals-like first-order phase transitions. This requires, as a necessary condition, the existence of some values of the charge for which the slope of Φ e (Q) vanishes. Thus, the search for zeroes of (64) is a first step in the identification of such transitions. A glance to figure 10 shows that these special values of the charge must be necessarily present for models belonging to families with A1 or UVD central behaviour (with −2 < p < 0) and B3 asymptotic behaviour. The same arises for models belonging to the UVD (with p < −2) and B1 families, or for those belonging to A2-B1 families with σ > 1. Other interesting cases are the models belonging to the A2-B2 family with σ > 1, when the condition holds. But, beyond this necessary condition, the confirmation of the presence of such phase transitions in a given model requires the simultaneous analysis of the Gibbs free energy, defined as for each set of extreme BHs. Similar methods have been applied to the search of this kind of phase transitions in the particular case of a Born-Infeld model in AdS background [33]. A more complete exploration of these issues is in progress.

The temperature T (r h , Q)
In analyzing the behaviour of the temperature (37) as a function of r h at constant Q let us introduce the function or, equivalently, where the last hand side of Eq.(16) has been used. The derivative of this function is easily obtained and, as a consequence of the positivity of ϕ, satisfies Thus η, at fixed Q, is a monotonically increasing function of r h for any admissible G-NED and, owing to the behaviour of ϕ and E as r h → ∞, approaches asymptotically the constant 1/4π. The behaviour of this function is plotted in figure 11 for representative examples of cases A1 (Euler-Heisenberg) and UVD (Reissner-Nordström) for a given value of Q. The parts of these curves in the negative η region correspond to the inner horizons and must be discarded for the thermodynamical analysis. The straight lines drawn from the origin to the points of these curves in the positive η region have slopes tan(θ) = T , which define the corresponding temperatures. Since these curves are continuous, monotonically increasing and approach asymptotically the constant value η = 1/4π for any Q, there is always an absolute maximum value of the angle θ. Thus, the temperature ranges from zero, on the cut points with the r h axis (corresponding to the extreme BH configurations), to a maximum value T max , corresponding to an angle θ max in Fig.11. This maximum (and, eventually, other local extrema of the curve) corresponds to BH configurations whose horizon radii are solutions of the equation which results from Eqs.(67)-(69). Moreover, from Eq.(67) and the asymptotic behaviour of the energy density in the B-cases it is easily seen that the temperature vanishes as for any large mass BH (this behaviour is related to the asymptotically Schwarzschild character of these space-times). The BH temperature as a function of r h at fixed Q is plotted in the small frame of figure 11. The slope of this curve is given by the derivative defined in Eq.(70). The temperature as a function of M has a similar shape, but the slope of the T − M curve, obtained from Eqs. (67) and (70), reads and diverges for the extreme BH configurations. Owing to the convex character of the curves η(r h , Q) in the particular cases of figure 11 there are at most two cut points with the radial lines for T < T max or, equivalently, there are at most two roots of Eq.(70) in these cases. However, in more general cases, the sign of the curvature of η(r h , Q) at constant Q can change many times as r h increases, leading to many cut points with the radial lines or, equivalently, to many solutions of equation (70). Let us illustrate this behaviour with a particular example. Consider a family of NEDs whose ESS solutions are of the form where R = r/ √ Q, the parameter Q being the charge associated to a given solution and α > 0, β > 0 and the integer n ≥ 2 being fixed constants for a given model of the family. These fields are monotonically decreasing functions of r if the parameters fulfill the (sufficient but not necessary) condition and correspond to the ESS solutions of charge Q associated to a family of admissible NEDs whose Lagrangian densities are parameterized by the constants α, β and n. The form of these Lagrangian densities can be directly obtained from Eqs. (6) and (73)  between extrema support stable phases. We conclude that for models whose ESS black hole solutions belong to cases UVD and A1 (with any Btype asymptotic behaviour) the temperature increases always from zero, for the extreme BH configurations, and vanish asymptotically for large mass BH configurations, as for the RN black holes. In the range of intermediate masses the temperature has always an absolute maximum and can exhibit a set of local extrema, determining a more or less involved phase structure, which is fully characterized by the solutions of Eq.(70). The behaviour of the temperature as a function of the charge at constant r h can be analyzed from its scaling law (R h = r h / √ Q) which can be directly deduced from Eq.(68) and the first integral (6). From this equation we obtain the derivative which is negative for any Q > 0 and vanishes at Q = 0, where the temperature reaches the Schwarzschild value T = (4πr h ) −1 . Figure 13 displays the temperature as a function of the charge at different r h values for the case of field solutions E(r) with monotonic slope (Born-Infeld here, as representative of A2-B2 cases).

The cases A2 and the black points
Although the validity of the analytic discussion of the precedent subsection is general, it has been applied there mostly to the thermal behaviour of BH solutions of UVD and A1 families only. Let us discuss here the cases A2, which lead to more involved (and richer) thermodynamic structures. The method used in determining the behaviour of the temperature as a function of the horizon radius (or of the mass) at constant Q still holds, but in this case the functions η(r h , Q) exhibit new characteristics which lead to configurations with qualitatively new thermodynamic properties. The main difference with the previous cases concerns the behaviour of these functions as r h → 0, which defines the black point configurations. To address this question let us go to the next order in the polynomial expansion (9) of the field solutions of a given A2 model around the center, which can be written as Here λ > σ and Q λ/2 c(Q) = c(Q = 1) = c 0 is a universal constant of the model. From Eq.(68) we see that the behaviour of η(r h ∼ 0, Q) can be obtained from Eq.(77) and the behaviour of ϕ(r h , Q) ≡ ϕ(E 2 (r h , Q)) near the black point configurations (r h ∼ 0). Using (77) and the first integral (6) we are led to if σ, λ = 2. When λ = 2, ϕ behaves as and if σ = 2 In these equations the constants ∆ are independent of the particular solution.
If σ > 2 they are given by ∆ = ϕ(a 2 ) > 0. Otherwise they can be explicitly calculated from the form of the Lagrangians, as functions of the exponents σ < 2, λ and the associated universal coefficients in the polynomial expansion of the field around the center (this may require to go beyond the order of the expansion in Eq.(77)). The corresponding expression for the A2 cases with σ < 2 and λ > 2 is while for σ = 2 we have The expression of ∆ becomes more involved when λ and some of the next-order exponents in the polynomial expansion of E(r ∼ 0) are smaller than two, but its calculation remains straightforward. Equations (77)-(82) replaced in (68) give the form of η for the small radius black holes in the different A2 cases. For the black points (r h = 0), η takes always the finite value which vanishes for extreme BPs (Q = Q c ) and is positive for non-extreme BPs (Q < Q c ). If Q > Q c , the small values of r h (< r he ) correspond to inner horizons, for which η(0, Q) < 0 does not define a meaningful temperature (see Fig.14). Obviously, for non-extreme BPs the temperature diverges in all cases, but the temperature of the extreme BPs is now dependent on the value of the parameter σ of each A2 model. Indeed, from the expansions of η at small r h in the different cases we are led to the several possible expressions for the temperature of the small-radius BHs: If σ = 2, we have and if λ = 2, In all these expressions the first term dominates if Q = Q c , in agreement with Eq.(83). For the extreme BPs the temperature is given by the slope of η(r h , Q c ) in r h = 0 and is null if σ > 1, finite if σ = 1 or divergent if σ < 1 (see Fig.14).
It is now evident that the temperature behaviour in case A2a (η(0, Q > Q c ) < 0) is similar to the ones of the cases UVD and A1 (see the small frame of Fig.14), exhibiting an extreme BH configuration with vanishing temperature Fig. 14 Small and large r h behaviour of η(r h , Q) at constant Q for the A2 cases. The three beams formed by continuous pieces of curves at small r h correspond to the three subcases defined by the values of the charge Q Qc. For each beam the three curves correspond to different values of the parameter σ 1, exhibiting different slopes (vanishing, finite, or divergent, respectively). The asymptotic behaviour of these curves is obtained from Eqs.(68), (8) and (11) and is also displayed with pieces of continuous lines at large r h . As in the previous cases (A1 and UVD), the form of the dashed lines (intermediate region) depends on the particular model and can exhibit more involved shapes than the ones displayed in this plot. The small plot exhibits the form of the temperature for the BI model (for which σ = 4) for several values of Q. Note that the temperature of the EBPs (cases A2c, with Q = Qc) is given by the slope of the corresponding curve at r h = 0 and vanishes in models with σ > 1, takes a finite value in models with σ = 1 and diverges in models with σ < 1. The temperature of the non-extreme BPs (A2b cases, for which Q < Qc) is always divergent.
at finite r h , vanishing temperature BHs for large r h , an absolute maximum temperature BH and (depending on the particular model) many possible BH configurations corresponding to local extrema of T for some values of r h (or M ), which define stable phases in some intermediate ranges. This behaviour is qualitatively independent of the values of the parameter σ.
In case A2b (η(0, Q < Q c ) > 0) there is a continuous spectrum of singlehorizon BHs, whose radii increase from r h = 0 (when M BP = ε(Q) = Q 3/2 ε(Q = 1), corresponding to a non-extreme black point; see figure 2) to r h → ∞ as M → ∞. In these cases the temperature diverges for the black points and vanishes asymptotically for large mass BHs. In the intermediate region, depending on the structure of the ESS fields, the temperature can exhibit local extrema defining ranges of stability, as in the example (73). Again, this behaviour is qualitatively independent on the values of the parameter σ.
In case A2c we have also a BH spectrum going from r h = 0 (corresponding now to an extreme black point) to r h → ∞. The mass of the extreme BPs equals the soliton energy and can be written as where the last term is obtained from Eqs. (21), (35) and (38), and is a universal relation for A2 models. The thermal behaviour of the extreme BPs for the different models is characterized by the values of the parameters σ and λ and can be directly obtained from Eqs.(84)-(86) and their r h derivatives: -For models with σ < 1 the temperature diverges for the extreme BPs, exhibiting a vertical asymptote at r h = 0. -For models with σ = 1 the temperature of the extreme BP takes the finite value (see Eq. (84)) which is a universal constant for each A2 model. The slope of the temperature curve at vanishing r h diverges for λ ≤ 2, being ∓∞ for c 0 ≷ 0, respectively. For λ > 2 this slope is 2∆, which is finite. Here the sign of ∆ determines the stability of the BPs. -For models with 1 < σ ≤ 2 the temperature of the BPs vanishes and its slope at r h = 0 is always positively divergent. -For models with σ > 2 the black point temperatures vanish and the slope of the temperature curves at the origin is 2∆ > 0.
A similar analysis of the second derivatives of the temperature with respect to r h , obtained from equations (84)-(86), give us the concave or convex character of the T − r h curves for configurations around the BPs. A summary of these results for the extreme BP solutions in this A2c cases is given in table 2.

The specific heats C Q (r h , Q)
The specific heat at constant charge of the BH configurations is defined as and we see that, as a function of r h or M , diverges on the relative extrema of the temperature function, exhibiting there vertical asymptote. From Eqs. (67) and (68) we obtain the general formula and eliminating the Lagrangian function ϕ in terms of the temperature and the field we are led to We see from this last expression and the asymptotic behaviour of T , given by Eq.(71), that the specific heat of large mass BH solutions of any admissible B-type model behaves as as should be expected from their asymptotically Schwarzschild character. The slope of C Q as a function of r h reads Using the asymptotic behaviours of the different functions in this formula for the B-cases we see that for large-r h BHs this derivative diverges as −4πr h , consistently with the asymptotic behaviour of C Q in Eq.(92). For the extreme BHs associated to the cases UVD, A1 and A2a we have η = T = 0 and ∂T ∂r h Q = 4ϕ > 0 (see Eq. (69)). Consequently C Q;extr = 0 for any extreme BH configuration of admissible models. The slope of the specific heat as a function of r h , given by Eq.(93), converges to the value for the extreme BHs of any admissible model. The slope of the specific heat as a function of the mass is obtained from Eq.(93) as and diverges in the extreme limit as 1/T .
Let us analyze now the specific heats of the BP configurations discussed above. Using Eq.(91) and replacing there the expressions (84) and (77) for the expansion of the temperature and the field around r h ∼ 0 we are led to and which is valid for σ, λ = 2. It is straightforward to obtain two similar expressions for the cases σ = 2 and λ = 2 from Eqs. (85) and (86). The analysis of these equations leads immediately to the following conclusions: -The specific heats vanish for any non-extreme BP (A2b cases, Q < Q c ), behaving as for small horizon radius BHs. -For extreme BPs (A2c cases, Q = Q c ) in models with σ > 2 the specific heats vanish, behaving as for small horizon radius BHs. -For extreme BPs in models with 1 < σ ≤ 2 the specific heats vanish, behaving as for small horizon radius BHs. -For extreme BPs in models with σ < 1 the specific heats vanish, behaving as for small horizon radius BHs.
-For extreme BPs in models with σ = 1 the dominant terms in the expression of the specific heat of small radius BH configurations reduce to when λ = 2 and when λ = 2. As easily seen, the specific heats of extreme BPs in these σ = 1 cases vanish, excepted for models with λ ≥ 3 and ∆ = 0, as we shall see at once. If λ = 2, C Q behaves for small horizon radius BHs as in Eq.(104), where the sign near zero is given by that of the coefficient c 0 . This behaviour is not dependent on the value of ∆. For λ > 2 and ∆ = 0, the specific heat of small-r h BHs behaves as where ∆ determines now the sign of the approach to zero. For 1 < λ < 2 the specific heat behaves as which vanishes with the opposite sign of c 0 and is not dependent on the value of ∆. Finally, for models with ∆ = 0 the expression (106) of C Q remains valid for any value of λ = 2. In these cases we can have extreme BP solutions with finite specific heat if λ = 3, besides models with divergent specific heat extreme BPs if λ > 3.
In the last line of table 2 we have included these results on the possible behaviours of the specific heats for extreme black point solutions associated to the family of A2-type gravitating NEDs. This discussion exhausts the possible thermodynamic structures of the G-ESS black hole solutions of admissible NEDs minimally coupled to gravity.

Some examples
In order to illustrate different points of this analysis let us consider some particular examples belonging to the different families of admissible models.
The thermodynamic properties of the Reissner-Nordström solution of the Einstein-Maxwell field equations are representative of the typical behaviour of the small-radius BH solutions of the UVD family. Although the gravitating counterpart of this NED has been extensively analyzed in the literature (see [14] and references therein) let us mention a few of its features, owing to this representative character. We have already make use of this example in the small frames of figures 3 and 11 representing the generic behaviour of the mass and temperature as functions of the BH radius for this kind of solutions. We have also displayed the temperature and the specific heat as functions of the mass and of the horizon radius in figure 15. The explicit expressions for these variables in the RN case are straightforwardly obtained from the formulae of the precedent analysis and lead to the expected thermodynamical behaviour for the small and large horizon radius BHs of the UVD-B2 cases. As a representative example of the A1 family let us consider the Euler-Heisenberg effective model of quantum electrodynamics [6,9], whose gravitating counterpart has been analyzed in Ref. [34]. As discussed in the introduction, this is an effective model accounting (to leading order) for screening effects of the Dirac vacuum on electromagnetic wave propagation. The interest on this model in the present context is related to the possibility of explore modifications induced by quantum corrections of the electrostatic field on the gravitational field generated by point-like charges (which, in absence of these corrections is described by the Reissner-Nordström BH solutions). The Lagrangian in this case reads [6] ϕ(X, Y ) = µ being a positive constant. The ESS field solution E(R) is the unique positive root of the equation (R = r/ √ Q). As easily seen, this field behaves at the center and asymptotically as and the model belongs to the A1-B2 family, whose ESS solutions in flat space are of finite energy [8]. The temperature of the BH solutions of this model, as a function of r h , can be obtained from Eq.(68) in terms of the root of (109) as and the expression for the specific heat results directly from Eq.(91). The mass as a function of r h for several values of Q has been plotted in figure 3 for this EH model (main frame), where also the behaviour of the RN model is displayed (small frame). The only qualitative difference between both models arises on the small-r h region, where the EH model leads to one-horizon BH solutions for masses over the total electrostatic energy (M > ε(Q)). In the thermodynamically pertinent regions in these figures, which lie between the (dotted) lines of extreme BHs and the (dashed) line of Schwarzschild BHs, the behaviour of the BH masses versus r h are qualitatively the same. Moreover, the behaviours of the functions η(r h ) and T (r h ) at fixed Q for these models (plotted in figure 11) as well as the behaviour of the specific heats (plotted in figure 15 for the RN model) are also similar. These results illustrate the qualitative similarity of the thermodynamic properties of the BH solutions of UVD and A1 families near the extreme BH configurations.
A first representative example of the A2 family is the well known Born-Infeld model [4]. As mentioned in the introduction, this model is a nonlinear generalization of classical electrodynamics which avoids the problem of the divergent self-energy of point-like charge solutions. Moreover, BI-like models arise as effective gauge field theories in the low-energy regime of string theory [3]. Besides their use in the study of black hole physics [20,21] they have found applications in cosmological contexts [5]. The corresponding Lagrangian density reads where 1/µ is the maximum field strength. When µ → 0 in Eq.(112), the Maxwell Lagrangian (ϕ = X/2) is recovered. The ESS solutions corresponding to model (112) take the form The expansion of these solutions around r = 0 gives the values of the maximum field strength a = 1/µ, the coefficients b 0 = a 3 /8, c 0 = 3a 5 /128 and the exponents σ = 4 and λ = 8 (see Eq. (77)). The constant ∆ in the expansion (78) of the Lagrangian density is now ∆ = 1/µ 2 = a 2 . The expression of the BH temperature as a function of r h and Q is obtained from the equations (68), (112) and (113) and reads which behaves as for small radius BHs. So we recover the results of the precedent analysis on the temperature of the BH configurations close to the extreme (Q = Q c ) and nonextreme (Q < Q c ) BPs, described by equation (84). For values of the charge Q > Q c (A2a cases) the radii of the extreme BH configurations as functions of Q, obtained from Eq.(31), take the form and, as expected, the temperature (115) vanishes for these configurations. Similar considerations can be done for the behaviour of the specific heat, which results directly from Eqs.(91), (113) and (115). For small radius BHs the behaviour of C Q is obtained from these formulae or, alternatively, from Eq.(99) in case A2b or from Eq.(100) in case A2c, once the characteristic parameters of the model have been fixed. Figure 16 shows the behaviour of the temperature, for these BI black holes, as a function of r h for several values of the charge, corresponding to the cases A2b (two upper curves), A2a (lower curve) and A2c. We have also plotted the behaviour of the specific heats for the two A2b cases and the A2c case.
The A2 models with σ = 1 lead to new thermodynamic behaviours. In particular, the temperature of the extreme black points becomes now finite. Let us give two examples of these cases, which are constructed by convenience, since they are illustrative of the different statements on the state variables discussed in section 4. The first one corresponds to a family of NEDs defined through the expression of their ESS solutions, which take the form and the other, obtained for lower values of Q, behaving nearly as the Schwarzschild BH temperature. The limit r h = 0 corresponds to the divergent-temperature non-extreme BPs. The curve starting at the origin corresponds to the critical value of the charge and exhibits a EBP there (case A2c). The lower curve, corresponding to Q > Qc (A2a cases) exhibits a behaviour similar to that of the UVD and A1 cases. The associated specific heats are plotted in the other panels: upper right and lower left panels correspond to the two A2b temperature curves. The lower right panel corresponds to the cases A2c, with the curve starting at r h = 0, but it is also representative of the cases A2a if the curve starts on a finite value of r h , corresponding to the EBH configurations. The vertical asymptote of C Q correspond to local extrema of the temperature.
where the exponents λ > 1 and µ > λ + 1, and the positive constants a, b 0 and c 0 are to be chosen in such a way that E(R) be positive and monotonically decreasing everywhere, as required for admissibility. From Eq.(118) and the first integral (6) we can obtain, through a quadrature, the expression of the admissible Lagrangian densities supporting these soliton-like elementary solutions. We shall make the choice a = 1, b 0 = 0.2, c 0 = 0.1, µ = 9/2 and two different choices for the exponent λ (λ = 3/2, corresponding to a B3 asymptotic behaviour or λ = 5/2, corresponding to a B2 asymptotic behaviour). Both choices lead to admissible models, for which the coefficients of the expansion of E(R) around R = 0 coincide with those of the numerator of (118). The two choices for the exponent λ lead to the same value of the extreme BP temperature but to different behaviours for the temperature of small-radius BHs (see table 2). For λ = 3/2 < 2 the temperature decreases from its BP value (= 2b 0 / √ πa) with infinite negative slope whereas for λ = 5/2 > 2 the temperature increases from the same value at r h = 0 with finite slope (= 2∆). In both cases the specific heats of the extreme BP configurations vanish, as for the model with λ = 3/2, and as for the model with λ = 5/2. Figure 17 shows the behaviour of the temperature as a function of r h at fixed Q = Q c for both models. Even though this behaviour is similar for finite horizon radius, the qualitative difference around r h = 0 (extreme BP configurations) is apparent. Finally, let us consider a family of models representative of the subclass of the A2 NEDs satisfying the condition σ = 1, λ = 3 and ∆ = 0. The assumed form of the Lagrangian for X ≥ 0 is where α is a positive constant parameterizing the family. It is easy to verify that, as a function of X, at Y = 0, this expression increases monotonically from ϕ(0, 0) = 0 and exhibits a vertical asymptote at X = α. Its slope at the origin takes the finite value and, consequently, this is a one-parameter family of admissible models belonging to the class A2-B2 (obviously, for X < 0 these Lagrangians should be properly extended taking into account the admissibility conditions). The general ESS solution of these models can be obtained explicitly in terms of the parameter α and the charge Q by using the first integral (6) and reduces to solve a quartic algebraic equation. Nevertheless, we are interested here in the behaviour of the field and the thermodynamic functions near the extreme BP configurations. By replacing the expression (77) and the derivative of (121) in the first integral we determine straightforwardly the parameters of the first terms of the polynomial expansion of the field around the center as or, equivalently, and we see that α must be interpreted as the maximum field strength. Moreover, using Eq.(81) we are led to the value ∆ = 0 for the integration constant of the expansion of the Lagrangian around r h = 0. Consequently, a glance to table 2 shows that the temperature of the extreme BP configuration in these cases is T EBP = 2/ √ πα, the slope of the temperature at r h = 0 vanishes and becomes negative as r h increases, and the specific heats attain, for the extreme BP configurations, the negative value C Q;EBP = −α 2 /40. Figure 18 displays these temperature and specific heat behaviours for the critical (Q = Q c ) configurations, which become atypical in the limit of extreme BP configurations.  Fig. 18 Behaviour of T (r h , Q C ) and C Q (r h , Q C ) for the model (121). As in the case of figure 17 the temperature of the EBP configuration is finite, but now its slope vanishes there.

Finite thermodynamic relations, scaling group and generating equations
Let us consider now some finite relations arising between the state functions introduced in the precedent section.

Generalized Smarr law
In recalculating the derivative of ε ex with respect to the charge at constant r h in equation (34) we can use the scaling law of the external energy function (20) and we are led to the relation or, alternatively, where the dependence in r h has been replaced by dependence in the entropy S. In the particular case of the RN black hole this expression boils down to which is the well known Smarr formula [35]. As another example, for the NED family defined in Ref. [18] by the Lagrangian density ϕ(X, Y = 0) = X p (particularized here to three space dimensions), p being an integer number, Eq.(126) provides the result which is precisely the Smarr-like law obtained in that reference. Moreover, it is straightforward to show that the equivalence of (126) and (127), which requires arises only for the ESS black hole solutions of the Einstein-Maxwell model. This condition is related to the correspondence between the energy content of the field created by a system of charges at rest and their electrostatic potential energy, which holds in flat-space Maxwell electrodynamics [36], but not for general NEDs. Thus, for any admissible G-NED the function V (r h , Q) gives the deviations from the exact Smarr law, induced by the effect of the nonlinear self-couplings of the electric field. Owing to the scaling properties (52) and (20) this function scales as the external energy (R h = r h / √ Q) and is fully determined by its characteristic expression for Q = 1.
We conclude that Eq.(126) is the natural generalization for admissible G-NEDs of the Smarr formula (127) of the Einstein-Maxwell theory.

Other finite relations
As already mentioned, the thermodynamical approach regards the ESS black hole solutions obtained from a given gravitating NED as different states of a unique system characterized by two state variables, the most immediate ones being the integration constants M and Q, which arise in the resolution of the Einstein equations (24). On the other hand, in solving the thermodynamic problem for a given model, we start with the first integral of the electromagnetic field equations (6) and we can obtain the state variables as functions of r h and Q. New state variables can be defined as α = α(r h , Q) and β = β(r h , Q), where these functions are assumed to be invertible in some pertinent ranges, leading to where the functions f and g can be multi-valuated. The universal scaling laws obtained in section 4 allow to build new finite relations between different thermodynamic variables, whose form is now dependent on the particular model. Let us write these scaling laws in a new (but equivalent) form for the main state variables introduced so far. First of all, the scaling law (7) for the electrostatic field leads to where θ is any positive parameter (which can be a state variable). From Eq.(75) we obtain for the temperature From Eq.(48) we obtain for the mass whereas for the electrostatic potential on the horizon we obtain, from Eq.(52), Let us denote as Γ (θ) (θ > 0) these scaling transformations. From Eqs.(132)-(135) it is easy to show that the product of two transformations satisfies Consequently the set of scaling transformations for the thermodynamical state functions associated to the BH solutions of any admissible G-NED has the structure of a multiplicative one-parameter group. From these laws and the identification of θ with different state variables we can obtain several finite relations whose generic forms are the same for all models, but involve particular expressions associated to each NED. Let us give two simple examples. Consider first the function Φ(r h , Q) obtained from the first integral (6) (after a quadrature) for a particular model. Using the scaling law (135) with θ = T we are led to where r h has been written in terms of the entropy. In the same way, the explicit expression of the temperature (67) and the scaling law (133), with θ = Φ, lead to These equations, together with the generalized Smarr law (126), allow for the elimination of two of the five state variables Φ, Q, M, T and S, leading to different expressions for the "equation of state" of the set of BHs associated to a given model.

Generating equations
The universality of the scaling laws allows for a more direct and enlightening procedure in generating the equations of state of general BH configurations. These laws lead to finite relations involving the derivatives of the state functions and are also "universal", in the sense that they are first-order differential equations whose form is independent of the particular gravitating NED. The characteristic curves of these equations generate, through properly chosen boundary conditions, surfaces in the space of state variables involved in the particular scaling law. The points of each one of these surfaces correspond to BH states associated to a given gravitating NED. As the structure of the beam of characteristic curves is also universal, this procedure allows for a separate analysis of some model-independent thermodynamic properties and those model-dependent behaviours induced by the particular boundary conditions corresponding to each NED. Let us present here a particular example illustrating the method and some of the expected results. Using the expression (53) of the Q-derivative of the electrostatic potential at constant r h and replacing the electrostatic field on the horizon in terms of the gradient of Φ we are led to This equation (which can also be obtained directly by derivation of Eq.(135) with respect to the parameter θ at θ = 1) involves the function Φ(r h , Q) and its first derivatives with respect to these arguments. It is remarkable because its "universality", since it must be satisfied by the solutions of any G-NED, as a direct consequence of the scaling law (135), which is a universal property for field theories associated to the action (22). On the other hand, it is independent of the general laws of thermodynamics, as can be easily checked for other thermodynamical systems. Thus, the particular forms of the "equations of state" of the systems of ESS black hole solutions associated to the different particular models, characterized by the two independent parameters r h and Q, are solutions of this equation, corresponding to surfaces in the Φ − r h − Q space. In order to explore further this interpretation, let us write Eq.(139) under the form where the horizon radius has been replaced by the entropy as independent variable. This expression suggests the introduction of the new variables in terms of which the equation takes the form This is a very simple first-order, linear, partial differential equation whose set of characteristics is the beam of straight lines parallel to the vector (1, 1, 1) in the (x, y, z) space. Consequently, the equations of state of the system of ESS black hole solutions associated to the G-NEDs are represented in this space by ruled surfaces generated by lines of this beam. In order to choose the particular surface associated to the system of BH solutions of a particular G-NED we can solve, for example, the first integral (6) for the characteristic (Q = 1) case and determine the form of Φ(r, Q = 1) through a quadrature. This defines a curve z = z(x = 0, y) in the y − z plane, which is a boundary condition for Eq.(142) (see figure 19) determining a unique surface solution z = z(x, y) in the (x, y, z) space, generated by the beam of characteristics. In the (Φ, Q, S) space the parametric equations for the beam of characteristics of (140) read where τ is the (positive) parameter and Q 0 , S 0 and Φ 0 are the coordinates of the points defining the different curves of the beam, corresponding to a particular BH configuration. These curves are parabolae having the origin as a common vertex, where the beam is singular. A similar analysis can be carried out from the scaling laws of other state variables. For example, the scaling law for the temperature (133) leads to the generating equation where x and y are the same independent variables defined in Eq.(141). This is a linear, inhomogeneous, first-order, partial differential equation, generating the state surfaces T = T (x, y) for the different G-NEDs, once the corresponding boundary conditions are specified. In the same way, the scaling law for the BH mass (134) leads to the generating equation whose set of solutions contains the state surfaces M = M (x, y) associated with all the admissible G-NEDs. These preliminary results require a more detailed investigation of the structures which lie at the ground of these scaling symmetries and their consequences for the BH thermodynamics. A systematic exploration of this issue is in progress, but goes beyond the scope of this paper.

Conclusions and perspectives
In this paper we have performed a generalization of the well known thermodynamic analysis of the RN black hole solutions of the Einstein Maxwell field equations to the BH solutions of any admissible gravitating NED. The exhaustive classification of the admissible NEDs in terms of the central and asymptotic behaviours of their elementary solutions in flat space, and the analysis of the geometric structure of their gravitating counterparts, have been carried out in references [10] and [11]. There, the asymptotic flatness of any gravitating ESS solution of these admissible NEDs has been established, regardless of the asymptotic behaviour of their integral of energy in flat space. When this integral diverges at large r (always slowly than r, owing to the admissibility conditions) we are led to the families of asymptotically IRD models, whose gravitating BH solutions exhibit geometric structures similar to those of the other models at finite radius, but approach flatness at large r slower than the Schwarzschild field. For such asymptotically anomalous models the standard thermodynamic analysis is meaningless, at least in the usual way (undefined ADM masses and asymptotic divergence of the electrostatic potentials).
For models whose ESS solutions have asymptotically convergent integrals of energy in flat space (B models) the gravitating elementary solutions behave as the Schwarzschild field at large r, no matter the asymptotic degree of convergence of these integrals. In such cases the ADM mass of the BH solutions is well defined, the electrostatic potential can be set to zero asymptotically and the thermodynamic laws can be consistently established. The thermodynamic properties of the BH solutions in these cases have been analyzed in terms of the behaviours of the integral of energy and the field around the center of the flat space ESS solutions. When this integral (and, hence, the field) diverge as r → 0 (UVD models) the geometric structures (two-horizons BHs, single-horizon extreme BHs or naked singularities) and the associated thermodynamical behaviours are qualitatively similar to those of the Reissner Nordström BHs, although the temperature and the specific heats as functions of the horizon radius (or of the mass) can exhibit more involved behaviours, related to the more or less involved structure of E(r) in each case. For models exhibiting finite-energy ESS solutions in flat space and whose electrostatic field diverges at the center (slower than r −1 , A1 models) the associated BHs can exhibit new geometric features, as compared with the UVD cases (mainly the existence of non-extreme single-horizon BHs), but the form of the mass-radius relations and the behaviours outside the event horizon are qualitatively similar in both cases and hence we are led to similar thermodynamic properties.
For models whose ESS solutions in flat space are finite-energy with finite maximum field strength (A2 family, to which the Born-Infeld model belongs) the set of associated gravitating ESS solutions exhibits a richer structure, including extreme and non-extreme black points as new possible configurations. These new configurations have thermodynamic properties which are related to the behaviour of the electrostatic field around the center. This relation has been extensively analyzed and, in particular, we have determined the conditions to be satisfied by these A2 models in order to support extreme black point solutions with zero, finite or divergent temperature. Concerning the nonextreme black points, their temperature is always divergent.
We have performed a detailed analysis of the properties of the extremal BH solutions for the different families. The thermodynamic behaviour of these configurations differ strongly among the different families and suggest the possible existence of models exhibiting van der Waals-like first-order phase transitions in the Φ − Q plane.
Once the possible structures and thermodynamic properties of the BH solutions of the different families of admissible NEDs have been determined and classified, the results of the present analysis allow for the explicit construction of specific models whose solutions exhibit prescribed behaviours. In this way we have built some new models as examples illustrating different aspects of the study carried out here.
From the invariance of scale of the general NED field equations in flat space, we have established universal scaling laws, satisfied by the thermodynamic functions defining the BH states of any model. These scaling laws lead to universal first-order partial differential equations involving sets of three state variables, whose solutions are surfaces in the three-dimensional spaces defined by these variables. Each one of these surfaces is generated by the beam of characteristics associated to the equation. From appropriate boundary conditions we can select the surface solution associated to a particular G-NED. The points of this surface correspond to the BH states which are the solutions of this model, and its analytic equation, relating the involved variables, defines a expression of the "equation of state" of this system of BHs.
To conclude, let us mention some perspectives and possible extensions of this study: First, in Ref. [8] we have established that the analysis of the ESS problem for NEDs in flat space can be extended to the same problem for generalized non-Abelian gauge field theories of compact semi-simple Lie groups. More explicitly, we have shown that the ESS solutions of a generalized gauge field theory, defined by a Lagrangian density of the form (1) (in terms of the usual gauge invariants with the ordinary trace definition in this non-Abelian case) have the same functional dependence in r as those of the NED defined by a Lagrangian density of the same form (in terms of the Abelian gauge invariants). In this way we have recently communicated some preliminary results on the extension of the analysis performed here to systems of non-rotating BHs of non-Abelian origin [37]. A more detailed report on this issue is in progress.
Second, the extensions of the present general analysis (geometric structures and thermodynamics) to the cases of NEDs coupled to more general actions for the gravitational field (e.g. with a cosmological constant term, or Gauss-Bonnet and Lovelock theories) is in progress, a context where some particular cases of NEDs have been already explored in the literature [26].
Finally, let us stress that a more complete analysis of the scaling symmetries and their consequences on the BH thermodynamics is necessary. In particular, the study of the structures of the beams of characteristics associated to different generating equations and of the boundary conditions selecting the equations of state of the BH systems associated with different NEDs, allows for the exploration of eventual phase transitions, an issue also under development.