Nonlinear Yang-Mills black holes

This paper is devoted to investigating the nonlinear non-abelian Yang-Mills black holes. We consider three Born-Infeld, exponential, and logarithmic nonlinear Yang-Mills theories with $SO(n-1)$ and $SO(n-2,1)$ semi-simple groups, which n is the dimension of spacetime, and obtain a new class of nonlinear Yang-Mills (NYM) black hole solutions. Depending on the values of dimension $n$, Yang-Mills charge $e$ and the mass $m$ and nonlinear parameters $\beta$, our solutions can lead to a naked singularity, a black hole with two horizons, an extreme or a Schwarzschild-type black hole. We also investigate the thermodynamic behaviors of the NYM black holes. For small charge values, the NYM solutions may be thermally stable in the canonical ensemble, if we consider an AdS spacetime with spherical $k=+1$ and hyperbolic $k=-1$ coordinates or a flat one with $k=+1$. However, there are no stable regions in the grand canonical ensemble in higher dimensions. For the NYM black hole, we observe a reentrant phase transition between large and small black holes in the BI-branch with small $\beta$, which cannot be visible for the nonlinear Reissner-Nordstrom AdS black hole in the higher dimension. For the limit $\beta\rightarrow\infty$, the critical ratio $\frac{P_{c} v_{c}}{T_{c}}$ tends to the constant value $3/8$ for each dimension $n$, while it depends on the dimension for the case of nonlinear electrodynamics black holes.


II. MAIN STRUCTURE OF THE NYM BLACK HOLE SOLUTIONS
The theory of the n-dimensional (n ≥ 4) nonlinear-Yang-Mills (NYM) black hole solutions originate from the action where R and Λ are respectively the Ricci scalar and the cosmological constant. We classify the nonlinear Lagrangian L(F ) for three Born-Infeld, exponential and logarithmic cases [54] L(F ) = where we have abbreviated the exponential and logarithmic nonlinear cases to EN and LN, respectively. If we consider a gauge group with N-parameters, so where Γ ab = C c ad C d bc is the metric tensor of the gauge group and detΓ ab is the related determinant. The indices a, b, c take values from 1 to N and C a bc 's are the structure constants of the gauge group theory. The gauge field tensor is defined as where e is the coupling constant of the non-abelian theory and A (a) µ 's represent the gauge potentials. For simplicity, we use the redefinition L(F ) = ζ 1 β 2 L(Y ), where and ζ 1 = +4, +4, −8 and Y = F 2 2β 2 , F 2 4β 2 , F 2 8β 2 are described for BI, EN and LN Yang-Mills theories, respectively. If we vary the action (1) with respect to the metric g µν and the gauge potential A (a) µ , then the gravitational and gauge field equations are obtained as follows where ζ 2 = −4, −2, +2 and ζ 3 = +4, +4, −8 have been specified for BI, EN and LN theories, respectively. To obtain a set of static and spherically symmetric non-abelian black hole solutions with spherical and hyperbolic horizons, we consider the following metric where dΩ 2 k denotes the line element of an (n−2)-dimensional hypersurface Σ with the constant curvature (n−2)(n−3)k. It is determined as where k = 1, −1 are devoted to the spherical and hyperbolic geometries, respectively, and θ ∈ [0, π 2 ]. If we introduce the coordinates x i 's sin(φ j ) , i = 2, ..., n − 2 and employ the Wu-Yang ansatz, then the gauge potentials are obtained from A (a) = e r 2 (x a dx n−1 − x n−1 dx a ) for a = 1, ..., n − 2 .., n − 3 , j = 2, ..., n − 2, and i < j (11) where b goes from (n − 1) to (n − 1)(n − 2)/2. The Lie algebra of the gauge potentials with k = +1 and −1 in Eq. (11) is isomorphic to SO(n − 1) and SO(n − 2, 1) gauge groups, respectively. It should be noted that n is equal to the spacetime dimension. We redefine γ ab in Eq. (3) as γ ab = ǫ a δ ab , no sum on a, where for SO(n − 1) gauge group, ǫ a = 1 , for a = 1 , ..., (n − 1)(n − 2) 2 , and for SO(n − 2, 1) gauge group For a better understanding, we have written the gauge potentials of the groups SO(3), SO(2, 1), SO(4) and SO(3, 1) in appendix (VIII A). The gauge potentials (11) can satisfy the gauge field equation (7). So, if we substitute these gauge potentials (11) and the metric (8) in Eq. (6), we can obtain the NYM black hole solutions where η = (n−2)(n−3)e 2 β 2 r 4 . The parameter m is an integration constant relating to the mass of the NYM black hole. It should be noted that for n = 4z + 1 where z ∈ N , the metric function in Eq. (15) can be written in terms of elementary functions. We have shown the function f (r) for n = 5 and n = 9 dimensions in appendix (VIII B). For n = 4z + 1, f (r) may be written as Eq. (16) shows that there is an equivalence between the four-dimensional NYM black hole solutions with SO(3) and SO(2, 1) gauge groups and a set of topological black hole solutions with k = 1 and k = −1 in nonlinear electrodynamics theory [54,55]. Therefore, we can deduce that there is a transformation between the non-abelian gauge fields and a set of abelian ones in n = 4 which satisfies the Yasskin theory [26]. However, for n > 4, we achieve a new class of solutions for the NYM black hole which is different from the nonlinear electrodynamics one [54,55]. So, the NYM solutions with n > 4 do not respect the Yasskin theorem. For large β, we assume that all three types of the metric functions in Eq. (15) reduce to the Einstein-Yang-Mills black hole solution as follows We choose r 0 = 1 for simplicity.

III. PHYSICAL BEHAVIORS OF THE NYM BLACK HOLE SOLUTIONS
In this section, we aim to study the physical structures of the NYM black hole solutions. If we calculate the Kretschmann scalar, R abcd R abcd , it goes to infinity as r → 0. Therefore, we can deduce an essential singularity located at r = 0 for the NYM black holes.
In the previous section, we concluded that the NYM and nonlinear electrodynamics black hole solutions are the same in n = 4. In Refs. [56][57][58][59], the authors have discussed the structure of the nonlinear electrodynamics black holes horizon in n = 4 and higher dimensions. To find the possible horizons, we investigate the behavior of the metric function f (r) near r = 0. We have shown the horizon structure for the Born-Infeld Yang-Mills (BIYM) black hole in n = 5, 6 and for k = 1 in Fig. (1). Considering k = 1, and Solving Eq. (15) for n = 5 the metric function is where C 5 is the integration constant for n = 5, which is related to the integral in Eq. (15). The integral in Eq. (15) is indefinite, and so we have considered the integration constant. We assume that the expansion of f (r) at large values of β (β → ∞) reduces to the Yang-Mills (YM) solution in Eq. (17) and so we find Now, the expansion of f (r) close to r = 0 in Eq. (18) takes the following form where A 5 is the 'marginal' mass for n = 5, which depends on the values of parameters β and e. As we observe in Fig.  (1(a)), independent of the parameters, f (r) goes to ∞ as r → ∞. However, for r → 0, we have the following cases: For the marginal case which is characterized by m = A 5 , the function f (r) in Eq. (20) has a finite value at r = 0 which is For m > A 5 , the function f (r) goes to −∞ as r → 0. Therefore, there is just one horizon and the AdS-BIYM black hole behavior is analogous to the Schwarzschild black hole (we abbreviate it to Schw-type). For m < A 5 , the solution goes to ∞ at the limit r → 0. Thus, the BIYM black hole has a similar behavior like the 'Reissner-Nordström' black hole (we abbreviate it to RN-type). In this case, the black hole may have zero (naked singularity), one (extremal black hole) or two horizons. Horizons (r + ) are the roots of the equation f (r + ) = 0. If the finite value of the metric function in Eq. (21) is positive,(i.e., for βe < √ 3 2 ), then the solution with m < A 5 leads to a naked singularity and so the only solution is the Swch-type with m > A 5 . However, when the finite value in Eq. (21) is negative (βe > √ 3 2 ), the solution with m < A 5 can describe a black hole with horizons for m ex < m < A 5 . The parameter m ex is the mass of the extremal black hole, which is determined from the conditions f (r = r ex ) = 0 and f ′ (r = r ex ) = 0. For n = 5, it is given by We have probed the horizon structure of the BIYM black hole for n = 6 In Fig. (1(b)). For n = 6, the expansion of the metric function in Eq. (15) around r = 0 becomes This shows that the marginal case in n = 6 happens only for m < 0. Therefore, the six-dimensional BIYM black hole has only one horizon when m > 0, which is the Schw-type. This behavior is a general feature of the nonlinear Yang-Mills black holes in some higher dimensions that the marginal mass is negative and there can be only one horizon and thus the Schw-type is the only solution.
In Figs. (2(a)) and (2(b)), we have investigated the horizon structures of the exponential and logarithmic nonlinear Yang-Mills black holes (we abbreviate them to ENYM and LNYM, respectively) in five dimensions. We observe the same behavior as the BIYM case when r → 0. Expanding the metric function, one can examine the horizon structure near the origin for ENYM and LNYM cases in the same way as the BIYM case. The expansions of the function f (r) in Eq. (15) around r = 0 for the ENYM and LNYM cases in n = 5 and for k = 1 are given by and respectively, where γ is Euler-Mascheroni constant. As we observe in Figs. (2(a)) and (2(b)) when m = A 5 we have f (r) = 1. For m > A 5 the behavior of the metric function is Schw-type. For m < A 5 depending on the parameters β and e, we may have zero(naked singularity), one(extremal black hole) or two horizons. In fact, the marginal mass, which depends on both e and β, is a boundary(or a margin) between two qualitatively different kinds of solutions shown in Fig. (1(a)) and Fig. (2). If the constant of integration m is larger than the marginal mass then the exact solution in Eq. (15) has only one horizon that is similar to the Schwarzchild black hole behavior, despite the fact that the black hole is charged. If m is smaller than the marginal mass then the exact solution (15) has two horizons, which is the same as Reissner Nordström behavior. The marginal mass forms the boundary between these two cases for n = 5 and some higher dimensions. We can also investigate the horizon structure of the solutions in higher dimensions. In general, the expansion of the metric function f (r) around r = 0 in Eq. (15) for the BIYM case is given by and for the ENYM and LNYM cases are and respectively, where C n is the integration constant for dimension n. One may obtain a value for C n by assuming that for large values of β, f (r) tends to the Yang-Mills solution.

IV. THERMODYNAMIC QUANTITIES OF THE NYM BLACK HOLE SOLUTIONS
According to the gauge/gravity duality, one can provide a relation between the strongly coupled gauge theories and the related weakly coupled string theories. From the holography viewpoint, the bulk string theory may inform the boundary gauge theory. Using AdS/CFT correspondence [60,61] which maps the conformal field theory to the asymptotically AdS spacetime with a higher dimension, thermodynamic properties of a black hole may reveal the dual physical system properties. For instance, the horizon of a black hole in the asymptotically AdS spacetime can give information about the finite temperature of its dual field theory. In this part, we would like to investigate the thermodynamic quantities of the NYM black hole and then check the first law of thermodynamics. The Hawking temperature of the NYM black hole may be obtained from where we have used f (r + ) = 0. From the so-called area law (which states that the entropy is one-quarter of the event horizon area of the black hole [62]), the entropy can be calculated as below where ω n−2 represents the volume of a (n − 2)-dimensional unit sphere and a (n − 2)-dimensional hypersurface with constant negative curvature for k = 1 and k = −1 respectively. To obtain the mass, we use the subtraction method of Brown and York [63]. For this purpose, let us write the metric (8) in the following form and choose a background metric with V 0 (r) is an arbitrary function that determines the zero of the energy to avoid the infinities of the mass. We again choose r 0 = 1. If we characterize σ ab as the metric of the spacelike surface Σ in ∂M, and n a and ξ b as the unit normal and the timelike killing vectors of this boundary, respectively, the mass of this black hole is calculated by where σ is the determinant of the metric σ ab and K 0 ab is the extrinsic curvature tensor of the background metric. As we use the limit r → ∞, the mass of the NYM black hole yields to where m is found from the equation f (r + ) = 0 in Eq. (15). The global Yang-Mills charge of the NYM black hole is obtained from the Gauss law We consider the mass M in Eq. (34) as a function of the entropy (30) and the charge (35), so the first law of thermodynamics is specified by where T = ∂M ∂S Q and Φ is the gauge potential, Φ = ∂M ∂Q S . Numerical calculations demonstrate that the evaluated T is in agreement with T + in Eq. (29). So, the first law of the black hole thermodynamics is satisfied, if we obtain the gauge potential of the solutions with n = 4z + 1 as below The gauge potential of the solutions with n = 4z + 1 can be derived exactly for each dimension, however it does not have a general form.

V. THERMAL STABILITY OF THE NYM BLACK HOLE SOLUTIONS
The thermal stability of a black hole is determined if one analyzes the behavior of the entropy S or the energy M with respect to the small variations of the thermodynamic coordinates around the equilibrium. In this section, we aim to consider S and Q as a set of thermodynamic variables and probe the NYM black hole thermal stability in the canonical and grand canonical ensembles. In the canonical ensemble, the parameter charge Q is fixed and so the black hole is thermally stable if the heat capacity is positive. It should be noted that in order to have physical solutions, the temperature should be also positive. So, a physically stable NYM black hole in the canonical ensemble is obtained, if the conditions T + > 0 and C Q > 0 are satisfied. In the grand canonical ensemble, both parameters S and Q are variables and so the positive value of the Hessian matrix determinant may lead to stable solutions. The Hessian matrix is defined as where we abbreviate the related determinant to det(H). In this ensemble, the two conditions ∂ 2 M ∂S 2 Q > 0 and ∂Q 2 S > 0 should be satisfied. If all three quantities, temperature T + , heat capacity C Q and det(H) are positive, then these two conditions are established spontaneously from Eqs. (38) and (39).
To find a physical stable region for the NYM black hole in both canonical and grand canonical ensembles, we have plotted the temperature T + , heat capacity C Q and Hessian-matrix determinant det(H) versus r + in Figs. (3)(4)(5)(6). As the third term in Eq. (29) is negative for the three BIYM, ENYM, and LNYM cases, a positive value for the temperature depends on the values of the first two terms. To reduce the effect of the third term in Eq. (29), we choose a small value for the parameter Q. We can also find from Eq. (29) that the temperature is positive just for the AdS (Λ < 0) solutions with k = ±1 , and also for dS(Λ > 0) and flat (Λ = 0) solutions with k = 1. We have probed the solutions with these features in Figs. (3)(4)(5)(6). We refuse to investigate the thermal stability of the dS solutions since our results have shown that it is not possible to obtain a positive region for both quantities T + and C Q .   5) and (6). In Fig. (3), the temperature is positive for all values of r + in dimensions n = 4, 5, 6. So, the thermal stability of these solutions in the canonical ensemble depends only on the heat capacity value. There is a r +min1 which the heat capacity is positive for r + > r +min1 and it increases as the dimension n increases. To obtain a physically stable region in the grand canonical ensemble, all the quantities T + , C Q and det(H) should be positive. We can conclude from Fig.  (3) that the stable region becomes smaller as the dimension n increases. For example, the four-dimensional BIYM solutions are stable for r + > r +min1 , while there is no stable regions for n = 6. For n = 5, there is just a small stable region for the BIYM black hole. Obviously, the obtained stable regions are different in these two ensembles, which is expected. However, one may choose the cosmological constant as a thermodynamics variable. It is argued in [64] that considering cosmological constant as a variable may lead to the same stable regions for the two ensembles.
For the AdS-ENYM black hole in Fig. (4), there is a r +min2 for each dimensions n = 4, 5, 6, which both T + and C Q are positive for r + > r +min2 . By increasing the dimension n, the positive value of det(H) decreases. The fourdimensional AdS-ENYM solutions with k = −1 are thermally stable for r + > r +min2 in the grand canonical ensemble, however, there is no thermal stability for n = 6.
In Fig. (5), we have probed the stability of the flat LNYM solutions with k = 1. As the heat capacity behavior is not clear in n = 4, 5, 6 for the range of 0 ≤ r + ≤ 1, so we have magnified it in Fig. (6). Figs. (5) and (6) show that there are two values r +min3 and r +max which both T + and C Q are positive for r +min3 ≤ r + ≤ r +max . As the dimension n increases, the values of r +min3 and r +max decrease. As the quantity det(H) is negative for the range r +min3 ≤ r + ≤ r +max , so we cannot find a stable region for the flat solutions in the grand canonical ensemble.

VI. CRITICAL BEHAVIOR OF THE NYM BLACK HOLE SOLUTIONS
In this section, we would like to study the critical behavior and phase transitions of the NYM black holes in the extended phase space. One can enlarge the thermodynamic phase space and consider the cosmological constant as a thermodynamic pressure [65][66][67]. Hawking and Page were the first who showed a phase transition for the Schwarzschild AdS black hole [68]. Recently, many studies about the critical behavior and phase transition of black holes have been done [69][70][71][72]. The critical behavior of the BI Maxwell (BIM) black hole in the AdS spacetime has been also investigated [56,58,73]. In Ref. [74], the critical behavior of the BI-dilaton black holes has been discussed. In the following, we first obtain a Smarr relation, then we get to an equation of state and a Gibbs energy to check out the phase transition of the NYM black hole. We also obtain the critical exponents of this black hole and compare them with the Van der Waals fluid.

A. Smarr relation
To obtain a Smarr relation, we should investigate the thermodynamics of the black hole in the extended phase space. For the NYM black hole, we consider the quantities S and Q, the dimensionful parameters Λ and β, and their conjugates as thermodynamic variables. In this way, we can write the first law of thermodynamics in the extended phase space where the pressure P is defined as If we consider the specific volume v = 4r+ n−2 and use Eq. (34), then the conjugate quantity of P is and the related conjugate of β in the n = 4z + 1 case is defined as It should be noted that the Smarr relation is not satisfied for the case n = 4z + 1. This is not unexpected as there are some other black hole solutions in the context of nonlinear electrodynamics for which the Smarr relation is not satisfied [75][76][77]. It was argued that the reason is that the trace of energy momentum tensor is not zero. The trace of energy momentum tensor is not zero for the NYM black holes. However, we have Smarr relation for n = 4z + 1. We guess it may originate from some properties of hypergeometric functions which have no explicit form for the case n = 4z + 1. We hope to investigate this issue in the future and find a physical reason for the Smarr relation not being satisfied.

B. Equation of state
To compare the critical behavior of the NYM black hole with that of the Van der Waals fluid, we first obtain the equation of state P (T, v, β) ≡ P , using equation (29). The critical points may be determined by using the following conditions We denote the volume, temperature and pressure of the critical points by v c , T c and P c , respectively. In the following, we will discuss the critical behavior of the NYM black hole for three types BIYM, ENYM, and LNYM, separately: Critical behavior of the BIYM solutions By substituting the relation (41) in Eq. (29), one can find the equation of state for the BIYM black hole, The critical behavior of NYM and nonlinear electrodynamics black holes are the same in four dimensions. The critical behavior of the BIM black hole in four and higher dimensions are in Refs. [56,57]. In this section, we study the critical behavior of NYM black holes in higher dimensions. If we consider the spherical case with k = 1 and impose the conditions (45) on the equation (46), we arrive at a cubic equation for the critical points where x, p and q are given by x is in terms of v c , therefor to have a positive value for v c , we have the condition To obtain an expression for the critical volume, we have to find the roots of equation (47). For equation (47), with p < 0 and real q, one can find one or three real roots. It has three real roots when 4p 3 + 27q 2 ≤ 0, which leads to We can write the roots in trigonometric form as where just x 0 and x 1 give a physical value for the critical volume v c in Eq. (48). x 0 and x 1 also satisfy the condition (49) which provides where . There is also one real root when β < β 0 , for which v c is negative and so the critical point can exist only for β ≥ β 0 . In the range of β 0 ≤ β ≤ β 2 , which we call the 'BI regime', there are two critical points corresponding to x = x 0,1 in equation (51). Although the critical temperature T c is positive for both critical points v c , one can find that for β ≥ β 1 = one of the critical points in the BI regime has negative pressure and thus it is unphysical. We have shown the P-v isotherms in the range of β 0 ≤ β ≤ β 1 and β 1 ≤ β ≤ β 2 for n = 5 in Figs. (8) and (9), respectively. For values of β greater than β 2 , there is also one critical point corresponding to x = x 1 in equation (51). In this range (β > β 2 ), which we call the 'YM regime', the critical behavior is identical to YM-AdS black hole. In other words, there is just one inflection point for T < T c and P-v isotherms are qualitatively identical to that of a Van der Walls fluid. We have depicted the critical behavior in this range (β > β 2 ) in Fig. (7) for n = 5. For a better understanding, we have brought the above results for n = 5 in Table (I). We can show that as β → ∞, the critical volume determined from x 1 reduces to the critical volume of the YM-AdS black hole. So, we name the branch determined from x 1 as the YM branch and the branch determined from x 0 as the BI branch, The behaviors of the critical values, T c , v c and P c with respect to the nonlinear parameter β are depicted in Fig.  (10) and Fig. (11). For large β, the critical values expand as We can observe that in the limit β → ∞, the critical values asymptote to those of the YM-AdS black hole. In this limit, the critical ratio tends to ρ c → 3/8, independent of the dimension n. It is in contrary to the abelian Maxwell theory in which the critical ratio, ρ c = Pcvc Tc depends on the dimension n [57].

Critical behaviors of the ENYM and LNYM solutions
The equations of state for the ENYM and LNYM solutions are defined respectively as and The critical behaviors of the ENYM and LNYM black holes are qualitatively the same as the BIYM case. However, the critical points and thus the branches mentioned in relation (53) cannot be determined analytically. Therefore, we solved Eq. (45) numerically in order to obtain the critical points. In the following, we will discuss the critical behavior just for the LNYM case. We have obtained the critical values of the six-dimensional LNYM solutions for different values of β in Table (II (12(b)). In the range of β 0 ≤ β ≤ β 2 , there are two critical points as we see in Fig. (12(a)), but for β > β 2 , there is only one critical point in Fig. (12(b)) and the P-v isotherm is analogous to the Van der Waals fluid.

C. Gibbs Energy
To get more information about the phase transition of the NYM black holes, we analyze the behavior of the Gibbs free energy, G. It can be achieved from the relation G(T, P ) = M − T S. We have plotted G versus T for the five-dimensional BIYM solution in Figs. (13) and (14). We can observe that the Gibbs free energy depends on the value of the nonlinear parameter β, as in the case of critical behavior. For β > β 2 , which we called the YM-branch, the phase transition is the same as the one for YM-AdS black holes or the RN-AdS black holes [56]. Namely, there is one critical point that corresponds to a phase transition from a small black hole to a large black hole when P < P c1 , see Fig. (13).
We have also shown the behavior of G in the BI branch in Figs. (14(a)) and (14(b)). In Fig. (14(a)) with the range of β ∈ (β 0 , β 1 ), there are two physical (with positive pressure) critical points described by T c1 and T c2 . The phase FIG. 14. G-T diagram of BIYM theory in n = 5. We have set e = k = 1. In the range of β0 ≤ β ≤ β1, there are two critical points with positive pressure. The Gibbs energy is not minimized at P = Pc2, and so the first order phase transition only occurs at P = Pc1. There is also a reentrant LBH/SBH/LBH transition for Pt ≤ P ≤ Pz. For β ∈ (β1, β2), there is only one physical critical point and so a first order small-large phase transition occurs for Pt ≤ P ≤ Pc1. We also observe a reentrant phase transition for the range of P ∈ (Pt, Pz).
transition at T = T c2 is not physical, since the Gibbs energy is not globally minimized at this point. However, there is a first-order phase transition from a small to a large black hole for T < T c1 which ends at T = T t . On the other hand, in a specific range of T ∈ (T t , T z ), one can observe a discontinuity in the global minimum of Gibbs energy which shows a reentrant LBH/SBH/LBH (large black hole/small black hole/large black hole) phase transition. Interestingly, this is a special significance for the higher-dimensional NYM solutions while there is no reentrant phase transition for the BI RN-AdS black holes in higher dimensions [56]. For the range of β 1 ≤ β ≤ β 2 in Fig. (14(b)), there is only one critical point with positive pressure at T = T c1 and a first order SBH/LBH phase transition occurs for T t ≤ T < T c1 . Similar to the previous case β ∈ (β 0 , β 1 ), in a specific range of T ∈ (T t , T z ), the global minimum of G is not continuous, which represents a reentrant LBH/SBH/LBH phase transition.
For the LNYM and ENYM cases, the Gibbs energy behavior is qualitatively the same as the BIYM case, so we do not probe them.

D. Critical exponents
Let us now calculate the critical exponents of the NYM black hole. These values can describe the physical quantities' behavior in the vicinity of the critical point. It should be noted that we use just the physical critical point that we named in the previous sections. There is only one physical critical point for each ranges β > β 2 , β 0 ≤ β ≤ β 1 and β 1 ≤ β ≤ β 2 . To calculate the exponent α, we should study the behavior of the heat capacity at constant volume. So, we write the entropy S as a function of T and V . If we use Eqs. (30) and (42), we can obtain where it shows that the entropy is independent of the temperature T and With defining the variables p, ν and τ as below the Eqs. (46), (58) and (59) can be rewritten as where ρ c = Pcvc Tc and Now if we expand Eq. (63) near the critical point, τ = t + 1 and ν = (1 + ω) 1/z , then we can obtain where B = 1 zρc , C = 1 and z = 3 for all dimensions. The derivative of Eq. (65) with respect to ω for t < 0 gets to dP = −P c (Bt + 3Cw 2 )dw.
By imposing the Maxwell's equal area law and considering a constant value for the pressure during the transition, we conclude that where ω s and ω l denote the volume of the small and large black holes. The non-trivial solution of Eq. (67) is obtained only for BCt < 0 as below As it is not possible to find an analytic result for the quantity C, so we gather the numeric results of A, B, and C for the BI case in Table(III). The results show that η = V c (ω l − ω s ) = 2V c ω l ∝ √ −t and so we get to the following result We can differentiate of Eq. (65) with respect to V and obtain the isothermal compressibility, t . This leads to the critical exponent γ = 1.
The critical isotherm t = 0 in Eq. (65) reduces to p − 1 = −Cω 3 , where it results to the "shape" The obtained results for α, β ′ , γ, and δ indicate that the critical exponents of the NYM black hole are exactly the same as those of the Van der Waals fluid [69]. In this section, we would like to investigate the dynamical stability of the nonlinear Yang-Mills black hole. Regge and Wheeler were the first who studied the stability by perturbation modes behaviors [78]. They decomposed the perturbations of 4-dimensional static and spherically symmetric background into odd-and even-parity sectors under a two dimensional rotation transformation. Dynamical stability of the nonlinear black hole solutions in Einstein gravity has been studied in Ref. [79]. Now we follow this paper to study the dynamical stability of the NYM black holes. We derive the dynamical stability of the four-dimensional NYM solutions. If one substitutesF ≡ 1 4 F 2 in the Lagrangian (2) and also define L(F ) ≡ −L(F )/4, then the Hamiltonian of the NYM Lagrangian is specified by H ≡ 2LFF − L where LF ≡ ∂L ∂F . It is proper to study the dynamical stability in the so called P frame where P ≡ L 2 FF . If H P (where H P ≡ ∂H ∂P ) vanishes nowhere outside the horizon, then the solutions are dynamically stable under odd parity perturbations. We obtain We evaluate H P in four dimensions. Using Eq.(3) in four dimensions, we get toF = e 2 2r 4 that results to the positive value for H P . We can conclude that the NYM solutions are dynamically stable under odd type perturbations since H P vanishes nowhere outside the horizon. For even-type perturbations, the solutions are dynamically unstable if the condition H xx > 0 is stablished that x = −2Q 2 P (Q is the Yang-Mills charge in Eq.(35)), and H xx = ∂ 2 H ∂x 2 . For the NYM black holes, H xx is obtained as follows We note that ENYM and LNYM black holes are unstable under even-type perturbations for respectively the regions β 2 < e 2 /r 4 and β 2 < e 2 /4r 4 .

VII. CONCLUDING RESULTS
In this paper, we attained a new n-dimensional analytic black hole solution for the non-abelian Yang-Mills gauge theory in the presence of three nonlinear Lagrangians in action, Born-Infeld, exponential and logarithmic Lagrangians. Using the Wu-Yang ansatz, we chose a set of gauge potentials with SO(n−1) and SO(n−2, 1) gauge symmetric groups and obtained the spherical and hyperbolic solutions. In four dimensions, the nonlinear Yang-Mills (NYM) solutions are similar to the nonlinear electrodynamics black hole solution in Maxwell theory and so the Yasskin theorem is established. We can also conclude that there is a transformation from the non-abelian gauge fields to a set of the abelian ones in n = 4. However, for n = 4, we could achieve a new set of nonlinear non-abelian Yang-Mills solutions. As expected, the nonlinear Yang-Mills solutions reduce to the Einstein-Yang-Mills ones when β → ∞.
We checked out the physical structure of the NYM black holes and observed that there is an essential singularity located at r = 0. Based on the behavior of the metric function in the limit r → 0, we found two different types of solutions for the horizon (Schw-type and RN-type) in n = 5. In some higher dimensions, the marginal mass is negative and so the Schw-type is the only solution.
We also calculated the thermodynamic quantities of the NYM black hole and probed the first law of thermodynamics. We analyzed the thermal stability of the solutions in both the canonical and the grand canonical ensembles. The stable regions happen for the NYM-AdS solutions with k = −1, 1 and for the flat ones with k = 1. The physically stable regions decrease as the dimension n increases. In the grand canonical ensemble with the mentioned parameters, stability emerges just for n = 4, and there are no stable regions for n = 6.
Furthermore, we investigated the critical behavior of the NYM black holes. We could access a Smarr-type formula for the dimensions n = 4z + 1, where z is an integer parameter. For the spherical solution with k = 1, we got to the exact critical points for the BIYM black hole, while we used a numeric method to obtain the critical points of the ENYM and LNYM black holes. There is an interesting range of β ∈ (β 0 , β 2 ), in which the critical behavior differs. One can find a discontinuity in the Gibbs energy for β ∈ (β 0 , β 2 ), which indicates a reentrant phase transition, and it is known as the zeroth-order phase transition. This reentrant phase transition only occurs in four dimensions for the nonlinear black hole in Maxwell theory [56]. However, in the case of NYM black holes, the reentrant phase transition occurs in four and higher dimensions. For β → ∞, the critical ratio goes to a constant value of 3/8, independent of the dimension n. This is one of the differences between the NYM black hole and the nonlinear electrodynamics one, in which the critical ratio depends on the dimension n [57]. We also calculated the critical exponents for the NYM black holes and found that the critical exponents are the same as those in the Van der Waals system. We also probed the dynamical stability of the NYM black holes. For odd-type perturbations, the NYM black holes are dynamically stable. Under even-type perturbations, we face instability for the ENYM and LNYM black holes in the regions β 2 < e 2 /r 4 and β 2 < e 2 /4r 4 , respectively.
B. The metric function f (r) for n = 5 and n = 9 For n = 5, the solution f (r) in Eq. (15) where E i (a, z) = z a−1 Γ(1−a, z) and Γ(a, x) and γ are respectively the gamma function and Euler-Mascheroni constant. The solution for n = 9 is It should be noted that we have determined the integration constants in Eq. (15) for n = 5 and n = 9 by assuming correspondence of BIYM, ENYM and LNYM theories with Yang-Mills theory for large values of β.