In the last two decades, the gravitational wave searches by the LIGO Scientific and Virgo collaboration for the merger of compact binaries [1–5] and the first measurement of the shadow image captured by the Event Horizon Telescope [6–11] have provided convincing evidence in terms of the existence of black holes (BHs) in the center of many galaxies, such as the elliptical M87 galaxy and spiral Milky Way [12]. Actually, the concept of a BH exactly originates from one of the most important predictions of general relativity and, in view of their causal structure, BHs have a surface to which any wave or particle passing through this surface cannot return; this surface is labeled as the event horizon.

It is now well known that BHs are abundantly present in the Universe. Generally speaking, according to the mass of BHs, there are four types of BHs available for discussion [13]: (1) miniature BHs with masses of at least a few Planck masses or more; (2) stellar mass BHs (M ∼ 10 M_⊙), formed after the death of some normal stars; (3) intermediate mass BHs M ∼ 10³ M_⊙; and (4) supermassive BHs in the centers of galaxies with masses M ∼ (10⁵ − 10⁹) M_⊙ (such as Sagittarius A* and M87). The first type has its own motivation for discussion, and perturbations of the latter three types of BHs are important for astrophysical observations. On the other hand, via the no-hair theorem, BHs can be characterized for three parameters in general relativity, i.e. charge, mass and angular momentum. Within this line of thought, there are also four types of BHs available for discussion [14]: (1) Schwarzschild BHs, characterized only by mass; (2) Reissner–Nordström BHs, characterized by mass and electric charge; (3) Kerr BHs, characterized by mass and angular momentum (such as the fact that the Sgr A* BH and M87 central BH can be described by Kerr BHs [15]); (4) Kerr–Newman BHs, characterized by charge, mass and angular momentum. These four types are vacuum/electrovacuum solutions obtained within the framework of general relativity, and the description of BHs and these solutions contain an unavoidable curvature singularity in their interiors [16]. However, general relativity cannot overcome the existence of singularity in a very high curvature regime. Thus, families of regular BH solutions of Einstein's gravity have been found where the singularity is eliminated. In fact, such solutions cannot be vacuum solutions of the Einstein equations, but necessarily contain an appropriately chosen additional field (such as nonlinear electrodynamics), which guarantees the violation of the energy conditions related to the existence of physical singularities [17]. Afterward, Bardeen proposed the Bardeen regular BH [18]. As the Bardeen solution presented a non-vanishing Einstein tensor, a new explanation emerged, i.e. the metric could be interpreted as a solution of Einstein equations coupled with nonlinear electrodynamics with a magnetic charge [19]. It is well known that BHs in the real world are not isolated and are not embedded in empty backgrounds. In astronomy, stellar BHs are expected to be surrounded by matter or fluids, for example, the fact that the solutions of BHs obtained by Kiselev just describe such spacetimes surrounded by anisotropic fluids [20, 21]. Originally, the fluid surrounding the BH was called quintessence, whose equation of state is given by p = wρ_q, where p is the pressure, ρ_q is the density of energy and w is the state parameter. As shown in [22], the fluid can be interpreted as a mix of fluids, which contain an electrically charged fluid, an electromagnetic and/or a scalar field, and it is not interpreted as quintessence. Different values of w are responsible for different matter surrounding the BH; specifically speaking, w = 0 represents an ideal gas, w = 1/3 denotes the ultrarelativistic particles and negative values correspond to exotic matter. Moreover, Dariescu et al showed that the Kiselev geometry is actually an exact solution of the Einstein equations coupled to nonlinear electrodynamics [23].

In conclusion, since the regular Kiselev BH solution family overcomes the singularity problem theoretically and the celestial model described by it is more in line with the actual situation, this kind of BH solution family has received considerable attention from scholars [24–26]. We are also very interested in this kind of BH spacetime. Thus, we study the stability of BH spacetime in terms of the Bardeen–Kiselev BH with the cosmological constant, motivated by [27]; its exact solution can be obtained by coupling the Einstein equation and nonlinear electrodynamics. We notice that this solution contains a cosmological parameter, which involves a major topic in cosmology, i.e. the accelerating expansion of our Universe. As an explanation for this phenomenon, a large number of dark-energy models have been proposed, among which the one with the cosmological constant is the most famous. This model is reasonable in physical theory and consistent with most observational results. It is fair to say that the exact solution (Bardeen–Kiselev BH with the cosmological constant) contains many kinds of theories, which should be considered as describing a toy astrophysical model of the real world. Meanwhile, this solution contains a magnetic charge, which cannot be found in the Universe. However, according to [28], magnetically charged BHs may be produced in the early Universe, and they are more likely to retain the magnetic charge, thus avoiding neutralization by ordinary matter accreting on them, unlike electrically charged counterparts. Therefore, it is worth studying the quasinormal modes (QNMs) of electrically or magnetically charged BHs to correctly capture the parameters of realistic BHs from observations. Recently, the research on the QNM of BHs has been one of the hot spots [29–31]. As the ”characteristic sound” of BHs, the QNM is considered to be strong evidence of direct detection of the existence of BHs in recent years, which has been observed in gravitational wave experiments. In-depth research and understanding of various BH QNMs is conducive to our detection of the existence of BHs through gravitational wave experiments, and helps us to further understand the characteristics of BHs in observation. Another motivation comes from the famous duality between supergravity in anti-de Sitter spacetime (AdS) and conformal field theory (CFT), AdS/CFT correspondence [13].

The stability of BH spacetime is probably one of the most interssting issues in the context of general relativity, since a BH under certain perturbations can help one to study the nature of the BH itself. In general, by studying the BH merger or the field evolution of the BH background, one can discuss the stability of BH spacetime under perturbation [32, 33]. It is well known that a BH in perturbation can emit gravitational waves, which are dominated by QNM, and this mode is a complex frequency of excited oscillation modes [34], where the real part and the imaginary part of QNM's frequency are the oscillation frequency and damping of the BH, respectively. For electrically or magnetically charged BHs, the Einstein–Maxwell theory is usually considered to be the standard theory and leads to the well known Reissner–Nordström solution [35]. Note that the QNMs may be subject to corrections due to some effects of electrodynamics beyond the Maxwell theory, and which actually exist, at least when we consider the quantum electrodynamics. Based on these considerations, we are interested to discuss the QNMs of magnetically charged BHs within a framework beyond the standard Einstein–Maxwell theory. Usually, a BH in perturbation can be classified into three distinct stages; the first consists of an initial outburst of a wave, which depends completely on the initial perturbing field. The second one is usually a long period of damping proper oscillations, which are dominated by the QNM; this stage contributes to gravitational wave detection. The final stage is power-law tail behavior at very late times. In this work, we focus on the second stage of the evolution of perturbations represented by QNM. Generally, the aim of the beginning of QNM calculation is to reduce the perturbation equations into a two-dimensional wavelike form with decoupled angular variables and, if the variables are decoupled, the equation for time and radial variables will transform into a Schrödinger-like form in a stationary background. Then, the corresponding potential function can be determined, which is exactly the key for the computation of the QNM frequencies. The QNM can be calculated using various methods, such as the continued fraction method [36], the Wenzel–Kramers–Brillouin (WKB) method [37], Pöschl–Teller approximation [38], the time domain method [39] and the expansion method [33]. In view of the accuracy of the WKB method in terms of the real and the imaginary parts of the dominant QNMs with n ≤ l, we use this method in the QNM frequencies calculation. In particular, we notice that this method has been applied to the numerical calculation method for QNM frequency in terms of regular Kiselev BHs [24, 25] and Kiselev BHs [40, 41]. In our work, we focus on the behavior of the dynamical response of the spherically symmetric and magnetical BH, which represents the exact solutions of coupled Einstein gravity and nonlinear electrodynamics (NED) to small electromagnetic and gravitational perturbations. In particular, we intend to identify whether it would distinguish the BHs related to the NED from those BHs related to the standard linear electrodynamics due to their response to electromagnetic and gravitational perturbations. Perturbations of BHs imply the exploration of stability of their spacetime, and the stability of different BHs involving NED has been considered in [42, 43]. In addition, outside the event horizon, the effective potential plays the role of a filter [44]. Specifically, some waves probably pass through the effective potential and transmit to infinity, while some waves are probably reflected back by the effective potential. Thus, if we consider an observer at infinity, obviously the radiation received by the observer is different from the radiation emitted from the event horizon [45]. This behavior in Hawking radiation is known as greybody factors [46], which can encode information about the horizon structure of BHs theoretically and modify the quasinormal spectra experimentally [47]. Also, to estimate the transmission probability of radiation from a BH event horizon to its asymptotic region, it is necessary to study the greybody factors of perturbations.

The authors in [27] discussed the Bardeen solution with the cosmological constant surrounded by Kiselev quintessence. They showed that this solution can be obtained by Einstein equations coupled with nonlinear electrodynamics, and it is not always regular based on what the conditions for regularity are. They also analyzed the thermodynamics associated with this type of solution by establishing the form of the Smarr formula and the first law of thermodynamics. In this work, we study the stability of this solution and the QNMs of perturbations; the reflection coefficient, greybody factor and absorption coefficient of perturbations are also involved. This paper is organized as follows. In section 2, we present a brief review of the Bardeen–Kiselev BH with the cosmological constant. In section 3, the gravitational and electromagnetic perturbations in the Bardeen–Kiselev BH with the cosmological constant are reported. The QNMs of the Bardeen–Kiselev BH with the cosmological constant for gravitational and electromagnetic perturbations are analyzed, and we compare the frequencies with the Reissner–Nordström-de Sitter BH surrounded by quintessence (RN-dSQ) [48] in section 4. In section 5, we focus on a discussion on the greybody factor, the reflection coefficient and the total absorption cross section. Finally, the work is summarized in section 6. In this paper, we consider natural units, where c = G = ℏ = 1, and the metric signature (+, − , − , − ).

Here, we review [27] and briefly discuss some of the calculation processes, i.e. we consider a spherically symmetric spacetimewhere the solution is magnetically charged and has the cosmological constant. The action that describes this theory iswhere R and λ denote the scalar curvature and the cosmological constant (a famous dark-energy model), respectively, and the function of the electromagnetic scalar F of the nonlinear Lagrangian of electromagnetic theory isHere, c is the normalization constant related to the density of Kiselev quintessence, ${\rho }_{q}=-\tfrac{c}{2}\tfrac{3w}{{r}^{w+3}}$, M is the Arnowitt–Deser–Misner mass of the BH and q is the magnetic charge. With regard to the L(F) function, the first term of the L(F) in our paper is the same as the Lagrange density used by Ayon-Beato and Garcia aiming to obtain the Bardeen metric [19]. For the second term of L(F), according to [23], the Kiselev quintessence parameter c and w are in relation to the electric or the magnetic charge, and they emerge in the power-Maxwell Lagrangian. Thus, the electromagnetic field and gravitation are connected with the Kiselev quintessence parameters c and w. Their connection changes the spacetime metric function; therefore, a significant impact will emerge in the subsequent perturbation study. Moreover, it should be noted that equation (3) can transform into ordinary electrodynamics in the weak field limit [27]. The electromagnetic scalar F is given by [49]where $F=\tfrac{2{q}^{2}}{{r}^{4}}$ and F_μν is the electrodynamic field tensor that can be expressed in terms of a gauge potential as F_μν = ∂_μA_ν − ∂_νA_μ. The electromagnetic field tensor implies that F_μν is anti-symmetric and it has only six independent components. The only non-zero component of F_μν for a spherically symmetric spacetime that is only magnetically charged is ${F}_{23}=q\sin \theta $.

Here, one can express the covariant equations of motion asif we consider the stress-energy tensor to nonlinear electrodynamicswith ${L}_{F}=\tfrac{\partial L(F)}{\partial F}$; then, by solving the Einstein equations, we obtain the metric functionand when c > 0 and −1 ≤ w ≤ 0, the anisotropic fluid fulfills the null energy conditions [22]. As stated before, different values of w are responsible for different matter surrounding the BH. Since we are interested in asymptotically dS-like spacetimes, in this work, we will focus on the cases −1 ≤ w ≤ −1/3. And with regard to equation (7), we find the Kiselev–(anti-)de Sitter solution with the cosmological constant for the limit q → 0 and the Bardeen–(anti-)de Sitter solution for c → 0 [50].

The study of BH perturbations was first proposed by Regge and Wheeler in terms of the investigation of the odd parity type of the spherical harmonics [51]. After this, Zerilli generalized it to the even parity case [52]. It is well known that generally there are two different kinds of perturbations of BHs in the context of the general theory of relativity: the test field in a BH background, and the perturbation of the metric. The former is achieved using the method of solving the dynamical equation for the given test field in the background of the BH; the second one is obtained by linearising the Einstein equation to derive the evolution equations, i.e. the gravitational perturbation. When compared with the strength of the external fields decaying in the vicinity of the BH, the gravitational radiation is the strongest one. In this section, we firstly focus on the gravitational perturbation of the Bardeen–Kiselev BH with the cosmological constant. The general procedure for gravitational perturbation is to introduce a small perturbation (h_μv ≪ 1) into the static background metric $({\tilde{g}}_{\mu v})$; then we assume the perturbed background metric g_μv as

The perturbations h_μv can be decomposed asThis formalism is similar to the axial decomposition in the Regge–Wheeler gauge. Note that it is not exactly the same gauge; we prefer to obtain the form of ξ(θ) through the field equations rather than by imposition of some kind of expansion in terms of spherical harmonics [53].

Here, we have considered the perturbation in the energy momentum tensor and obtain Einstein's equation asandIt should be noted that [54] studied the QNMs of gravitational perturbation around a regular Bardeen BH surrounded by quintessence in vacuum, regardless of the perturbed energy momentum tensor; we may view that our consideration is more interesting and general. The perturbed Einstein equations then lead to the equalities in equation (12) with T_θφ = 0 and T_rφ = h₁Lξ(θ)where l is the multipole number, k² = 8π and $\xi (\theta )={P}_{{\ell }}(\cos \theta )$(the Legendre polynomials). By using the definition $\varphi (t,r)=\tfrac{f}{r}{h}_{1}(t,r)$, equation (12) can be expressed aswhere ′ denotes the derivative with respect to the radial coordinate r. To facilitate this procedure, we change the variable r to the tortoise coordinate r_* with the definition dr = f(r)dr_*, and when r → r_c, r_* → ∞, when r → r₊, r_* → − ∞, where r_c is the cosmological horizon and r₊ is the event horizon; then we havewith the effective potentialThis formalism is the same as was found in [55] in the limit λ → 0, but there is a difference in the numerical factor in the effective potential in terms of the coefficient of L and f.

Here, if we consider the temporal dependence as φ ∼ $\Psi$e^−iωt, we can get the master equation for gravitational perturbations of the BH aswhere the parameter ω denotes the dissipative modes in time. Thus, the BH will oscillate after the perturbation and then go back to a stable state, which is known as quasinormal frequencies.

Next, the behavior of the dynamical response of the spherically symmetric BH for electromagnetic perturbations is our topic, and we assume that electromagnetic perturbations do not alter the geometry of spacetime [49]. We then decompose the four-vector potential of the electromagnetic field as the unperturbed background potential ${\widetilde{A}}_{\mu }$ and perturbed part δA_μIn view of the static and spherically symmetric background, we consider the unperturbed four-vector potential of the magnetically charged BH asIn a spherically symmetric background, the perturbation in the vector potential can be expressed as a superposition of vector spherical harmonics, i.e.where Y_ℓm(θ, φ) denotes scalar spherical harmonics.

The non-vanishing covariant components of the electromagnetic field tensor of the four-potential equation (17) with the perturbation equation (19) are given byBy the relation F^μν = g^μαg^νβF_αβ, all non-zero contravariant components of the electromagnetic field tensor are as follows:In view of the fact that the field strength remains the same at the zeroth order but has components in the first order for the total four-vector potential, by combining equations (20) and (21), we present the electromagnetic field strength F aswhere the first term corresponds to the unperturbed electromagnetic field strength $\tilde{F}$, and the second term represents the contribution of the perturbation to the field strength δF, which leads to $F=\tilde{F}+\delta F$.

Near $\tilde{F}$, by using the Taylor series up to the first-order term for the expansion of L_F, we havewhere ${\tilde{L}}_{\tilde{F}\tilde{F}}={\partial }_{\tilde{F}}^{2}\tilde{L}={\partial }_{\tilde{F}}{\tilde{L}}_{\tilde{F}}$; note that $\tilde{F}$ and ${\tilde{L}}_{\tilde{F}}$ depend on r. By combining equation (21) into equation (5) we haveFor u = θ and u = φ, the above equation transforms intowhere ${\tilde{L}}_{\tilde{F}}^{{\prime} }$ denotes the first derivative of ${\tilde{L}}_{\tilde{F}}$ with respect to r. For convenience, we take a transformation $\psi (t,r)\,={\left({\widetilde{L}}_{\tilde{F}}\right)}^{-\tfrac{1}{2}}\varphi (t,r)$, and by considering the tortoise coordinate, we obtain the Schrödinger-like wave equationwith the effective potentialThis formalism is the same as that found in [32] and [49]. Also, there is a difference in the numerical factor in the effective potential.

For the gravitational and electromagnetic perturbations, we present the effective potentials V(r) for the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ in figure 1. Specifically, the overall nature of both potentials in the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ are the same, i.e. they are positive definite between the event and cosmological horizons, and they all have a single maxima. Besides, for a fixed set of parameters in electromagnetic perturbation, the height of the RN-dSQ potential is larger than the Bardeen–Kiselev one, which implies that the RN-dSQ has a smaller absorption coefficient than the Bardeen–Kiselev BH with the cosmological constant. Meanwhile, the opposing observation appears in gravitational perturbation. The figure shows that the potential decreases with increasing parameters ∣w∣ or λ, which indicates that the smaller value of ∣w∣ or λ suppresses the emission modes for gravitational and electromagnetic perturbations.

In this section, our purpose is to discuss the QNMs and the stability of the perturbations in the Bardeen–Kiselev BH with athe cosmological constant spacetime, especially for the QNMs of the Bardeen–Kiselev BH with the cosmological constant for gravitational and electromagnetic perturbations. QNMs for a perturbed BH spacetime are the solutions to the wave equation given in equations (16) and (26), and to derive these solutions, it is necessary to impose proper boundary conditions:However, it is tricky to analytically solve the time-independent second-order differential equation (such as equation (16)) with the potential (such as equation (15)) for a nonlinear magnetically charged BH with the cosmological constant surrounded by Kiselev quintessence. The WKB method can be used for the effective potential, which configurates the form of a potential barrier and takes constant values at the event horizon and spatial infinity. Specifically, this method is based on matching the asymptotic WKB solutions at spatial infinity and the event horizon with a Taylor expansion near the top of the potential barrier through two turning points. Schutz and Will first applied this method to the scattering problem around BHs [56]. Subsequently, Iyer and Will extended it to the third-order WKB approximation [37], Konoplya reported the sixth order [57] and Matyjasek and Opala recently provided the thirteenth order [58]. Now we use the WKB approximation method to numerically calculate the QNM frequencies of the electromagnetic and gravitational perturbations in the Bardeen–Kiselev BH with the cosmological constant. By using the considered potential functions, one can obtain the QNM frequencies via a sixth-order WKB method and, as seen in [59], this method is the most accurate one for finding the quasinormal spectrum with lower overtones. The BH potential V(r) in the present of the sixth-order formula iswhere, among them, V₀ is the maximum effective potential of V(r) at the tortoise coordinate r_*, n is the overtone number (we study the case n = 0) and the correction term Λ_i can be obtained in [37]. And generally speaking, the quasinormal frequencies ω take the form ω = ω_R − iω_I, where the real and the imaginary parts of ω denote actual field oscillation and damping of the perturbation, respectively.

We have numerically obtained the QNM frequencies for the gravitational and electromagnetic perturbations. In figures 2, 3 and 4, we have exploited the sixth order WKB approximation for calculating the frequencies of the gravitational perturbations of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ; these frequencies can be obtained by individually varying q, c and λ. In figure 2, for the same parameter space, one can see that the oscillation frequency of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ increases with the increasing magnetic charge q. For the Bardeen–Kiselev BH with the cosmological constant, the damping rate decreases as the charge increases, which implies that the decay of the modes is slower, and for the RN-dSQ, $-{Im}$ w increases with the increasing q. Therefore, for the same parameter space, with the increasing magnetic charge q, the RN-dSQ is more stable than the Bardeen–Kiselev BH with the cosmological constant. In figures 3 and 4, we observe similar qualitative behavior of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ, i.e. both the oscillation frequency and the damping rate decrease with increasing values of c and λ.

Figures 2 to 4 show that the oscillation frequency and the damping rate of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ increase with the increasing l and w. And by considering the effect of the state parameter w and the normalization constant c on the oscillation frequency and the damping rate of the Bardeen–Kiselev BH with the cosmological constant, we can conclude that, due to the presence of Kiselev quintessence, the gravitational perturbations of the Bardeen–Kiselev BH with the cosmological constant damp and oscillate more slowly, and these results are consistent with [24]. Moreover, by varying the BH parameters, we analyze the behavior of both the real and imaginary parts of the Bardeen–Kiselev quasinormal frequencies and compare these frequencies with the RN-dSQ. Interestingly, it should be noted that for the gravitational perturbations, the response of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ in terms of the imaginary part of w are different only when the charge parameter q is varied; this behavior can be used to understand nonlinear and linear electromagnetic fields in curved spacetime separately.

In tables 1, 2 and 3, we report the quasinormal frequencies in electromagnetic perturbations of the Bardeen–Kiselev BH with the cosmological constant and compare frequencies with the RN-dSQ. Specifically, in table 1, we present the QNM frequencies of electromagnetic field perturbation in the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ in terms of the changeable magnetic charge q. It shows that the response of the QNM frequencies of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ under electromagnetic perturbations are the same as the gravitational perturbation when the charge parameter q is changeable. We present the QNM frequencies of electromagnetic perturbations in the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ in terms of the changeable normalization factor c in table 2. It shows that the oscillation frequency of the two BHs decrease with the increasing values of c while, for the Bardeen–Kiselev BH with the cosmological constant, the damping rate increases with the increasing values of c; this implies that the decay of the modes in the Bardeen–Kiselev BH with the cosmological constant is faster, and the damping rate of the RN-dSQ decreases with the increasing values of c. In table 3, we present the QNM frequencies of electromagnetic field perturbation in the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ in terms of the changeable cosmological constant λ, which shows that both the oscillation frequency and the damping rate decrease with the increasing values of λ, and this behavior is the same as the gravitational perturbations. From these tables, for the same BH parameters, when the state parameter w increases, the real parts of the quasinormal frequencies of the Bardeen–Kiselev BH with the cosmological constant increase while the absolute values of the imaginary parts decrease, while the real parts of the quasinormal frequencies and the absolute values of the imaginary parts of the RN-dSQ increase. Also, when considering the effect of the state parameter w and the normalization constant c on the oscillation frequency and the damping rate of the Bardeen–Kiselev BH with the cosmological constant, we can conclude that, in the presence of Kiselev quintessence, the electromagnetic perturbations of the Bardeen–Kiselev BH with the cosmological constant damp much faster and oscillate more slowly. This behavior is different to the gravitational perturbation.

With regard to these figures and tables, by varying the BH parameters, we analyze the behavior of both the real and imaginary parts of the Bardeen–Kiselev quasinormal frequencies and compare these frequencies with the RN-dSQ. Interestingly, it shows that the responses of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ under electromagnetic perturbations are different when the charge parameter q, the state parameter w and the normalization factor c are varied; however, for the gravitational perturbations, the responses of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ are different only when the charge parameter q is varied. Therefore, we may deduce that compared with the gravitational perturbation, the electromagnetic perturbation can be used to understand nonlinear and linear electromagnetic fields in curved spacetime separately.

Finally, with regard to the QNM's oscillation shape, we use the finite difference method to study the dynamical evolution of the nonlinear electrodynamics field perturbations and gravitational perturbations in the time domain [60, 61] and examine the stability of the Bardeen–Kiselev BH with the cosmological constant. Figure 5 shows the evolution of perturbation in log scale. We can see that the QNM oscillation of electromagnetic perturbations in the Bardeen–Kiselev BH with the cosmological constant is sensitive to the smaller normalization factor c. Also, when the normalization factor c takes smaller values, the QNM oscillation of gravitational perturbations in the Bardeen–Kiselev BH with the cosmological constant and the absence of Kiselev quintessence are almost the same; therefore, they are not sensitive to the smaller c values. Meanwhile, for the gravitational perturbations, as the normalization factor c increases, the oscillation frequency will decrease; this result is consistent with the result obtained by the WKB method. Moreover, in view of the finite value of the cosmological horizon r_c, our integration domain is limited in the range of −r_* to r_*; it is tricky to obtain numerical values of φ at a late time. Thus, there is no power-law tail in the dynamics and, if we observe equation (1), it is found that as λ → 0, r_c → ∞; therefore, numerically, the domain of integration also becomes large enough and one can finally obtain very late time dynamics, but based on the fact that this small value of λ in numerical computation has its own challenges.

Hawking [62] predicted that when the temperature of the BH is proportional to its surface gravity, the BH will emit particles, which behave almost like a blackbody. Basically, BHs can be interpreted as the thermal systems that have an associated temperature and entropy; therefore, BHs produce radiation, and such radiation is called Hawking radiation [63]. In the literature, the WKB method has been used to obtain reflection and transmission coefficients (greybody factors) [64–68], such as in the context of braneworld models and wormholes. In this section, we intend to discuss the frequency-dependent reflection R(ω) and transmission coefficient T(ω) for a scattering process of the gravitational and electromagnetic perturbations from the Bardeen–Kiselev BH with the cosmological constant. According to the Hawking radiation, at the event horizon, the emission rate of a BH in a mode with frequency ω iswhere α represents the inverse of the Hawking temperature and the symbol ± corresponds to fermions (bosons). However, not all the radiation is able to reach the distant observer, i.e. a part of the radiation will be tunneled through the potential barrier and will reach the distant observer, while the other part will be reflected back toward the BH, and the radiation recorded by the distant observer will no longer appear as a blackbody. Therefore, the emission rate measured by an observer at infinity for a frequency mode ω can be expressed aswhere γ_l is the greybody factor, which is defined asThe formalism of the asymptotic behavior of the wave after scattering off of the effective potential can be expressed in the tortoise coordinate as In equations (32) and (33), the reflection and transmission coefficients are functions of the oscillation frequency ω of the wave. The reflection coefficient in the presence of the WKB approximation is given bywhere η isNote that the values of Λ_i can be found from equation (28); thus by conserving probability, we have

Next, using the WKB method, we will describe the calculation of R(ω) and T(ω). Besides, if r₀ is the value of r, where the the potential V(r) is the maximum, as seen in [69], there are three cases to explain the relation between ω² and V(r₀), and we we will focus on ω² ∼V(r₀), because the WKB approximation has high accuracy.

By using the effective potentials from section 3, we numerically plot the variation of the reflection and transmission coefficients on the frequency of electromagnetic and gravitational perturbations with various parameters (ℓ, q, c, λ, w) in figures 6 to 9. In these figures, we can see that the value of the greybody bound is zero when the frequency is minimal, and the value of the greybody bound turns out to be 1 when the frequency is large enough, which implies that the wave is basically totally reflected when the frequency is small. Meanwhile, in view of the tunneling effect, a partial wave could pass through the potential barrier when the frequency increases, or the wave will not be reflected when the frequency reaches a certain level. Moreover, note that according to equations (31) and (36), for the high-frequency and low-frequency limits, the properties of the reflection and transmission coefficients should be opposite; in fact, these results can be verified from figures 6 to 9.

The electromagnetic perturbations (s = 1) and gravitational perturbations (s = 2) in figure 6 exhibit similar behavior, i.e. for a fixed frequency, the greybody factor will decrease with the increasing multipole number ℓ; the response of the greybody factor for different ℓ under gravitational perturbation is larger than the case under electromagnetic perturbation.

Similar electromagnetic and gravitational perturbation behavior appears in figure 7. For a fixed frequency, the greybody factor will decrease with the increasing magnetic charge q; this result is different to [32], and it may be the reason that the existence of Kiselev quintessence changes the metric (equation (1)). Besides, for the same spacetime parameters, the response of the greybody factor for different q under electromagnetic perturbations is larger than the case under gravitational perturbations.

Figure 8 shows that under electromagnetic and gravitational perturbations, for a fixed frequency, the greybody factor will increase with the increasing normalization factor c, and the response of the greybody factor for different c under electromagnetic perturbations is larger than the case under gravitational perturbations.

Figure 9 shows similar behavior of electromagnetic and gravitational perturbations, i.e. for a fixed frequency, the greybody factor will slightly increase with the increasing cosmological constant λ.

In these figures, for the same parameters in terms of electromagnetic and gravitational perturbations, the greybody factor will decrease with the increasing state parameter w and, considering the effect of the normalization factor c on the greybody factor in figure 8, we may deduce that due to the presence of Kiselev quintessence, the transmission coefficient will increase. As mentioned above, the properties of the reflection and transmission coefficients are opposite; thus, we will not present a specific description of reflection coefficients.

Finally, we will study the total absorption cross section in the context of electromagnetic and gravitational perturbations for different parameter spaces in the Bardeen–Kiselev BH with the cosmological constant background. The total absorption cross section is given by

Figure 10 shows the total absorption cross section under electromagnetic (s = 1) and gravitational (s = 2) perturbations in the Bardeen–Kiselev BH with the cosmological constant. For convenience, we have summed over ℓ = 2 to ℓ = 8 modes to determine σ. It shows that for a fixed frequency, the absorption cross section is always larger for gravitational perturbations with respect to electromagnetic ones, and this result is consistent with [32]. Moreover, as the transmission coefficient approaches 1 at some critical value of ω, whether it is with regard to the electromagnetic perturbations or gravitational perturbations and regardless of the BH parameters, the total absorption cross section falls off as $\tfrac{1}{{\omega }^{2}}$; therefore, we can find the fall-off region in this figure.

In this work, we studied the gravitational and electromagnetic perturbations in the Bardeen–Kiselev BH with the cosmological constant in terms of QNMs and compared it with the RN-dSQ. We present the effective potentials under the two perturbations. It is shown that both potentials are positive definite between the event and cosmological horizons, and they all have a single maxima, which indicates that for a fixed set of parameters, the potentials decrease with increasing parameters ∣w∣ or λ; basically, the smaller value of ∣w∣ or λ suppresses the emission modes for gravitational and electromagnetic perturbation. Next, using the sixth-order WKB method, we obtain the quasinormal frequencies of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ under gravitational perturbations in figures 2, 3 and 4. We show that the response of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ in terms of the imaginary part of ω are different only when the charge parameter q is varied. Besides, in figure 3, we can see that the absolute values of the imaginary parts as well as the real parts of the quasinormal frequencies of the Bardeen–Kiselev BH with the cosmological constant with Kiselev quintessence (c ≠ 0) are smaller compared to those without Kiselev quintessence (c = 0) and, as the parameter c increases, they will decrease. Also, figures 2 to 4 show that the oscillation frequency and the damping rate of the Bardeen–Kiselev BH with the cosmological constant increase with the increasing w. Meanwhile, increasing c or decreasing w implies increasing the density of Kiselev quintessence. Therefore, we can remark that due to the presence of Kiselev quintessence, the gravitational perturbations of the Bardeen–Kiselev BH with the cosmological constant damp and oscillate more slowly; this behavior increases when increasing the density of Kiselev quintessence, which implies that Kiselev quintessence reduces the dissipative effect of the BH on its neighborhood.

In tables 1 to 3, we calculated the quasinormal frequencies in electromagnetic perturbations of the Bardeen–Kiselev BH with the cosmological constant and compared these frequencies with the RN-dSQ. We showed that due to the presence of Kiselev quintessence, the electromagnetic perturbations around the Bardeen–Kiselev BH with the cosmological constant damps faster and oscillates slowly, and this behavior is different to the gravitational perturbations. Therefore, from these tables, the responses of the Bardeen–Kiselev BH with the cosmological constant and the RN-dSQ with the cosmological constant under electromagnetic perturbations are different when the charge parameter q, the state parameter w and the normalization factor c are varied.

Next, with regard to the QNM oscillation shape, we used the finite difference method to study the dynamical evolution of the nonlinear electrodynamics field perturbation and gravitational perturbation in the time domain.

Moreover, we studied the greybody factors of the Bardeen–Kiselev BH with the cosmological constant under gravitational and electromagnetic perturbations and observed that:Finally, we investigated the total absorption cross section under electromagnetic and gravitational perturbations in the Bardeen–Kiselev BH with the cosmological constant. It is shown that for a fixed frequency, the absorption cross section is always larger for gravitational perturbations with respect to electromagnetic ones.