Thermodynamic description and quasinormal modes of adS black holes in Born-Infeld massive gravity with a non-abelian hair

We construct a new class of asymptotically (a)dS black hole solutions of Einstein-Yang-Mills massive gravity in the presence of Born-Infeld nonlinear electrodynamics. The obtained solutions possess a Coulomb electric charge, massive term and a non-abelian hair as well. We calculate the conserved and thermodynamic quantities, and investigate the validity of the first law of thermodynamics. Also, we investigate thermal stability conditions by using the sign of heat capacity through canonical ensemble. Next, we consider the cosmological constant as a thermodynamical pressure and study the van der Waals like phase transition of black holes in the extended phase space thermodynamics. Our results indicate the existence of a phase transition which is affected by the parameters of theory. Finally, we consider a massless scalar perturbation in the background of asymptotically adS solutions and calculate the quasinormal modes by employing the pseudospectral method. The imaginary part of quasinormal frequencies is the time scale of a thermal state (in the conformal field theory) for the approach to thermal equilibrium.

On the other hand, the existence of some limitations in the Maxwell theory motivates one to consider nonlinear electrodynamics (NED) [29][30][31][32][33][34][35][36][37]. Moreover, it was shown that NED can remove both the big bang and black hole singularities [38][39][40][41][42][43]. In addition, the effects of NED are important in superstrongly magnetized compact objects [44][45][46]. Considering GR coupled to NED attracts attention due to its specific properties in gauge/gravity coupling. Besides, NED theories are richer than the linear Maxwell theory and in some special cases, they reduce to the Maxwell electrodynamics.
One of the most interesting NED theories has been introduced by Born and Infeld [47,48] in order to remove the divergency of self energy of a point-like charge at the origin. The Lagrangian of Born-Infeld (BI) nonlinear gauge field is given by where β is BI nonlinearity parameter, F M = F µν F µν is the Maxwell invariant, F µν = 2∇ [µ A ν] is the Faraday tensor, and A ν is the gauge potential. Using the expansion of this Lagrangian for a large value of nonlinearity parameter leads to the Maxwell linear Lagrangian in which we receive the Maxwell Lagrangian at β → ∞. BI NED arises in the low-energy limit of the open string theory [49][50][51][52][53][54]. From the AdS/CFT correspondence point of view, it has been shown that, unlike gravitational correction, higher derivative terms of nonlinear electrodynamics do not have effect on the ratio of shear viscosity over entropy density [55]. Besides, NED theories make crucial effects on the condensation of the superconductor and its energy gap [56,57]. GR in the presence of BI NED has been investigated for static black holes [58][59][60][61][62][63][64][65][66][67][68][69], wormholes [70][71][72][73], rotating black objects [74][75][76][77][78][79], and superconductors [57,[80][81][82]. In addition, black hole solutions and their van der Waals like behavior in massive gravity coupled to BI NED have been studied in [14,83].
On the other hand, in addition to the Maxwell field, one can consider the non-abelian Yang-Mills (YM) field as matter source coupled to gravity. The presence of non-abelian gauge fields in the spectrum of some string models motivates us to consider them coupled to GR. In addition, the YM equations are present in the low energy limit of these models. Considering the YM field coupled with gravity violates the black hole uniqueness theorem and leads to hairy black holes. In such a situation, the field equations become highly nonlinear so that the early attempts for finding the black hole solutions in YM theory were performed numerically.
The purpose of this paper is obtaining the exact black hole solutions of Einstein-Massive theory in the presence of YM and BI NED fields (which is a more general solution compared with Reissner-Nordström, Einstein-Born-Infeld [112], Einstein-Yang-Mills [84], Einstein-massive gravity [13], Einstein-Born-Infeld-massive gravity [14] and etc.), and also, studying the thermal stability and phase transition of these black holes. Besides, we consider the massless scalar perturbations in the background of asymptotically adS solutions and calculate the quasinormal modes by employing the pseudospectral method. We investigate the effects of the free parameters on the quasinormal modes and dynamical stability. We also show that how the free parameters affect the time scale that a thermal state in conformal field theory (CFT) needs to pass to meet the thermal equilibrium.

II. FIELD EQUATIONS AND BLACK HOLE SOLUTIONS
Here, we consider the following (3 + 1)-dimensional action of EYM-Massive gravity with BI NED for the model where L BI (F M ) and F Y M = Tr F (a) µν F (a)µν are, respectively, the Lagrangian of BI NED (1) and the YM invariant. In addition, m is related to the graviton mass while f refers to an auxiliary reference metric which its components depend on the metric under consideration. Moreover, c i 's are some free constants and U i 's are symmetric polynomials of the eigenvalues of 4 × 4 matrix K µ ν = √ g µσ f σν which have the following forms where the rectangular bracket represents the trace of K µ ν . It is easy to obtain three tensorial field equations which come from the variation of action (3) with respect to the metric tensor g µν , and the gauge potentials A µ and A (a) µ as where ∧ D µ is the covariant derivative of the gauge field. The energy-momentum tensor of electromagnetic and YM fields, and also, χ µν can be written as In addition, the YM tensor F (a) µν has the following form µ is the YM potential and the symbols f (a) (b)(c) 's denote the real structure constants of the 3-parameters YM gauge group SU (2) (note: the structure constants can be calculated by using the commutation relation of the gauge group generators).
In order to obtain the spherically symmetric black hole solutions of EYM-Massive theory coupled to BI NED, we restrict attention to the following metric with the following reference metric ansatz [11] f µν = diag 0, 0, c 2 , c 2 sin 2 θ , where c is an arbitrary positive constant. Using the metric ansatz (12), U i 's reduce to the following explicit forms [11] Considering the field equations (5) with the following radial gauge potential ansatz one can obtain the following differential equation where E(r) = −h ′ (r) and prime refers to d/dr. Solving Eq. (15), we obtain where q is an integration constant which is related to the total electric charge of the black hole. It is clear that in the limit β → ∞, Eq. (16) tends to q/r 2 , and therefore, the Maxwell electric field will be recovered.
Hereafter and for the sake of simplicity, we use the position dependent generators t (r) , t (θ) , and t (ϕ) of the gauge group instead of the standard generators t (1) , t (2) , and t (3) . The relation between the basis of SU (2) group and the standard basis are t (r) = sin θ cos νϕt (1) + sin θ sin νϕt (2) + cos θt (3) t (θ) = cos θ cos νϕt (1) + cos θ sin νϕt (2) − sin θt (3) t (ϕ) = − sin νϕt (1) + cos νϕt (2) , and it is straightforward to show that these generators satisfy the following commutation relations In order to solve the YM field equations (6), just like the electromagnetic case, it is required to choose a gauge potential ansatz. Here, we are interested in the magnetic Wu-Yang ansatz of the gauge potential with the following nonzero components [97,101] A (a) where the magnetic parameter ν is a non-vanishing integer. It is easy to show that the chosen Wu-Yang gauge potential (19) satisfies the YM field equations (6). Using the YM tensor field (10) with Wu-Yang ansatz (19), one can show that the only non-vanishing component of the YM field is Considering the metric (11) with the electromagnetic (16) and YM fields (20), one can show that the only two different components of the field equations (4) are Since there is one common unknown function in both e tt and e θθ equations, it is expected to find that the mentioned field equations are not independent. After some manipulations, one can obtain the second order field equation by a suitable combination of first order one as and therefore, the solutions of e tt with an integration constant satisfy e θθ equation, directly. Solving Eq. (21), we can obtain the following metric function where is a hypergeometric function and m 0 is the only integration constant which is related to the total mass of black hole. Considering the obtained f (r), one finds that the fourth term is related to the magnetic charge (hair), the fifth term is related to the massive gravitons, and finally, the last term comes from the nonlinearity of electric charge. Now, it is worthwhile to investigate the asymptotic behavior of the nonlinearity parameter β on the solutions. Expanding the metric function (24) for large β one can recover the Maxwellian limit of the solutions. Therefore, for the massless graviton, m = 0, and linear electrodynamics, β → ∞, the metric function (25) reduces to the EYM solution with Maxwell field, as we expected.
On the other hand, for small values of the nonlinearity parameter (highly nonlinear solutions), we have which shows that the black hole is neutral at the highly nonlinear regime (β → 0). Considering Eq. (24), it is clear that the asymptotical behavior of the solutions is adS (or dS) provided Λ < 0 (or Λ > 0). In order to find the singularity of the solutions, one can obtain the Kretschmann scalar as which by inserting (24), it is straightforward to show that the Kretschmann scalar has the following behavior Equation (28) shows that there is an essential singularity located at the origin, r = 0. Moreover, the asymptotical behavior of the Kretschmann scalar for the large enough r confirms that the solutions are asymptotically (a)dS. Moreover, this singularity can be covered with an event horizon (for Λ < 0), and therefore, one can interpret the singularity as a black hole (Fig. 1). As a final point of this section, we should note that the metric function can possess more than two real positive roots which this behavior is due to giving mass to the gravitons (see [14,16] for more details).

A. Conserved and thermodynamic quantities
Here, we first obtain the conserved and thermodynamic quantities of the black hole solutions, and then examine the validity of the first law of thermodynamics.
The Hawking temperature of the black hole on the event (outermost) horizon, r + , can be obtained by using the definition of surface gravity, κ, where χ = ∂ t is the null Killing vector of the horizon. Thus, the temperature is obtained as It is worthwhile to mention that fourth term of RHS of Eq. (30) does not depend on the horizon radius, and therefore, one can regard it as a constant background temperature, T 0 = m 2 cc1 4π . As a result, we can investigate the solutions by using an effective temperature,T = T − T 0 .
The electric potential Φ, measured at infinity with respect to the horizon r + , is obtained by Since we are working in the context of Einstein gravity, the entropy of the black holes still obeys the so-called area law. Therefore, the entropy of black holes is equal to one-quarter of the horizon area with the following explicit form In order to obtain the electric charge of the black hole, we use the flux of the electric field at infinity, yielding It was shown that by using the Hamiltonian approach, one can obtain the total mass M in the context of massive gravity as [13] where m 0 comes from the fact that f (r = r + ) = 0. In the limit that the parameter β is small, one may think that the last term in Eq. (26) should contribute to the total mass along with m 0 (since they have the same radial dependence). But it is not correct since in order to calculate the total energy (mass) of spacetime, we have to regard the related action for large values of r. It is notable that the asymptotic behavior of the metric function for both limits r −→ ∞ and β −→ ∞ is the same. So, we should consider Eq. (25) to find the mass term, as a coefficient of r −1 in four dimensions. However, the functional form of the mass term should be proportional to r −1 for arbitrary values of r.
We should also note that the total mass is related to the geometrical mass, which is an integration constant of the gravitational field equation. If we regard the last term of Eq. (26) as a (piece of) mass-term, two problems are appeared; the first one is related to the higher-order series expansion of Eq. (26), in which they will be related to higher orders of mass with the same dimensional analysis, but they are not proportional to r −1 at all. The second one is related to the Smarr relation. It is straightforward to check that the Smarr relation is valid only for the mass related to the geometrical mass (m 0 ). Regarding the mentioned additional term, the Smarr relation is violated.
In addition, since the considered gravitational configuration has a time-like Killing vector, it is straightforward to calculate the energy (mass) of the system as the corresponding conserved quantity. Using the Arnowitt-Deser-Misner (ADM) method [113,114], one finds the mentioned conserved quantity for general solutions, f (r), is related to the geometrical mass, m 0 . Now, we are in a position to check the validity of the first law of thermodynamics. To do so, we use the entropy (32), the electric charge (33), and the mass (34) to obtain mass as a function of entropy and electric charge where We consider the entropy (S) and electric charge (Q E ) as a complete set of extensive parameters, and define the temperature (T ) and electric potential (Φ E ) as the intensive parameters conjugate to them Using Eqs. (32) and (33), one can easily show that the temperature (36) and electric potential (37) are, respectively, equal to Eqs. (30) and (31). Thus, these quantities satisfy the first law of thermodynamics On the other hand, the obtained black holes enjoy a global YM charge as well. In order to find this magnetic charge, we use the following definition In order to complete the first law of thermodynamics in differential form (38), one can consider the YM charge as an extensive thermodynamic variable and introduce an effective YM potential conjugate to it as an intensive variable which satisfies the first law of thermodynamics in a more complete way Regarding the differential form of the first law, it is worth mentioning that this equation may be completed by other additional terms, such as V dP in the extended phase space. In order to check the validity of the existence of such terms, one should check the first law in a non-differential form, the so-called Smarr relation. After some manipulations, one can find that where which confirm that the existence of additional terms and leads to a more complete form of the first law of thermodynamics It is worthwhile to mention that although it is possible to add C 2 dc 2 to the first law of thermodynamics (44) mathematically, we are not allowed due to the fact that all intensive and extensive thermodynamic parameters should appear in the Smarr formula (42). Therefore, we considered c 2 as a constant (not a thermodynamic variable) since it did not appear in the Smarr formula.

B. Thermal stability
In this section, we use the heat capacity for investigating the thermal stability of the obtained black hole solutions.
In this regard, one should consider the sign of heat capacity (its positivity and negativity) to study the stability conditions. The root of heat capacity (or temperature) represents a bound point. This point is a kind of border which is located between physical black holes related to the positive temperature and non-physical ones with a negative temperature. On the other hand, in our case, both divergence points of the heat capacity indicate one thermal phase transition point where black holes jump from one divergency to the other one. Besides, the heat capacity changes sign at such divergence points. So, one can conclude that the divergence point is a kind of bound-like point which is located between unstable black holes with negative heat capacity and stable (or metastable) ones. Therefore, it is logical to say that the physical stable black holes are located everywhere that both the heat capacity and temperature are positive, simultaneously.
Here, we study the thermal stability of the asymptotically adS solutions with Λ < 0. The heat capacity at constant electric and YM charges is given by where T has been obtained in Eq. (30). Considering (32), (33) and (35), one can easily show that the denominator of heat capacity is We recall that thermal stability criteria are based on the sign of heat capacity and it may change at root and divergence points. Therefore, it is necessary to look for the root and divergence points of the heat capacity at the first step. But unfortunately, because of the complexity of Eq. (45), it is not possible to obtain the root and divergencies of the heat capacity, analytically. So, we adopt the numerical analysis to obtain both bound and thermal phase transition points.
Before applying the numerical calculations, we are interested to clarify the general behavior of the heat capacity and temperature for the small and large black holes. For the fixed values of different parameters, there could exist two special r + 's, say r + min and r + max (see Fig. 2). The small black holes and large black holes are located before r + min and after r + max , respectively. The region of r + min < r + < r + max belongs to the intermediate black holes. Using the series expanding of (45), one obtains Considering Eq. (47), it is clear that for sufficiently small r + , the heat capacity and temperature are negative, and therefore, we have an unstable and non-physical black hole. Whereas from Eq. (48), we find that for large r + , both heat capacity and temperature are positive and there exists stable and physical black hole. In other words, Eqs. (47) and (48) confirm that the small black holes (r + < r + min ) are unstable and non-physical, whereas the large black holes (r + > r + max ) are physical and enjoy thermal stability. It is notable that in a special case there is just one specific horizon radius, r +s . In this case, we have unstable black holes for r + < r +s and stable ones for r + > r +s . However, it is not possible to identify this last property analytically, but we show it in Fig. 2 (see continues line). Now, we back to the numerical analysis of the heat capacity. Although we studied the general behavior of the heat capacity for the small and large black holes, the numerical calculations help us to classified the intermediate black holes (r + min < r + < r + max ). However, we are not going to study all possible behaviors of the heat capacity (because they contain different cases due to lots of free parameters) and just take some interesting ones. Figure 2 shows some different possibilities for the heat capacity. Clearly, this figure confirms that the small black holes are unstable (Eq. (47)) and large black holes are stable (Eq. (48)). According to the numerical analysis, we find that the heat capacity contains (i) only one bound point, (ii) one bound point and two divergencies, and (iii) three bound points and two divergencies. In the first case, we have unstable and non-physical black holes before the bound point (r +s ), but after this point, stable and physical black holes are presented. It is worthwhile to recall that from Eqs. (47) and (48), we expected such behavior. In the second case, we have stable and physical solutions between the bound point and smaller divergency. There are physical and unstable black holes between two divergencies. It is notable to mention that the large black holes are stable and physical as well. As for the last case, there are stable and unstable solutions respectively before and after the larger divergency.     In addition, we investigate the effects of different parameters on the bound points and divergencies of the heat capacity in tables I − III. From the table I, we find that the specific horizon radius, r +s , increases as the electric (magnetic) charge of black hole increases too. This could happen when the black hole absorbs electric (magnetic) charge. As a result, the region of unstable black holes increases. When the nonlinearity parameter increases and the nonlinear theory tends to the Maxwell case (2), the critical horizon radius increases. On the contrary, r +s is a decreasing function of the graviton mass (m). So, by increasing m, the region of unstable black holes decreases. Considering table II, it is clear that the smaller root (r + min ) and smaller divergency are decreasing functions of m, but the larger divergency (r + max ) increases as the massive parameter increases. In addition, we have found the same effects for β, q, and ν, but opposite behavior is seen for m. Table III shows that the smaller root (r + min ), the smaller divergency, and middle root decrease as the massive parameter increases, whereas the larger divergency and the larger root (r + max ) are increasing functions of m. Like case (ii), one can see the same behavior for β, q, and ν. The smaller divergency, r + min , and r + max are increasing functions of these parameters (m, β, q, and ν), but the middle and larger divergency are decreasing functions of them. Based on these three tables, we conclude that the qualitative effects of β, q, and ν on the heat capacity are quite the same.

C. P − V criticality in the extended phase space
It is well known that the most black holes can undergo a van der Waals like phase transition when one considers the cosmological constant as a thermodynamic pressure. In this section, we employ this analogy between the cosmological constant and pressure in the canonical ensemble (fixed Q E , Q Y M , β, and c 1 ) to investigate the P − V criticality and study phase transition of obtained black holes in extended phase space. Using the temperature given in Eq. (30) and the relation of P = −Λ/8π, it is straightforward to show that the equation of state is given by  and we made this choice in order to have a unique critical temperature (see appendix for more details). The thermodynamic volume is an extensive parameter which is conjugated to the pressure and has the following form where H is the enthalpy of the system. In this perspective, the total mass of black hole plays the role of enthalpy instead of internal energy due to the fact that the cosmological constant is not a fixed parameter anymore and it is actually a thermodynamic variable. Therefore, the thermodynamic volume is calculated as Hereafter, we use r + instead of V as a thermodynamic variable since it is proportional to the specific volume [115,116]. In order to study the phase transition of the black holes, we need to obtain the Gibbs free energy. In this extended phase space, one can determine the Gibbs free energy by using the following definition where H 1+ = H 1 (r = r + ). In addition, using the properties of inflection point and after some manipulations, we obtain the following equation Considering this equation, we find that it is not possible to obtain the critical horizon radius, r +c , analytically. As a result, we will not be able to calculate, analytically, the other critical parameters as well. So, we use the numerical analysis in order to study the van der Waals like phase transition of the black holes. In addition, we use such numerical analysis for investigating the effects of different parameters on the critical quantities.
Paul Ehrenfest has categorized the phase transition of thermodynamical systems based on the discontinuity in derivatives of the Gibbs free energy. The order of a phase transition is the order of the lowest differential of the Gibbs free energy that shows a discontinuity at the phase transition point. Thus, in a first order phase transition, there    exists a discontinuity in the first derivative of G (the entropy or volume). Next, in a second order phase transition, the entropy or volume becomes a continuous function and the heat capacity which is given by shows a sharp spike. Clearly, Fig. 3 confirms that the black holes under consideration enjoy the first order phase transition for temperatures and pressures less than their critical values and they undergo a second order phase transition at the critical point. For instance, we plot P − r + isotherm, G −T , and P −T diagrams for some fixed parameters to show the general phase transition behavior of the solutions (Fig. 4). Considering Fig. 4, we find that the obtained black holes have a van der Waals like phase transition between small black holes (SBH) and large black holes (LBH), and therefore, they enjoy a first order SBH-LBH phase transition. In this figure, P − r + isotherms show SBH area on the left, SBH+LBH coexistence area in the middle, and LBH area on the right. The dotted curve is a boundary between the regions of SBH, SBH+LBH, and LBH in the P − r + diagram. For temperatures above the critical temperature, there is no physical distinction between SBH and LBH phases, and this area is denoted as the supercritical region. In addition, in the G −T diagram, the phase transition point is located at the cross point, where SBH+LBH are presented, and black holes always choose the lowest energy. Moreover, the P −T diagram indicates the coexistence line between SBH and LBH which terminates at the critical point. The critical point is located at the topmost of the coexistence line with P = P c , r + = r +c , andT =T c . If black hole crosses the coexistence line from left to right or top to bottom, the system goes under a first order phase transition from SBH to LBH. Above the critical point, SBH and LBH are physically indistinguishable which is denoted by supercritical region.
From the left panel of Fig. 5, one can see that the red dashed (solid green) line corresponds to the negative (positive) heat capacity at constant pressure, C P , in the right panel. In addition, the divergencies of C P is indicated by two small black points A and B in the G −T diagram. The path bounded by these points is unconditionally unstable, but the paths A − C and B − C are metastable. Equivalently, in C P diagram, the region between point C 1 (C 2 ) and smaller (larger) divergency is metastable, and SBH-LBH phase transition does occur between C 1 and C 2 . This figure shows that during the phase transition from SBH to LBH, the heat capacity of the system increases. Moreover, this figure confirms that in order to have SBH-LBH phase transition, a local instability in the heat capacity is required.
In addition, Fig. 6 shows that the generalization of Einstein-Maxwell black holes into massive gravity and YM theory has a significant effect on the Reissner-Nordström black holes. In this theory, the region of SBH and LBH increases, and therefore, there is van der Waals like phase transition for higher temperatures and pressures compared with Reissner-Nordström black holes. In order to study the effects of different parameters on the critical points, we take table IV based on the numerical analysis. It is worthwhile to mention that by increasing the critical temperature and pressure, the region of SBH and LBH increases, and therefore, the region of phase transition increases too. From table IV , we find that the critical horizon radius is a decreasing function of the massive parameter and the critical temperature and pressure are increasing functions of this parameter. Considering table IV , one can see opposite behavior for the other parameters such as q, β, and ν. In other words, the critical horizon radius is an increasing function of these parameters, whereas the critical temperature and pressure are decreasing functions of them.

IV. ADS/CFT CORRESPONDENCE
In this section, we are going to point out two applications of the obtained solutions in the context of the adS/CFT correspondence. The adS/CFT correspondence relates string theory on asymptotically adS spacetimes to a conformal field theory on the boundary [117]. It is well-known that this holographic correspondence between a quantum field theory and a gravitational theory can be extended to explain some aspects of nuclear physics [118]. In addition, some phenomena like the Nernst effect [120,121], superconductivity [122], Hall effect [119] and the decaying time scale of perturbations of a thermal state in the field theory [123] have dual descriptions in gravitational theory.

A. Holographic superconductors
Here, we give some tips regarding the holographically dual superconductors of the Lagrangian (3). First of all, one should note that at the boundary (r → ∞), the metric function (24) tends to and we find that the nonlinearity parameter β does not play a significant role in the conductivity. Therefore, the BI NED can be replaced by Maxwell electrodynamics and the proper Lagrangian takes the following form In this case, due to the presence of Maxwell and YM fields, there are two options to investigate the holographic superconductors based on perturbing either Maxwell field or YM field. If we perturb the Maxwell (YM) field, the YM (Maxwell) field can be considered as an extra filed that is added to the Lagrangian as a matter source. If one wants to choose the Maxwell field to investigate the holographic superconductors, the case will be very similar to [11] (except the extra YM field they are the same) and it can be followed. Otherwise, if the YM field is preferred to describe the conductivity, the SU (2) gauge group should break down to the gauge symmetry U (1) 3 generated by the third component t 3 of the gauge field SU (2) [124] (see also [125,126]). Thus, the electromagnetic U (1) gauge symmetry is identified with the abelian U (1) 3 subgroup of the SU (2) group. Therefore, U (1) 3 is interpreted as the gauge group of electromagnetism which is considered in the boundary theory and Maxwell electrodynamics is an extra field.

B. Quasinormal modes
In terms of the adS/CFT correspondence, a large black hole in adS spacetime corresponds to an approximately thermal state in conformal field theory. Scalar perturbations of the black hole correspond to perturbations of this state. Thus, the decay of the scalar field describes the decay of perturbations of this thermal state. Therefore, we can calculate the time scale for the approach to thermal equilibrium by calculating the quasinormal modes (QNMs) of a large static black hole in asymptotically adS spacetime. Here, we shall obtain the QNMs of constructed black hole solutions to find the stability time scale of the corresponding thermal state. The other advantage of calculating the QNMs is investigating the dynamical stability of obtained black hole solutions undergoing scalar perturbations.
In order to calculate the QNMs, one can follow either Horowitz-Hubeny approach [123] or pseudospectral method [127]. The first one is based on Fröbenius expansion of the modes near the event horizon and forcing the differential equation to obey the boundary condition at the horizon. The second method replaces the continuous variable by a discrete set of points and solves the resulting generalized eigenvalue equation. However, we follow the pseudospectral method and use a public code presented in [128] to calculate the QN modes.
We now consider the fluctuations of a massless scalar field in the background spacetime of obtained black holes. In order to use the pseudospectral method, it is convenient to obtain the master equation in Eddington-Finkelstein coordinates. In these coordinates, the background line element takes the form , L is the adS radius related to the cosmological constant by Λ = −3/L 2 , and u = 1/r. Thus, u = 0 corresponds to the boundary and u = 1 represents the horizon. The equation of motion for a minimally coupled scalar field is governed by the Klein-Gordon equation It is convenient to expand the scalar field eigenfunction Φ in the form where Y lm (θ, ϕ) denotes the spherical harmonics. Substituting the scalar field decomposition (61) into (60) leads to the following second-order differential equation for the radial part in which ℓ is the multipole number and ω = ω r − iω i is the QN frequency with an imaginary part ω i giving damping of perturbations and a real part ω r giving oscillations. Therefore, in terms of the adS/CFT correspondence, τ = 1/ω i is the time scale that the thermal state needs to pass to meet the thermal equilibrium. On the other hand, the negativity of the imaginary part guarantees the dynamical stability of the black hole [123]. Otherwise, the perturbations increase in time and the spacetime becomes unstable. Causality requires ingoing modes at the event horizon and finite modes at spatial infinity that results in a discrete spectrum of frequencies ω. In order to analyze the behavior of modes ψ (u) near the horizon and the spacial infinity, we set r + = 1 and replace the value of M by considering f (r + ) = 0 without loss of generality. Starting with the horizon, by substituting an ansatz ψ (u) = (1 − u) p in (62), we find two solutions as ψ in (u) ∝ Const and ψ out (u) where Ω = ω/ (2πT ). By considering the time dependence e −iωt , the ψ out (u) behaves as In order to keep a constant phase, 1 − u has to increase as t increases, and thus u should decrease which means that this solution is outgoing. Therefore, we must consider just the ingoing solution ψ in (u) ∝ Const. There are two solutions near the event horizon; a normalizable mode ψ (u) ∝ u 3 and a non-normalizable one ψ (u) ∝ Const. If we rescale ψ (u) = u 2ψ (u), then the normalizable mode tends to zero linearly, whereas the non-normalizable mode diverges as ∼ u −2 . Doing this redefinition, the wave equation (62) becomes where Now, the normalizable mode behaves smoothly at the boundary and it should be considered, while we discard the other solution. The wave equation (64) is an input for the code and one can fix the free parameters and the event horizon radius r + = u −1 + to calculate the QN modes. In the previous section, it was shown that the small black holes are unstable and non-physical, whereas the large black holes are physical and enjoy thermal stability (see Eqs. (47) and (48)). On the other hand, the large black holes correspond to the thermal states in CFT. Thus, we shall focus on the QNMs of large black holes (r + >> L) and discard the small ones (r + << L) for L = 1 as the adS radius.
Here, we set the free parameters as q = 1, c = 1, c 1 = −1, c 2 = 2, and ℓ = 0, and evaluate the QNMs for different values of m, β, ν, and r + . In table V , we list the QNM frequencies for the fundamental mode (n = 0) and the first overtone (n = 1) of intermediate black holes (r + = 5, 10) and large ones (r + = 50, 100). From this table, one can see that as the overtone number and the event horizon radius increase, both the real and imaginary parts of frequencies increase as well. But an opposite behavior is seen for increasing in the graviton mass. Besides, the real (imaginary) part of the frequencies decreases (increases) when the magnetic charge increases. We recall that the nonlinearity parameter β does not play a significant role at the boundary r → ∞ (u → 0), and thus increasing/decreasing in β does not change the value of QNMs as it can be seen from the table. Therefore, the BI NED can be replaced by Maxwell electrodynamics when we want to investigate the applications of the solutions in the context of the adS/CFT correspondence. We should mention that as the imaginary part of frequencies increases, the corresponds thermal state meets the stability faster. In addition, the obtained black hole solutions undergoing massless scalar perturbations are dynamically stable since all the frequencies have a negative imaginary part. m β ν r + = 5 r + = 10 r + = 50 r + = 100 It is worthwhile to mention that as r + increases, changing in ν does not affect the QNMs significantly (compare the first line and last line for r + = 50, 100 in table V ). But this is not correct in the case of m (compare the first line and second line for r + = 50, 100). In order to explain this fact, one may consider the temperature (30) for large black holes at the first step and secondly, look at the relation between the QNMs and this temperature illustrated in Fig. 7. As one can see, both the real and imaginary parts of frequencies increase linearly with increase in the temperature (69). Therefore, changing in ν does not affect the QNMs since it is absent in (69), whereas m is present. From (69), we can find that increasing in c and c 1 leads to increasing in QNMs, but q and c 2 do not change the QNMs in the case of large black holes, as ν did not. The points in Fig. 7, representing the QNMs, lie on straight lines through the origin. For the real part, the lines are given by ω r = 7.747T, n = 0 ω r = 13.236T, n = 1 for m = 2, ω r = 7.752T, n = 0 ω r = 13.230T, n = 1 for m = 3, Re Ω , n 0 Im Ω , n 0 Re Ω , n 1 Im Ω , n 1 Re Ω , n 0 Im Ω , n 0 Re Ω , n 1 Im Ω , n 1 while for the imaginary part we have ω i = 11.158T, n = 0 ω i = 20.594T, n = 1 for m = 2, ω i = 11.151T, n = 0 ω i = 20.596T, n = 1 for m = 3.
In terms of the adS/CFT correspondence, τ = 1/ω i is the time scale for the approach to thermal equilibrium. Therefore, Eqs. (72) and (73) are the main results of this subsection. One may note that both the real and the imaginary parts of the frequencies are linear functions of r + since the temperature of large black holes is a linear function of r + . Interestingly, the same result was found for the Schwarzschild-adS black hole [123].

V. CONCLUSIONS
In this paper, we have obtained Einstein-Massive black hole solutions in the presence of YM and BI NED fields. We have also studied the geometric properties of the solutions and it was shown that there is an essential singularity at the origin which can be covered with an event horizon. In addition, we have calculated the conserved and thermodynamical quantities, and it was shown that even though the YM and BI NED fields modify the solutions, the first law of thermodynamics is still valid.
Moreover, we have studied the thermal stability of the obtained black holes and investigated the effects of different parameters on the stability conditions. We have found that the large black holes (r + > r + max ) are physical and stable, whereas the small black holes (r + < r + min ) are non-physical (T < 0). Furthermore, we have classified the medium black holes (r + min < r + < r + max ) in Fig. 2 and investigated the effects of different parameters on thermal stability of these black holes in tables I − III.
In addition, we have considered the cosmological constant as thermodynamical pressure and it was shown that the obtained black holes enjoy the first order SBH-LBH phase transition. Also, we have studied this kind of phase transition in the heat capacity diagram and specified the unstable and metastable phases of obtained black holes related to the negative and positive heat capacities, respectively. It was shown that during the phase transition from SBH to LBH, the heat capacity of the system increases. We have seen that the generalization of Reissner-Nordström solutions into massive gravity and YM theory increases the critical temperature and pressure, and as a result, the region of SBH and LBH increases. Moreover, we have investigated the effects of different parameters on the critical points, and we found that the parameters q, β, and ν have opposite effect on the critical points compared with the massive parameter, m.
Besides, we have considered massless scalar perturbations in the background of obtained black holes in asymptotically adS spacetime. We also have calculated the QN frequencies by using the pseudospectral method in order to investigate the dynamical stability of the black holes, the effects of different parameters on the QNMs, and obtain the time scale of the thermal state for the approach to thermal equilibrium in CFT. It was seen that the obtained solutions are dynamically stable and BI NED generalization does not affect the frequencies. Furthermore, it was shown that increasing in r + , c, c 1 , and m lead to increase in both the real and imaginary parts of the frequencies. It is worthwhile to mention that this result depends on the sign of c and c 1 (through the text, we considered a negative value for c 1 , and therefore, increasing in m has led to decrease in the QNMs). In addition, we have found that ν, q, and c 2 do not affect the QNMs in the case of large black holes. Since a static large black hole in adS spacetime corresponds to an approximately thermal state in conformal field theory, ν, q, and c 2 have no effect on the time scale of the thermal state. Just like the Schwarzschild-adS black holes [123], both the real and imaginary parts of frequencies for the large black holes were linear functions of the temperature.
As a final remark, it is worth mentioning that although we consider the ADM mass in the context of black hole thermodynamics, there is another extension of mass (so-called hairy mass) for hairy black holes which is related to the calculation of the null circular geodesic (photon-sphere) [129,130]. Such a hairy mass is not related to our discussion in this paper and it can be considered as a new work with photon-sphere concentration.