Entanglement entropy and Page curve of black holes with island in massive gravity

By applying the island rule proposed recently, we compute the entanglement entropy of Hawking radiation and study the Page curve for the eternal black holes in massive gravity. We investigate for both the neutral and charged black holes which the corresponding results of Schwarzschild and Reissner–Nordström black holes are restored in the limit of massless graviton. We show for the neutral and non-extremal charged black holes that the island is not formed at the early times of the evaporation and hence the entanglement entropy increases linearly in time. However, for the extremal charged black hole, the calculation of the entanglement entropy at the early times without the island is ill-defined because the metric is divergent at the curvature singularity. This implies that new physics in the UV region must be taken into account to make the metric behaving smoothly at the very short distances. At the late times, with the emergence of one island near the event horizon, the entanglement entropy is saturated by the Bekenstein–Hawking entropy of black holes. In addition, we analyze the impact of massive gravity parameters on the size of island, the entanglement entropy, the Page time, and the scrambling time in detail.


Introduction
The behavior of back holes as the thermodynamic objects with the temperature and entropy (which are determined in terms of the surface gravity and the horizon area, respectively) provides a deep connection between the research areas of general relativity, thermodynamics, and quantum mechanics [1][2][3][4][5]. In addition, the black hole thermodynamics may offer indispensable insights into quantum gravity. By taking into account the quantum effects of the matter fields near the event horizon, S. Hawking showed that black holes can emit the radiation with the nearly thermal spectrum [5]. Suppose that black holes are formed by the gravitational collapse of a e-mail: nam.caohoang@phenikaa-uni.edu.vn (corresponding author) the matter with the initial state in a pure (quantum) state, after black holes evaporate completely the final state of the system would be in a mixed (thermal) state. In this way, the quantum information of the initial state is not conserved in the evaporation process of black holes, leading to the so-called information loss paradox [6]. On the other hand, the time evolution in this way is contradictory to one of the fundamental principles of quantum mechanics, namely the unitarity principle which requires that the final state must be the pure state if the system starts from the same kind of the state.
If a black hole is formed from a pure state, then its entropy is zero. Later, the black hole would evaporate due to emitting the Hawking radiation. During the early time of the black hole evaporation, the entanglement entropy of Hawking radiation should increase in time because more and more Hawking quanta are emitted and entangled with the remaining black hole. Whereas, the thermodynamic or Bekenstein-Hawking entropy of the black hole would decrease due to its horizon area shrinking. The behavior of the entanglement entropy of Hawking radiation changes at the time when the Bekenstein bound is violated. The Bekenstein bound implies that the fine-grained entropy of the black hole should not be larger than the Bekenstein-Hawking entropy of the black hole. As the degrees of freedom of the remaining black hole together with the outgoing radiation as a whole is pure state, the finegrained entropy of the remaining black hole should be equal to the entanglement entropy of the Hawking radiation. Hence, at the moment at which the entanglement entropy of Hawking radiation is equal to the Bekenstein-Hawking entropy of the remaining black hole, called the Page time, the entanglement entropy of Hawking radiation must decrease and drop down to zero as the black hole evaporates completely. This means that the quantum information of the initial black hole is encoded in the Hawking radiation and the black hole evaporation is consistent with the unitarity evolution. In this way, the unitarity evolution of the black hole evaporation corresponds to that the entanglement entropy of Hawking radia- tion follows the Page curve [7,8]: the entanglement entropy of Hawking radiation increases monotonically at the early period of the black hole evaporation, then reaches a maximum value at the Page time, and finally has to go to zero at the end of the evaporation process, as indicated in Fig. 1. It is expected that the Page curve is obtained in a theory of quantum gravity. Therefore, reproducing the Page curve for the time evolution of the black hole evaporation is an important step towards not only resolving the information loss paradox of black holes but also understanding completely quantum gravity.
The recent works have shown a significant progress in deducing the Page curve for the entanglement entropy of Hawking radiation by taking into account the configuration with the islands [9][10][11][12]. The islands are some regions I which are completely disconnected from the region R of Hawking radiation which is assumed to be far away from black holes such that the backreaction of Hawking quanta on the spacetime geometry is negligible. The boundaries of the islands extremize the generalized entropy functional and hence they are called the extremal surfaces. A density matrix relating to the states of Hawking radiation is normally calculated by taking the partial trace over the states in the complementary part of the radiation region. However, by using the quantum extremal surface technique, it was found that the islands appear in the complementary region of R at the late times of the black hole evaporation process and hence the states in the islands should be eliminated from the states which are traced out. According to the quantum extremal surface prescription, the entanglement entropy of Hawking radiation is obtained as the minimum value of the generalized entropy functional [9,[13][14][15][16][17] S(R) = min ext A(∂ I ) where G N is the Newton gravitational constant, ∂ I refers to the island boundaries, A(∂ I ) is the total area or volume of The island formula (1) means that the entanglement entropy of Hawking radiation is obtained as the minimum value of the generalized entropy functional over all possible extremal surfaces corresponding to all possible locations of the island. The entanglement entropy of Hawking radiation calculated in this method contains two contributions which come from the area term of the island and the von Neumann entropy of the quantum fields on the union of the radiation and island regions. At the early time of the black hole evaporation, there is no island forming and since the black hole entanglement wedge contains all the black hole interior. The entanglement entropy of Hawking radiation thus is the von Neumann entropy of quantum fields in the black hole exterior and it increases as more and more Hawking quanta are emitted. However, at the late time, there is the emergence of the island whose boundary is very close the black hole horizon and which extends almost through the whole black hole interior. Consequently, the partners of Hawking quanta which fall inside the black hole are almost contained in the island. This means that the von Neumann entropy of the quantum fields on the union of the radiation and island regions is small. On the other hand, the dominant contribution for the entanglement entropy of Hawking radiation is from the area term of the island. When the black hole shrinks due to the evaporation, this term would decrease. Therefore, the entanglement entropy of Hawking radiation would go to zero when the black hole evaporates completely. This means that the behavior of the entanglement entropy of Hawking radiation with the island method follows the Page curve, as depicted in Fig. 2.
Interestingly, it was shown that the configuration with the island can be emerged from the gravitational Euclidean path integral using the replica trick [18,19]. In the gravitational replica method, the different replica sheets with the boundary conditions kept fixed are connected together by the so-called replica wormholes which are new saddle points in the gravitational Euclidean path integral. In the semi-classical limit, they are the dominant contributions to the effective action for the geometry which is a sum of the gravitational action and the partition function of the quantum fields. In this way, the presence of the replica wormholes yields the same island formula for the entanglement entropy of Hawking radiation as using the quantum extremal surface technique.
The island rule was initially considered for the twodimensional gravitational systems where the explicit computations for the entanglement entropy of Hawking radiation and the Page curve can be easily performed by using the semi-classical method due to the presence of the high symmetries. For the two-dimensional black holes in the context of Jackiw-Teitelboim (JT) gravity, the islands are emerged at the later times of the black hole evaporation and hence their presence makes the entanglement entropy of Hawking radiation remaining finite at the late stage of the evaporation process [12,18]. The island consideration in the twodimensional models was extended to study the asymptotically flat 2d dilaton black holes [20], the two-sided Janus black holes [21], and the evaporating black holes [22][23][24].
Although the four-and higher-dimensional gravitational systems are complicated due to the lack of the symmetric analyses, the recent works have demonstrated that the island rule is applicable to calculate the entanglement entropy of Hawking radiation and reproduce the Page curve consistent with the unitarity time evolution. For the four-and higherdimensional eternal Schwarzschild black holes, the authors in [25] showed the emergence of an island whose dominant contribution leads to the finiteness of the entanglement entropy of Hawking and the Page curve incorporating the unitarity principle. Here, the boundaries of the island are located in the outer vicinity of the event horizon of the Schwarzschild black holes. For considering the backreaction of Hawking radiation, which is natural in the context of black holes in a cosmology supported by the radiation, the existence of islands was pointed out in the cosmological braneworlds [26]. These results have also been confirmed for the Reissner-Nordström (RN) black holes [27,28], the charged/neutral dilaton black holes [29][30][31], the Kaluza-Klein black holes [32], and the black holes including the higher derivative terms [33]. Additionally, the islands corresponding to the left/right entanglement of a conformal defect was studied in Randall-Sundrum braneworld model involving weakly gravitating bath [34]. In this direction, there are also the investigations about the entanglement of purification and complexities for multi-boundary wormhole models of islands [35,36].
As mentioned, the island rule has been extended to calculate the entanglement entropy of Hawking radiation for various black hole geometries which are well-known in the literature for the four-and higher-dimensional cases. But, graviton is massless in the gravity frameworks which are equipped to derive these black hole geometries. In addition, it was argued that the island proposal coupling the gravitating bath induces a mass for the bulk graviton [37]. Whereas, the authors in [38] argued that the islands might not constitute the consistent entanglement wedges in the gravity theories with massless graviton. These results imply that calculating the entanglement entropy of Hawking radiation and the Page curve for black holes using the island method would be in the context of massive gravity.
Considering a nonzero mass of graviton is one of the infrared (IR) modifications of gravity, which has the cosmological consequences of which the lately accelerating expansion of the universe can be naturally explained without invoking dark energy. 1 There are the first attempts to construct a gravity theory accompanied with a nonzero mass of graviton, such as Fierz-Pauli (linear) massive gravity [40] or nonlinear massive gravity with the Vainshtein mechanism [41]. However, these massive gravity theories suffer from the pathological problems which are well-known as the van Dam-Veltman-Zakharov discontinuity [42][43][44] (the predictions in the massless limit do not coincide with those of Einstein gravity) and the Boulware-Deser (BD) ghost [45]. These pathological problems thus became the obstacle in establishing a consistent theory of massive gravity. Until recent years, a significant progress has been made by de Rham, Gabadadze and Tolley (dRGT) who proposed a nonlinear massive gravity which is free of the BD ghost and gives the same prediction in the massless limit as Einstein gravity does [46,47]. It was shown that the ghost-free potential structure of massive gravity does not change under the quantum corrections at one-loop [48]. In the metric formulation of massive gravity, it is difficult to write down the interactions of different helicity modes of massive spin-2 field due to the square root structure of potential. However, using the vierbein language one can obtain the decoupling limit of massive gravity which consists of the full interactions of the helicity-1 and helicity-0 modes [49]. For reviews about the recent progress in massive gravity, see [50,51].
It was also pointed that massive gravity has the interesting implications for cosmology. The authors in [52] argued that massive gravity allows open Friedmann-Robertson-Walker (FRW) cosmological solutions contrary to the no-go result [53] which does not extend to the open FRW universes. Perturbations around FRW background solutions were studied with the extra metric which is promoted to be dynamical [54]. New massive gravity theories proposed in [55] contain stable cosmological solutions without the requirement of infinite strong coupling. The approaches towards healthy cos-mological solutions were proposed in [56]. In addition, various black hole solutions and their thermodynamic properties have been extensively investigated in massive gravity [57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74]. The novel thermodynamic phenomena of black holes in the extended phase space with the cosmological constant identified as the thermodynamic pressure have been pointed out in the context of massive gravity like van der Waals phase transition of black holes [75], the efficiency of heat engine provided by black holes [76,77], and Joule-Thomson expansion of black holes [78]. The entanglement entropy in the van der Waals phase transition in the context of massive gravity was studied in [79]. For another massive object, the the structure of neutron star was investigated in the context of massive gravity [80]. Furthermore, the effects of massive gravity on s-/p-wave holographic superconductors have been explored in [81][82][83].
Because of many interesting aspects of massive gravity, it is worth to study the entanglement entropy of Hawking radiation and the Page curve by applying the island rule for the black hole geometries in the context of massive gravity and analyze the effects due to the nonzero mass of graviton on them. We consider the four-dimensional black hole solutions in massive gravity which were found in [58]. The presence of graviton mass in massive gravity leads to two new terms in the black hole solutions. The first term corresponds to a linear term in the radial coordinate r which would play the role of the quintessence matter with the proper coupling parameter of massive gravity. Whereas, the second term corresponds to the global monopole solution which was obtained in the presence of topological defect as a result of the spontaneous symmetry breaking in the early universe [84]. Interestingly, the presence of the linear term leads to the asymptotically non-flat behavior of black hole geometries. The entanglement entropy of Hawking radiation and the Page curve are explored by the island method for various black hole geometries corresponding to the asymptotically flat (AdS/dS) behavior. Hence, it is worth checking whether the island method is still able to be applied to the black hole geometries with the asymptotically non-flat behavior, which is important to show the wide range of applicability of the island method. In addition, the linear term and global monopole term modify the dynamic properties and the evaporation of black holes. Thus, these terms would lead to the changes in the entanglement entropy of Hawking radiation and the Page curve. We find the expressions for the entanglement entropy of Hawking radiation, the Page time, and the scrambling time in the presence of linear term and global monopole term and analyze their impacts on these quantities. We point to that in the limit of massless graviton the expressions for the entanglement entropy of Hawking radiation, the Page time, and the scrambling time would reduce to the corresponding expressions in the context of Einstein gravity. This shows a continuity of dRGT massive gravity with Einstein gravity in the limit which the graviton mass approaches zero and thus supports theory of dRGT massive gravity as the consistent theory of modified gravity. This paper is organized as follows. In Sect. 2, we introduce the four-dimensional black hole solutions in massive gravity which were found in [58]. Then, we determine various related coordinates and rewrite the metric in the Kruskal coordinate, which are all basis for the evaluations in the next sections. In Sect. 3, we calculate the entanglement entropy of Hawking radiation for the neutral black hole without and with the island which corresponds to the early and late times of the black hole evaporation, respectively. In addition, we derive the Page and scrambling times and study the effects of massive gravity on these quantities, the location of the island boundaries, and the entanglement entropy of Hawking radiation. In Sect. 4, we repeat the calculations of Sect. 3 for the charged black hole in massive gravity for the non-extremal and extremal cases. With respect to the non-extremal case, the effects of massive gravity on the Page time and the scrambling time are qualitatively the same as the neutral black hole when the distance between the event and Cauchy horizons is not so close, but in the contrary they drastically change. The results found for the extremal case are almost different from the one of the neutral and non-extremal charged black holes due to their causal structure. Finally, we conclude in Sect. 5.

Black holes in massive gravity
In this section, we review briefly the black hole solutions in the framework of massive gravity in four dimensions, found in [58]. Then, we introduce the Kruskal coordinate and rewrite the line element in this coordinate for various black hole geometries in massive gravity.
The action of massive gravity coupled to the Maxwell field in four dimensions is given by where R is the scalar curvature of spacetime, m g is the graviton mass, c i are the coupling parameters, f is the reference (fiducial) metric which is not dynamical, U i are symmetric polynomials which are written in terms of the eigenvalues of the 4 × 4 matrix K μ ν = g μλ f λν as , The corresponding equations of motion read where With the choice of gauge-fixed ansatz for the reference metric as where c 0 is a constant set to be one without loss of generality, the spherically symmetric black hole solution is found as [58] ds where with M and Q to be the ADM mass and the electric charge of the black hole. We see that this black hole solution would reduce to the Schwarzschild black hole in the massless limit m g → 0. The asymptotic behavior of the black hole solution is dependent on the sign of the coupling parameter c 1 . For c 1 < 0, the term m 2 g c 1 r/2 in the metric function f (r ) plays the role of the quintessence matter with the quintessential state parameter ω q = −2/3 [85,86]. In this situation, besides the solutions of the black hole horizons, the equation f (r ) = 0 leads to the solution of the cosmological horizon. In this work, we only consider the positive sign of c 1 which corresponds to the absence of the cosmological horizon.
In the following, we introduce the Kruskal coordinate and rewrite the line element in this coordinate for two cases: (i) Neutral black hole and (ii) Charged black hole.

Neutral black hole
As the electric charge of the black hole is zero, i.e. Q = 0, we obtain the neutral black hole solution. In this case, the black hole only has one (event) horizon which is denoted by r h given by We can rewrite the metric function f (r ) in the event horizon r h as follows where By defining the tortoise coordinate as where κ refers to the surface gravity at the event horizon given by we introduce the Eddington-Finkelstein coordinate constructed relying on the paths of the radially incoming and outgoing photons as follows Then, the Kruskal coordinate is defined as The line element is rewritten in terms of the Kruskal coordinate as where the conformal factor W 2 (r ) is given by

Charged black hole
From the behavior of the black hole mass function in terms of the horizon radius as depicted in Fig. 3, we see that if the black hole mass is larger than a minimum value M min , the black hole would possess two different horizons which are the event (outer) horizon r + and the Cauchy (inner) horizon r − , corresponding to the non-extremal case. In this situation, by parameterizing the metric function f (r ) given in Eq. (7) in terms of the horizons r ± , f (r ) is rewritten as where the ADM mass M and the electric charge Q of the black hole are expressed in terms of the horizons r ± as The tortoise coordinate is defined as where κ ± are the surface gravity at the horizons r ± and are given by In the Kruskal coordinate defined as the line element of the non-extremal charged black hole is rewritten as where the conformal factor W 2 (r ) reads Now we consider the case that the black hole mass is equal to the minimum value M min , corresponding to the extremal case. In this situation, two horizons r ± coincide together and hence the black hole possesses only one horizon at The metric function f (r ) in Eq. (17) thus becomes which corresponds to the tortoise coordinate as Then, in terms of the Kruskal coordinate defined as we rewrite the line element of the extremal charged black hole as where the conformal factor W 2 (r ) is given by

Entanglement entropy for neutral black hole
In this section, we are interested in evaluating the entanglement entropy of Hawking radiation which is emitted by the black hole and identified as the matter sector coupled to gravity. We consider the contribution to the entanglement entropy from the configurations without and with the islands. For the case of the island configuration, we will restrict to the calculation in the presence of one island for simplification without loss of generality. We will indicate that the configuration of one island is sufficient to reproduce the finiteness for the entanglement entropy at the late times and thus it would lead to the Page curve consistent with the unitarity principle. represent the boundaries (cutoff surfaces) of R − and R + , respectively. I refers to the island, extending between the left and right wedges, whose boundaries are denoted by a ± The Penrose diagram of the eternal neutral black hole (which is in the thermal equilibrium with a heat bath) without the island in massive gravity is depicted in the left panel of Fig. 4 corresponding to the maximally extended spacetime The region of the Hawking radiation is given by the union of two regions R − and R + which are located in the left and right wedges of the Penrose diagram, respectively. The boundaries of the regions R − and R + , which are introduced to avoid the IR divergence, are denoted by b − and b + , respectively. The (t, r ) coordinates where β is the inverse of the Hawking temperature of the black hole. For the region of Hawking radiation which is sufficiently far away from the black hole, we can ignore the backreaction of the matter sector on the spacetime geometry. With the distance between two boundaries of b + and b − to be large enough compared to the scale of the size of these boundaries, we can consider the approximation of the twodimensional conformal field theory (CFT) to compute the entanglement entropy [25]. We also assume that the initial state of the system is in the pure state and hence the entanglement entropy in the radiation region R − ∪ R + is equal to that in its complement which is one interval [b − , b + ]. In the absence of the island, the entanglement entropy is given by [10,87] where c is the central charge of the two-dimensional CFT and d(b − , b + ) is the geodesic distance between b − and b + which is calculated for the metric of the following form ds 2 = −W 2 dU dV .
In the presence of an island, the corresponding Penrose diagram is shown in the right panel of Fig. 4. The boundaries of the island located in the left and right wedges of the Penrose diagram are represented by a + and a − whose (t, r ) coordinates are (t a , a) and (−t a + iβ/2, a), respectively. In the island construction, the generalized entropy is given as where the first term is the Bekenstein-Hawking entropy coming from the contribution of two island boundaries and the entropy of the matter sector S mat (R − ∪ R + ∪ I ) is calculated by the formula for two intervals as [88] According to the prescription of the quantum extremal surface, the dominant contribution for the entanglement entropy comes from the configuration which minimizes the generalized entropy. Therefore, we shall compute the entanglement entropy by extremizing the generalized entropy over all possible boundary surfaces of the island and then take the minimal value.

Entanglement entropy without island
First we study the behavior of the entanglement entropy of Hawking radiation in the case of the no island configuration. From Eq. (30), one obtains an explicit expression of the entanglement entropy for the neutral black hole in massive gravity in the absence of the island as where we have used r h b in third line due to the boundaries of the radiation region assumed to be far away from the black hole. At the early time approximation t b 1/κ or t b r h , we have Here, the first term is the entanglement entropy of the radiation region at the initial state up to a constant and the second term manifests the quadratic growth of the entanglement entropy with time at the early times. For t b 1/κ corresponding to the late times, we find This expression implies that the entanglement entropy of Hawking radiation increases linearly in time and becomes infinite as t b → ∞, which corresponds to the fact that the distance between two regions R − and R + is very large at the late times. The infinitely large value of the entanglement entropy in the limit t b → ∞ means that the information about the initial-state matter which collapses into the black hole or the information about the particles falling into the black hole cannot be retrieved from the Hawking radiation.
In this sense, the contribution of the no island configuration to the entanglement entropy leads to the non-unitary time evolution of the black hole in the evaporation process. In the next subsection, we will show that this conflicting issue with the unitarity of quantum mechanics can be resolved with the island configuration which emerges at the late times of the evaporation process.

Entanglement entropy with an island
We calculate the entanglement entropy of Hawking radiation with including the contribution of the configuration with an island. From Eqs. (31) and (32), we can write the generalized entropy in the presence of one island for the metric (15) as where the first term comes from the two-sided area of the boundaries of the island and the second and third terms are the contributions of the matter fields on the union of the radiation and island regions.
In the presence of the island, let us study the behavior of the entanglement entropy at the early and late times of the black hole evaporation: Early times At the early times of the black hole evaporation, the entanglement entropy of Hawking radiation is small and hence the island should be lie inside the black hole. In addition, we assumed the boundaries of the radiation region far away from the event horizon of the black hole. These lead to the following approximation [25,27] As a result, we can properly neglect the third term in the expression (36) and thus we obtain where the ellipses refer to the terms which are independent on a and t a . In order to find the position of the island boundaries which extremize the generalized entropy S gen , we need to extremize S gen given in Eq. (38) over all possible (a, t a ) locations of island according to the island method [9,[13][14][15][16][17], which means that we need to solve the following extremizing equations Note that, on the black hole interior the radial coordinate r is actually the timelike coordinate. Hence, at the early times we have a/r h 1. Using this approximation and cG N /r 2 h 1, we obtain the location of the island boundaries as where l P ≡ √ G N is the Planck length. This result indicates that the size of the island at the early times is about the Planck length. However, the upper cutoff length in the derivation of the island formula should be far above the Planck length where the Planck scale physical degrees of freedom are integrated out. This can be realized from the assumption of the replica symmetry, which is broken by the effects of quantum gravity, to derive the island formula for the entanglement entropy [18,19]. In this sense, the island would not emerge at the early times. Thus, the entanglement entropy at the early times should be determined by the geometry configuration without island and since it grows with time, as analyzed in the previous subsection.
Furthermore, we note that massive gravity does not make a significant effect on the size of the island given in Eq. (40). This is realized from the fact that the nonzero mass of graviton can be considered as the IR or long-length corrections which hence do not affect significantly on the physics in the Planck scale region.
Late times We turn to consider the late time behavior of the entanglement entropy. In this situation, the approximation is taken as [25,27] which leads to With this approximation, the time-dependent component of the generalized entropy is approximated as From this expression, it is straightforward to find that extremizing the generalized entropy with respect to t a leads to t a = t b . Substituting this result into the generalized entropy within the approximation (41), we obtain the following approximate expression This expression suggests that the entanglement entropy is no longer dependent on time at the late times. We consider the situation that the island is located near the event horizon, i.e. a = r h + +O( 2 ) with 1. By solving perturbatively the extremizing condition ∂ S gen /∂a = 0, we find Clearly, the second term is positive and is suppressed by the second power of cG N /r 2 h , which is really small as expected. This implies that the location of the island boundaries is slightly outside the event horizon of the black hole. We also observe that the nonzero mass of graviton affects directly and indirectly on the the location of the island boundaries through the terms relating to δ and the modification of the event horizon radius r h , respectively. However, due to the fact that the second term in Eq. (45) is small, the direct effects of massive gravity on the location of the island boundaries are negligible compared to its indirect effects. From the expansion of the event horizon radius in terms of the graviton mass m g for the small m g situation (which is consistent with the fact) as we realize that the nonzero mass of graviton which makes the location of the island boundaries shifting inside or outside the black hole depends on the sign of the term c 2 + c 1 M. If the black hole mass is larger than the ratio −c 2 /c 1 or the event horizon radius of the corresponding Schwarzschild black hole (in Einstein gravity) is larger than the ratio −2c 2 /c 1 , the location of the island boundaries would be closer to the event horizon of the black hole compared to Einstein gravity. On the contrary, the location of the island boundaries would be pushed further.
In addition, as the location of the cutoff surface approaches the event horizon, the second term in Eq. (45) becomes large and hence it is no longer considered as the correction. This means that in this situation the boundaries of the island are not close to the event horizon. Consequently, some part of the island would lie inside the region of Hawking radiation and hence the concept of the island here does not make sense. On the other hand, the calculation of the entanglement entropy using the island formula is invalid if the radiation region is close to the event horizon.
By substituting the location of the island boundaries given in Eq. (45) into the approximate expression of the generalized entropy at the late times given in Eq. (44), we determine the entanglement entropy as Here, the first term is twice of the Bekenstein-Hawking entropy of the black hole, which comes from the area of two-sided island and it is the dominant contribution to the entanglement entropy. On the other hand, the entanglement entropy which is observed in a single side of the Penrose diagram is approximately the Bekenstein-Hawking entropy. The second term is the logarithm correction for the entanglement entropy, which comes from the quantum nature of the matter fields. Other correction terms are strongly suppressed by the powers of cG N /r 2 h and thus they are very small. In this way, the presence of the island which becomes the dominant configuration at the late times leads a finite saturation value for the entanglement entropy which is approximated at the leading order to be twice of the Bekenstein-Hawking entropy. The expansion (46) implies that at the leading order the nonzero mass of graviton reduces(raises) the bound of the entanglement entropy of Hawking radiation compared to the massless case if the black hole mass is larger(smaller) than the ratio −c 2 /c 1 .
In summary, we have calculated the entanglement entropy of Hawking radiation by extremizing the generalized entropy over all possible boundary surfaces of the island from which we determine the location of the island boundaries and then we obtain a minimum value. At the early stage of the black hole evaporation, no island is emerged and the configuration without the island minimizes the generalized entropy. As a result, in during this moment the entanglement entropy is approximately a linearly increasing function in time. However, this behavior of the entanglement entropy changes drastically at the late stage of the black hole evaporation as the island appears with its boundaries located slightly outside the event horizon. The configuration with the island leads to a minimum value of the generalized entropy which is twice of the Bekenstein-Hawking entropy of the black hole at the leading order. Interestingly, transition between the linear growth and time-independent constant behaviors of the entanglement entropy can be realized from Page's argument that the entanglement entropy can be approximately given by the thermal entropy of the subsystem if it is sufficiently small compared to the total system [7,8]. In the beginning of the black hole evaporation where the amount of Hawking radiation emitted by the black hole and entering into the boundaries of the radiation region is still small, the radiation region is substantially smaller than the total system. Thus, the entanglement entropy can be approximated by the thermal entropy of the radiation which grows with time due to the increasing of the amount of the radiation. But, as the amount of the radiation becomes more and more at the late stage of the black hole evaporation, the subsystem would be replaced by the black hole and hence the entanglement entropy stops the growth with time and reaches the bound of the Bekenstein-Hawking entropy of the black hole.
We can check that our calculations for the location of the island boundaries and the entanglement entropy with the island configuration would reduce to the corresponding results of the four-dimensional Schwarzschild black hole [25] in the limit of that the graviton mass goes to zero, i.e. m g → 0. By using the fact that m g → 0 corresponds to δ → ∞ and the following limit we find the results in the case of the four-dimensional Schwarzschild black hole as follows where r h = 2M.

Page time and scrambling time
The results of the previous subsections allow us to sketch the behavior of the entanglement entropy in time as shown in Fig. 5. We observe that, at the first stage of the black hole evaporation, the entanglement entropy grows linearly with time due to the dominant contribution of the configuration without the island. At the late stage of the black hole evaporation, the island emerges near the event horizon and becomes the preferred configuration. As a result, the linear growth of the entanglement entropy reaches maximum and is replaced by a constant. The Page time is the moment that the entanglement entropy reaches the maximal value corresponding to the transition between two configurations without and with the island. From Fig. 5, we see that the Page time can be approximately calculated from the crossing of the solid blue line (the early times without the island) and the solid red line (the late times with the island). On the other hand, by equating Eqs. (35) and (47), we derive the Page time as where the first term is the Page time for the Schwarzschild black hole in the context of Einstein gravity (corresponding to the massless graviton) and the second and higher order terms are the contributions due to the nonzero mass of graviton. The expression (51) suggests that the nonzero mass of graviton makes the evaporation of the neutral black hole more quickly to reach the Page time compared to the massless case if the black hole mass M is larger than the ratio −4c 2 /5c 1 . Of course, this happens the contrary for M < −4c 2 /5c 1 . The reduction (or the increasing) of the Page time in massive gravity with the proper coupling parameters c 1,2 can be understood as follows. First, with M > −c 2 /c 1 the Bekenstein-Hawking entropy of the black hole or the bound of the entanglement entropy in massive gravity is smaller than that in Einstein gravity. This means that the system in massive gravity takes a shorter duration to reach the entropy bound if the Hawking temperature of the black hole in two theories of gravity is the same. Second, by expanding the Hawking temperature of the black hole in massive gravity as where the first term is the Hawking temperature of the conventional Schwarzschild black hole and the remaining terms come from the presence of the graviton mass, we find that the Hawking temperature of the neutral black hole in massive gravity is larger than that in Einstein gravity in the region of M > −2c 2 /3c 1 . The larger temperature means that the black hole emits the Hawking radiation more rapidly and hence it yields the appearance of the island more early. Combining the behavior of both the Bekenstein-Hawking entropy and the Hawking temperature in terms of the coupling parameters of massive gravity leads the intermediate value −4c 2 /5c 1 which lies between −c 2 /c 1 and 2c 2 /3c 1 . The presence of the island reproduces the behavior of the entanglement entropy following the Page curve. Thus, we can consider the scrambling time which is defined as the minimum time for the recovery of the information which can be retrieved from the Hawking radiation after falling into the black hole according to the Hayden-Preskill protocol [89]. Recall that, in the entanglement wedge construction according to the prescription of the island, the density matrix of the Hawking radiation is represented by the union of the radiation region R − ∪ R + and the island. This implies that the information of the signal which is thrown into the island is not contained in the black hole but it could be decoded by the Hawking radiation. If an observer sends a signal from the cutoff surface, it would reach the island after an earliest duration t scr identified as the scrambling time defined as The leading order contribution for the scrambling time is proportional to the inverse of the Hawking temperature and the logarithm of the black hole entropy, which is consistent with the result obtained in Ref. [90]. In this way, the expression of the scrambling time at the leading order is universal. Furthermore, we observe that the scrambling time is very small compared to the Page time. In order to see the explicit influence of the nonzero mass of graviton on the scrambling time, we expand the expression for t scr as where the first term is the prediction of Einstein gravity. This expression implies that the presence of the graviton mass would reduce the scrambling time if the ratio of the massive gravity coupling parameters c 1,2 satisfies

Entanglement entropy for charged black holes
In this section, we shall compute the entanglement entropy of Hawking radiation with the absence and presence of the islands, the Page time, and the scrambling time for the nonextremal and extremal charged black holes in massive gravity. The arguments and the method which are presented in the case of the neutral black hole can still be applied for the situation of the nonzero electric charge. There are some differences for the extremal charged black hole, which are due to its Penrose diagram to be a one-sided geometry, rather than the two-sided geometry like the Penrose diagram of the neutral and non-extremal charged black holes.

Non-extremal case
First, we consider the non-extremal charged black hole whose metric in the Kruskal coordinate is given by (22). The Penrose diagram without and with the island configuration is depicted in the left and right panels of Fig. 6, respectively. For the case of the configuration with no island, the entanglement entropy is obtained as where we have used the approximation r h b. In analogy to the neutral black hole which we study in the previous section, the entanglement entropy for the non-extremal charged black hole with no island grows linearly with time in the limit of 1/κ + t b and thus it becomes infinite at the late stage of the evaporation. Of course, this is inconsistent with the unitarity time evolution. Therefore, we expect the emergence of the island configuration which minimizes the generalized entropy by which the entanglement entropy stops the linear increasing and reaches a saturation value. Now we arrive at the case of the configuration with one island. The generalized entropy with an island at the early times reads where the ellipses refer to the terms independent on the temporal and spatial location of the island boundaries. By extremizing S et gen with respect to a, we find which implies that the size of the island in the beginning of the black hole evaporation is in order of the Planck length.
On the other hand, no island emerges at the early times and thus the behavior of the entanglement entropy is governed by the configuration with no island. Next we study the late time behavior of the entanglement entropy in the presence of the island. The generalized entropy in this situation is found as where we have used the fact that t a = t b extremizes the generalized entropy. In order to find the spatial location of the island boundaries, we extremize the generalized entropy S lt gen with respect to a where the boundaries of the island are located in the vicinity of the event horizon. The solution for a = r h + + O( 2 ) to the subleading order approximation is derived as Here, we see that in analogy to the case of the neutral black hole the subleading term of the location of the island boundaries is suppressed by (cG N /r 2 + ) 2 and the higher order terms should be strongly suppressed by the further powers of cG N /r 2 + . Then, by inserting the location of the island bound- Fig. 6 The Penrose diagram of the eternal non-extremal charged black hole in massive gravity with the no island configuration (left panel) and the island configuration (right panel). The radiation region is the union of two parts R ± whose boundaries are denoted by b ± . I refers to the island with the boundaries represented by a ± aries just obtained into S lt gen , the entanglement entropy for the non-extremal black hole is finally determined as This result indicates that the entanglement entropy at the late times would reach a saturation value whose leading order is the twice of the Bekenstein-Hawking entropy of the black hole, i.e. S E E 2πr 2 + /G N = 2S BH . Beside, it obtains the logarithm corrections coming form the quantum matter. Other correction terms are suppressed by the powers of cG N /r 2 + . We expand the leading contribution of the entanglement entropy S E E for the non-extremal black hole in terms of the graviton mass as where r 0± = M ± M 2 − Q 2 are the the radii of the event and Cauchy horizons of the RN black hole corresponding to Einstein gravity or the case of massless graviton. Similar to the neutral black hole, if the event horizon radius of the corresponding RN black hole (in Einstein gravity) is larger than the ratio −2c 2 /c 1 , the nonzero mass of graviton would reduce the entanglement entropy compared to the massless case. The Page time for the non-extremal charged black hole is easily derived as t Page = 3S BH /π cT H where T H = κ + /2π from equating the linear growth of the entanglement entropy at the early times with the asymptotic constant value 2S BH . We expand the Page time in terms of m 2 g as where the first term is the Page time associated with the RN black hole and the second and higher order terms are the contributions arising due to the nonzero mass of graviton. We observe that if the following relation is satisfied, the entanglement entropy in massive gravity would take a shorter time to reach the saturation value, and on the contrary it would take a longer time for the emergence of the island. As discussed in the case of the neutral black hole, this is because massive gravity with the coupling parameters satisfying (64) reduces the bound of the entanglement entropy or/and increases the Hawking temperature causing the black hole radiating more rapidly. However, there is an important difference between the neutral black hole and the non-extremal charged black hole in massive gravity. With respect to the neutral black hole, the Page time in massive gravity is always smaller than that in Einstein gravity for the coupling parameters c 1,2 which are both positive. But, with respect to the non-extremal charged black hole, this only happens if r 0+ (or the black hole mass) is sufficiently far from r 0− (or the electric charge of the black hole). Whereas, when r 0+ is sufficiently close to r 0− , the second term in (63) is positive which means that the Page time for the non-extremal charged black hole in massive gravity is larger than that in Einstein gravity.
We compute the scrambling time for the non-extremal charged black hole in massive gravity, which is obtained as This expression implies that the scrambling time is very small compared to the Page time and its behavior in the parameters of massive gravity is basically the same as the Page time for the non-extremal charged black hole discussed above.
In the limit of the vanishing electric charge or the Cauchy horizon radius approaching zero, one can easily see that the size of the island, the entanglement entropy, the Page time, and the scrambling time reproduce the corresponding results for the neutral black hole as obtained in the previous section. Also, in the limit of the massless graviton m g → 0, our calculations for these quantities would reduce the corresponding computations for the RN black hole reported by the authors in [27].

Extremal case
For the extremal charged black hole, the metric in the Kruskal coordinate is given in Eq. (28) and the Penrose diagram is shown in Fig. 7. Compared to the non-extremal case, the radiation region is only given by one region R + . In the absence of the island, the entanglement entropy is computed as where b 0 = (t b , 0) is a reference point which represents the touch of the Cauchy surface at the singularity. Unfortunately, the conformal factor W 2 (0) is ill-defined because this function is divergent at r = 0. As a result, the computation of the entanglement entropy as well as the Page time in the absence of the island would lead to the ill-defined result in the extremal case. However, we hope that the UV corrections at the very short distances such as the nonperturbative renormalization group in the quantization of Einstein gravity [91][92][93], string T-duality [94] or noncommutative geometry [95] would make the conformal factor W 2 (r ) in the Kruskal coordinate behaving smoothly as r approaches zero. And, since it would be able to yield a finite result for the computation of the entanglement entropy and the Page time. This would be further investigated in our future works.
In the presence of an island, the singularity at r = 0 as discussed above can be avoided because in this situation the entanglement entropy is proportional to the logarithm of the geodesic distance between the boundary of the island and the cutoff surface. More specifically, the generalized entropy for the extremal case with the contribution of one island is computed as |a − r e ||b − r e | ab ×(a + 2r e + δ) where  Fig. 7 The Penrose diagram of the eternal extremal charged black hole in massive gravity without island (left panel) and with one island (right panel). The radiation region is represented by R + with the boundary b + . I denotes the island extending from the singularity to its boundary a + It is straightforward to see that the extremal condition ∂ S gen /∂t a = 0 leads to t a = t b . In addition, we assume that the location of the island boundary is slightly outside the horizon r e of the black hole, i.e. a ≈ r e , as well as the cutoff surface of the radiation region is far away from the horizon r e of the black hole. With this assumption b r e ≈ a, we have the following approximation 1.
By extremizing the generalized entropy within this approximation with respect to a, we find the location of the island boundary as where the subleading term is suppressed by square root of cG N /r 2 e , rather than the second power as in the cases of the neutral and non-extremal charged black holes. Then, the entanglement entropy reads where the terms relating to the first order of c are the logarithm corrections which are small compared to the subleading terms. We observe that the entanglement entropy for the extremal charged black hole becomes a finite constant at the late times of the evaporation, like the neutral and non-extremal charged black holes. However, there here is an important difference that the asymptotic constant value for the entanglement entropy with respect to the extremal charged black hole is approximately the Bekenstein-Hawking entropy (rather than the twice of the Bekenstein-Hawking entropy), i.e. S E E πr 2 e /G N = S B H . This is due to the essential difference in their causal structure: the Penrose diagram of the neutral and nonextremal charged black holes is the two-sided geometry, whereas the Penrose diagram of the extremal charged black hole is the one-sided geometry. This difference also implies that one cannot derive the entanglement entropy for the extremal charged black hole by taking the continuous extremal limit r + → r − of the non-extremal charged black hole.

Conclusion
There are the arguments [37,38] which implied that calculating the entanglement entropy of Hawking radiation and the Page curve for black holes using the island method would be considered in the situation of massive graviton. This is one of the reasons that we consider the entanglement entropy of Hawking radiation and the Page curve using the island method in the framework of dRGT massive gravity which is considered as a candidate for a consistent theory of gravity accompanied with a nonzero mass of graviton. Furthermore, an interesting aspect of black holes in the context of dRGT massive gravity is that the behavior of black hole geometries is asymptotically non-flat. For various black hole geometries corresponding to the asymptotically flat behavior, it was pointed to that the island is emerged at the late time of the black hole evaporation where the boundary of the island is very close the black hole horizon and the island extends almost through the whole black hole interior. In this work, we also show the emergence of the island at the late time of the evaporation the same as the cases of the asymptotically flat black hole geometries. This result and the previous confirmations thus support the emergence of the island at the late time of the evaporation as a universal feature of the semiclassical description of the black hole evaporation: the computations based on the semiclassical description of the black hole evaporation can create replica wormholes or the island where the information is stored; since the black hole evaporation follows the Page curve consistent with the unitarity evolution.
More explicitly, we calculated the entanglement entropy of Hawking radiation emitted by the eternal black holes, the corresponding Page curve, and the scrambling time in the context of massive gravity whose free parameters are the graviton mass and two coupling parameters c 1,2 in four dimensions. In our calculations, we first employ the island rule which was recently obtained by using the quantum extremal surface technique as well as from the gravitational Euclidean path integral using the replica trick. According to the prescription of the quantum extremal surface, the entanglement entropy is derived as the minimum value of the generalized entropy which is a sum of the Bekenstein-Hawking entropy of the island boundaries and the von Neumann entropy of the matter fields on the union of the radiation region and the island. In order to find the minimum value of the generalized entropy, we need to extremize the generalized entropy over all possible boundary surfaces of the island. Second, we use the approximation of the two-dimensional conformal field theory where the contribution of the matter sector to the entanglement entropy is easily computed as the logarithm of the disjoint intervals, when the distance between the cutoff surfaces of the radiation region is sufficiently large compared to the scale of the size of the cutoff surfaces.
For the neutral and non-extremal black holes, we indicate that the island does not appear at the early times of the black hole evaporation or in other words the no island configuration is the dominant contribution to the entanglement entropy at the early stage of the evaporation. Hence, the entanglement entropy grows linearly with time. However, as the amount of Hawking radiation becomes sufficiently large at the late times, the island emerges slightly outside the event horizon of the black holes and becomes the preferred configuration by which the growth of the entanglement entropy reaches a finite saturation value which is the twice of the Bekenstein-Hawking entropy of the black hole at the leading order. Furthermore, we calculate the Page time and scrambling time which are approximately given by 3S B H π cT H and log S B H 2π T H , respectively, which is universal for various black hole geometries [25,[27][28][29][30][31][32]. Whereas, for the extremal charged black hole, the entanglement entropy at the early times without the island is ill-defined because the conformal factor of the metric in the Kruskal coordinate is divergent at r = 0 corresponding to that the Cauchy surface hits the curvature singularity. This implies that we need to consider new physics in the UV region which would make the metric behaving smoothly at the very short distances, which will be further investigated in our future works. At the late times when the island is formed, the entanglement entropy for the extremal charged black hole reaches the saturation value which is the Bekenstein-Hawking entropy. The differences between the extremal charged black hole and the non-extremal charged (and neutral) black hole are due to their Penrose diagram: the Penrose diagram of the neutral and non-extremal charged black holes is the two-sided geometry, whereas it is the onesided geometry for the extremal charged black hole. We also show that the corresponding results of the Schwarzschild and Reissner-Nordström black holes are restored in the limit of that the graviton mass goes to zero. In addition, we study the impact of the parameters of massive gravity on the size of the island, the entanglement entropy, the Page time, and the scrambling time. The direct impact of massive gravity can be ignored compared to its indirect impact which modifies the horizon radius of the black hole. An analytic investigation can be performed in the region of the small graviton mass which is consistent with the presently experimental constraint. If the black hole mass M for the case of the zero electric charge (or a combination of the mass and the electric charge of the black hole for the case of the nonzero electric charge) is larger than −c 2 /c 1 , the bound of the entanglement entropy in massive gravity is lower than that in Einstein gravity, and this happens the contrary for M < −c 2 /c 1 . This suggests that, with the nonzero mass of graviton and the proper coupling parameters, it takes a shorter duration in order for the entanglement entropy reaching the saturation value. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical work, since there is no data associated with this article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .