1 Introduction

The connection between thermodynamics and gravity has attracted great attention in recent decades. In 1994, Jacobson showed that the Einstein equations could be derived from fundamental thermodynamic relation which hold for local Rindler causal horizons [1]. Verlinde put forward the interesting but incomplete viewpoint that gravity could be explained as entropy force [2]. Meanwhile, a series of papers proved that the maximum entropy principle for a perfect fluid holds in different theories. This principle shows that if the constraint equation and some thermodynamic conditions are satisfied, the gravitational field equations could be derived by the extrema of the total entropy [3,4,5,6,7,8,9]. Recently, Jacobson proposed the “causal diamond” structure and considered entanglement equilibrium to imply the Einstein equations [10]. All of these studies suggest a general and solid connection between thermodynamics and gravity. Moreover, the emergence has also been widely discussed in the past years [11, 12]. It was considered that gravity may not be the fundamental assumption describing nature. However, most of these discussions of the connection between gravity and thermodynamics are focused on establishing the relation between the first order variation of the thermodynamical properties and the gravitational equations. This manuscript will reveal their connection in higher order variation by investigating the stability of a perfect fluid in static background spacetime.

In general there are two methods of testing the stability of a static configuration of a fluid [13]. One is the dynamical method, and the other is the thermodynamical method. In the dynamical method, one requires that the physical quantities deviate only slightly from equilibrium state. Assuming that the motions are adiabatic and reversible, then the variation separation approach makes the perturbational fields \(\psi ^a\) take the form \(\psi ^a(r)e^{-i\omega t}\). Then the equation for \(\psi ^a\) could be transformed to a Sturm–Liouville eigenvalue problem, with \(\omega \) the eigenvalue. The stability is tested by ascertaining whether or not all the \(\omega ^2\) are positive. The stability problem in general relativity was first discussed by Chandrasekhar [14, 15]. After that, Friedman investigated the dynamic stability of relativistic stars with respect to perturbations that arise in the Lagrangian displacement framework [16,17,18]. Recently, Wald and Seifert presented a general dynamical method for the analysis of the stability of static, spherically symmetric solutions to spherically symmetric perturbations [19]. This method could be used in an arbitrary diffeomorphism covariant Lagrangian field theory in which the field equations have at most second order derivatives of the metric [20]. In the thermodynamical method, describing an isolated system, the negative of the second variation of the total entropy corresponding to thermodynamical stability was discussed by Cocke [13] and Sorkin [3]. Recently, Roupas proved that the maximum of the total entropy for a perfect fluid gives the same criterion as Yabushita’s result with some additional conditions added [21,22,23].

However, using the dynamical method to solve the stability problem in static spacetime always will be limited to spherical symmetry and radial perturbations. For the cases without spherical symmetry, or with more general perturbations, the dynamical method is very hard to deal with. We believe that the thermodynamical method is promising as we try to solve this difficulty. In this manuscript, by extending Roupas’ result, we obtain a general formula for the second variation of the total entropy for a perfect fluid, as the criterion for thermodynamical stability. It should be noted that this criterion could be applied in a general static background without spherical symmetry.

It is worth noting that recently Wald et al. presented a very comprehensive discussion of the equivalence of dynamical stability and thermodynamical stability [24]. The results in [24] seem to be similar to a part of our manuscript; however, our assumptions and arguments are different from [24]. For instance, with the definition of the ADM mass, a crucial assumption in [24] is that the spacetime should be asymptotically flat, while our derivations apply to any region, imposing no global conditions on spacetime.

The rest of this paper is organized as follows. In the next section, we first of all give a formula for the thermodynamical stability criterion for a perfect fluid in static background spacetime. Then this criterion is applied specifically to the case of radial perturbations of static, spherical symmetric perfect fluid in Einstein gravity. It is found that our thermodynamical criterion is consistent with the dynamical criterion found by Wald. In Sect. III, the explicit expression of the thermodynamical stability criterion of a perfect fluid under non-radial perturbation in general static background is obtained. Finally, we summarize this paper with some comments and discussions.

Throughout our discussion, units will be used in which \(c=G=1\). The letters (a, b, c) denote the abstract index. We also ignore the factor \(\kappa =8\pi \) in the Einstein equations, hence \(G_{\mu \nu }=T_{\mu \nu }\).

2 Criterion for the thermodynamical stability

First, we briefly review the proof of the maximum entropy principle for a perfect fluid in general static spacetime [5]. We assume that all the quantities are measured by static background observers which are orthogonal to the hypersurface \(\Sigma \), and we consider the fluid over any selected region C on \(\Sigma \) to satisfy the ordinary thermodynamic relations and Tolman’s law, which gives \(T\chi ={\text {const.}}\), where T and \(\chi \) are the temperature of the fluid and the redshift factor, respectively. Without loss of generality, we take \(T\chi =1\). In [5] we have shown that if the constraint equation and some boundary conditions are satisfied, and the total particle number is fixed, the variation of the total entropy is

$$\begin{aligned} \delta S = \int _C\frac{{\sqrt{h}}}{T}\delta \rho +\frac{p+\rho }{T}\delta {\sqrt{h}}\,, \end{aligned}$$
(1)

where \(\rho \) and p represent the energy density and the pressure of the fluid, respectively, and h is the determinant of the induced metric \(h_{ab}\) on \(\Sigma \). Then we proved that the extrema of the total entropy \(\delta S=0\) gives the components of gravitational field equations [5]. The maximum entropy principle implies that the gravitational field equations may be replaced by a constraint equation and thermodynamic relations.

An isolated star in thermodynamical equilibrium is said to be thermodynamically stable if \(\delta ^2S<0\). When discussing the second variation of the total entropy, we consider the Tolman law also to be valid, since the state deviates slightly from the equilibrium state. Applying the local first law of thermodynamics,

$$\begin{aligned} T\delta s=\delta \rho -\mu \delta n, \end{aligned}$$
(2)

and the integrated form of the Gibbs–Duhem relation

$$\begin{aligned} p+\rho =Ts+\mu n, \end{aligned}$$
(3)

with the fact that \(\mu /T=\mathrm{{const.}}\) [5], we obtain

$$\begin{aligned} \delta p=s\delta T+n\frac{\mu }{T}\delta T=\frac{p+\rho }{T}\delta T \,, \end{aligned}$$
(4)

where \(\mu \) and n denote chemical potential and particle number density, respectively. So the second variation of the total entropy could be written as

$$\begin{aligned} \delta ^2S= & {} \int _C\frac{1}{T}\bigg [2\delta \rho \delta {\sqrt{h}}+{\sqrt{h}}\delta ^2\rho \nonumber \\&+(p+\rho )\delta ^2{\sqrt{h}}-\frac{\delta p\delta \rho }{p+\rho }{\sqrt{h}}\bigg ]\,. \end{aligned}$$
(5)

\(\delta ^2S<0\) means the system is thermodynamically stable. Hence, using Eq. (5) one could directly obtain the specific form of the stability criterion. A natural question is whether the thermodynamical stability is equivalent to dynamical stability or not.

As a concrete example, we investigate whether this equivalence is valid for spherically symmetric perturbations of a static, spherical symmetric perfect fluid in Einstein gravity. For such a spacetime, the metric takes the form [25]

$$\begin{aligned} \mathrm{d}s^2=-e^{2\Phi (t,r)}\mathrm{d}t^2+e^{2\Lambda (t,r)}\mathrm{d}r^2+r^2\mathrm{d}\Omega ^2. \end{aligned}$$
(6)

According to Chandrasekhar’s procedure, under the perturbation, the four-velocity becomes [15]

$$\begin{aligned}&u_0=-e^{\Phi }, u^0=e^{-\Phi }, \nonumber \\&u_1=e^{2\Lambda -\Phi }\mathbf {v}, u^1=e^{-\Phi }\mathbf {v}, \end{aligned}$$
(7)

where \(\mathbf {v}=\mathrm{d}r/\mathrm{d}t=\partial \xi /\partial t\). Here \(\xi \) is the radial “Lagrangian displacement”, which describes the radial displacement of each fluid element from its “equilibrium position”. Then the tr component of the Einstein equations, i.e. \(T_0^1=G_0^1\), gives

$$\begin{aligned} -(p+\rho )e^{2\Lambda }\frac{\partial \xi }{\partial x^0}=\frac{2}{r}\frac{\partial \Lambda }{\partial x^0}. \end{aligned}$$
(8)

Direct integration yields

$$\begin{aligned} \delta \Lambda =-\frac{r}{2}e^{2\Lambda }(p+\rho )\xi . \end{aligned}$$
(9)

With the first variation of constraint equation, \(\delta G_{00}=\delta T_{00}\), we have

$$\begin{aligned} \delta \rho =-\frac{1}{r^2}\frac{\partial }{\partial r}[r^2(p+\rho )\xi ]. \end{aligned}$$
(10)

Therefore, the second variation of \(\rho \) could be written as

$$\begin{aligned} \delta ^2\rho = -\frac{1}{r^2}\frac{\partial }{\partial r}[r^2(\delta p+\delta \rho )\xi +r^2(p+\rho )\delta \xi ]. \end{aligned}$$
(11)

Meanwhile, the variation of the induced metric could be written as

$$\begin{aligned} \delta {\sqrt{h}}= e^{\Lambda }r^2\sin \theta \cdot \delta \Lambda = -\frac{1}{2}e^{3\Lambda }r^3(p+\rho )\sin \theta \cdot \xi \end{aligned}$$
(12)

and

$$\begin{aligned} \delta ^2{\sqrt{h}}= & {} \frac{3}{4}r^4e^{5\Lambda }(p+\rho )^2 \sin \theta \cdot \xi ^2 \nonumber \\&-\frac{1}{2}e^{3\Lambda }r^3(\delta p+\delta \rho )\sin \theta \cdot \xi \nonumber \\&-\frac{1}{2}e^{3\Lambda }r^3(p+\rho )\sin \theta \cdot \delta \xi . \end{aligned}$$
(13)

From Eq. (6), \(1/T=\chi =\sqrt{-g_{00}}=e^{\Phi }\), using integration by parts and dropping the boundary terms, the first three terms in the right hand side of Eq. (5) could be calculated one by one. Note that in this case \(\int _C\) becomes \(4\pi \int _r\mathrm{d}r\). Explicitly, the first term becomes

$$\begin{aligned}&4\pi \int _r \mathrm{d}r\frac{2}{T}\delta \rho \delta {\sqrt{h}}= 4\pi \int _r \mathrm{d}r \frac{e^{\Phi +3\Lambda }}{2r}\frac{\partial }{\partial r}[r^2(p+\rho )\xi ]^2 \nonumber \\&\quad = 4\pi \int _r \mathrm{d}r\bigg [-\frac{3}{4}e^{\Phi +5{\Lambda }}r^4(p+\rho )^3\xi ^2\nonumber \\&\qquad +\frac{1}{2}e^{\Phi +3\Lambda }\bigg (\frac{2\Phi '}{r}+\frac{1}{r^2}\bigg )r^4(p+\rho )^2 \xi ^{2} \bigg ], \end{aligned}$$
(14)

the second term becomes

$$\begin{aligned}&4\pi \int _r \mathrm{d}r\frac{1}{T}{\sqrt{h}}\delta ^2\rho = 4\pi \int _r \mathrm{d}r e^{\Phi +\Lambda }r^2 \cdot \frac{1}{r^2}\frac{\partial }{\partial r}\nonumber \\&[r^2(\delta p+\delta \rho )\xi +r^2(p+\rho )\delta \xi ] \nonumber \\&\quad = 4\pi \int _r \mathrm{d}r e^{\Phi +3\Lambda }\bigg [\frac{1}{2}r^3(p+\rho )(\delta p+\delta \rho )\xi \nonumber \\&\qquad +\frac{1}{2}r^3(p+\rho )^2\delta \xi \bigg ], \end{aligned}$$
(15)

and the third term becomes

$$\begin{aligned}&4\pi \int _r \mathrm{d}r\frac{1}{T}(p+\rho )\delta ^2{\sqrt{h}}=4\pi \int _r \mathrm{d}r\bigg [\frac{3}{4}r^4e^{\Phi +5\Lambda }(p+\rho )^3\xi ^2\nonumber \\&\quad -\frac{1}{2}e^{\Phi +3\Lambda }r^3(\delta p+\delta \rho )(p+\rho )\xi \nonumber \\&\quad -\frac{1}{2}e^{\Phi +3\Lambda }r^3(p+\rho )^2\delta \xi \bigg ]. \end{aligned}$$
(16)

Substituting these results into Eq. (5) we obtain

$$\begin{aligned} \delta ^2S= & {} 4\pi \int _r \mathrm{d}r\bigg \{\frac{1}{2}e^{\Phi +3\Lambda }\bigg (\frac{2\Phi '}{r}+\frac{1}{r^2}\bigg )r^4(p+\rho )^2\xi ^2\nonumber \\&-e^{\Phi +\Lambda }r^2\frac{\delta p\delta \rho }{p+\rho }\bigg \}. \end{aligned}$$
(17)

Note that we have chosen the energy density \(\rho \) and the particle number density n as independent variables on obtaining the first variation of the total entropy [5]. But Chandrasekhar has chosen \(\rho \) and p as independent variables [14, 15]. It is not obvious how to compare our result with Chandrasekhar’s. However, using the dynamical method, Seifert and Wald give a stability criterion of the star with a “barotropic” equation of the state of the form \(\rho =\rho (n)\); in this situation there is only one thermodynamical variable. We will show that our stability criterion would exactly reduce to Wald’s result. For this purpose, we start with the Lagrangian for the perfect fluid used in [19],

$$\begin{aligned} \mathcal {L}_{\text {mat}}=-\varrho (n). \end{aligned}$$
(18)

Wald showed that there exists an identification:

$$\begin{aligned} \rho \rightarrow \varrho , ~\ p \rightarrow \frac{\partial \varrho }{\partial n}n-\varrho . \end{aligned}$$
(19)

It is thus easy to obtain

$$\begin{aligned} \delta \rho = \frac{\partial \varrho }{\partial n}\delta n, ~\ \delta p = \frac{\partial ^2\varrho }{\partial n^2}n\delta n. \end{aligned}$$
(20)

So Eq. (17) becomes

$$\begin{aligned} \delta ^2 S= & {} 4\pi \int _r \mathrm{d}r\bigg [r^2e^{\Phi +\Lambda }\bigg (2\frac{\partial \Phi }{\partial r}+\frac{1}{r}\bigg )\bigg (\frac{\partial \Lambda }{\partial r}\nonumber \\&+\frac{\partial \Phi }{\partial r}\bigg )(p+\rho )\xi ^2-e^{\Phi +\Lambda }r^2\frac{\partial ^2\varrho }{\partial n^2}(\delta n)^2\bigg ]. \end{aligned}$$
(21)

Note that the variation of particle number density n could be written as [19]

$$\begin{aligned} \delta n= & {} n\bigg [\frac{\partial \xi }{\partial r}+\bigg (\frac{1}{\nu }\frac{\partial \nu }{\partial r}+\frac{\partial \Lambda }{\partial r}+\frac{2}{r}\bigg )\xi \nonumber \\&\quad -\,\bigg (\frac{\partial \Phi }{\partial r}+\frac{\partial \Lambda }{\partial r}\bigg )\xi \bigg ]=\frac{e^{\Phi }}{r^2}\frac{\partial }{\partial r}(r^2e^{-\Phi }n \xi ). \end{aligned}$$
(22)

We obtain

$$\begin{aligned} \delta ^2S= & {} 4\pi \int _r \mathrm{d}r \bigg [r^2e^{\Phi +\Lambda }\bigg (2\frac{\partial \Phi }{\partial r}+\frac{1}{r}\bigg )\nonumber \\&\times \, \bigg (\frac{\partial \Lambda }{\partial r}+\frac{\partial \Phi }{\partial r}\bigg )(p+\rho )\xi ^2 -\frac{e^{3\Phi +\Lambda }}{r^2}\frac{\partial ^2\varrho }{\partial n^2} \bigg (\frac{\partial }{\partial r}(r^2e^{-\Phi }n\xi )\bigg )^2\bigg ].\nonumber \\ \end{aligned}$$
(23)

Note that the thermodynamical stability requires only the second variation of the total entropy to be negative. The above expression, Eq. (23), agrees with the result in Seifert and Wald. So far, we have shown that our criterion for the thermodynamical stability is consistent with the dynamical stability criterion for spherical perturbations of a static, spherical symmetry star with a barotropic equation of state. In the next section, we will give the stability criterion in more general cases.

3 Thermodynamical stability criterion in general cases without spherical symmetry

The criterion for the thermodynamical stability applied not only for the above particular case, but also for more general case. In static spacetime, for the cases without spatial symmetry, or with non-radial perturbations, the method of dynamical stability is hard to deal with. However, thermodynamical stability seems easier to handle these cases. Assuming that the metric of a perfect fluid star in a background spacetime could be written as

$$\begin{aligned} \mathrm{d}s^2=g_{00}\mathrm{d}t^2+g_{ij}\mathrm{d}x^i\mathrm{d}x^j, \end{aligned}$$
(24)

and the perturbation fields to be \(\delta g_{\mu \nu }\), we show how to get the explicitly form of the second variation of the total entropy \(\delta ^2S\) in a general static spacetime. In this section, \(h_{ab}\) denotes the induced metric of the \(t=\mathrm{{const.}}\) slice \(\Sigma \). We use \(A_a\), \(D_a\) and \(\Box \) to denote the four-acceleration of the observer, the 3-dimensional covariant derivative and 3-dimensional Box operator \(D_aD^a\) in \(\Sigma \), respectively.

The extrinsic curvature of \(\Sigma \) is

$$\begin{aligned} K_{ab}=h_a^c\nabla _c u_b. \end{aligned}$$
(25)

The relation between the ordinary curvature R and the 3-dimensional curvature \(R^{(3)}\) on \(\Sigma \) is given by

$$\begin{aligned} R^{(3)}=R+2R_{ab}u^au^b-\frac{1}{2}(K_a^a)^2+\frac{1}{2}K_{ab}K^{ab}, \end{aligned}$$
(26)

which yields

$$\begin{aligned} \rho =\frac{1}{2}R^{(3)}+\frac{1}{2}(K_a^a)^2-\frac{1}{2}K_{ab}K^{ab} \,. \end{aligned}$$
(27)

Note that \(K_{ab}|_0=0\) in static background spacetime. We obtain

$$\begin{aligned} \delta \rho =\frac{1}{2}\delta R^{(3)}\end{aligned}$$
(28)

and

$$\begin{aligned} \delta ^2\rho= & {} \frac{1}{2}\delta ^2R^{(3)}+h^{ac}h^{bd}(\delta K_{ac}\delta K_{bd}-\delta K_{ab}\delta K_{cd}). \end{aligned}$$
(29)

The perturbation of the induced metric is

$$\begin{aligned} \delta {\sqrt{h}}=\frac{1}{2}{\sqrt{h}}h^{ab}\delta h_{ab}. \end{aligned}$$
(30)

Then we calculate the variation of extension curvature, \(\delta K_{ab}\),

$$\begin{aligned} \delta K_{ab}= & {} \delta (h_a^c\nabla _cu_b) \nonumber \\= & {} \nabla _cu_b\delta h_a^c+h_a^c\nabla _c\delta u_b-h_a^c\delta \Gamma ^{\mu }_{cb} u_{\mu }. \end{aligned}$$
(31)

Using

$$\begin{aligned} \delta u_a=-\frac{1}{2}u_au^bu^c\delta g_{bc}, \end{aligned}$$
(32)

and the fact \(\nabla _au_b=-A_bu_a\) in a static background, we have

$$\begin{aligned} h^{ac}h^{bd}\delta K_{ac}\delta K_{bd}=h^{ab} h^{cd}u_{e}u_{f}\delta \Gamma ^{e}_{ab}\cdot \delta \Gamma ^{f}_{cd} \end{aligned}$$
(33)

and

$$\begin{aligned} -h^{ac}h^{bd}\delta K_{ab}\delta K_{cd}=-h^{ac}h^{bd}u_{e}u_{f}\delta \Gamma ^{e}_{ab}\cdot \delta \Gamma ^{f}_{cd}. \end{aligned}$$
(34)

Substituting Eqs. (33) and (34) into Eq. (29), the expression of \(\delta ^2\rho \) could be obtained.

Now we calculate all terms in Eq. (5) one by one. Noting that the induced metric and its derivatives are fixed on the boundary of the selected region C, we could use integration by parts and drop the boundary terms. The first term of Eq. (5) could be written as

$$\begin{aligned} {\mathcal {G}}_{1}= & {} \int _{C} \frac{2}{T}{\delta }{\rho }{\delta }{{\sqrt{h}}} \nonumber \\= & {} \int _{C}\frac{\chi }{2}{{\sqrt{h}}} \bigg ( h^{cd}\delta h_{cd}\cdot D^{a}D^{b}\delta h_{ab} - h^{ab}h^{cd}\delta h_{cd} \nonumber \\&\cdot D^{e}D_{e}\delta h_{ab}-h^{cd}R^{(3) {ab}} \delta h_{ab}\delta h_{cd} \bigg ). \end{aligned}$$
(35)

The second term of Eq. (5) becomes

$$\begin{aligned} \mathcal {G}_2= & {} \int _C \frac{{\sqrt{h}}}{T} \delta ^2\rho \nonumber \\= & {} \int _C \frac{\chi {\sqrt{h}}}{2}\delta \bigg (h^{ab}\delta R^{(3)}_{ab}+R^{(3)}_{ab}\delta h^{ab}\bigg )\nonumber \\&+\,\chi {\sqrt{h}}(h^{ab}h^{cd}-h^{ac}h^{bd}) u_{\mu }u_{\nu }\delta \Gamma ^{\mu }_{ab}\cdot \delta \Gamma ^{\nu }_{cd} \nonumber \\= & {} \int _C \frac{\chi {\sqrt{h}}}{2}\delta \bigg [h^{ac}h^{bd}D_c D_d\delta h_{ab}\nonumber \\&-\,h^{ab}D^cD_c\delta h_{ab}+R^{(3)}_{ab}\delta h^{ab} \bigg ] \nonumber \\&+\,\chi {\sqrt{h}}(h^{ab}h^{cd}-h^{ac}h^{bd}) u_{e}u_{f}\delta \Gamma ^{e}_{ab}\cdot \delta \Gamma ^{f}_{cd} \nonumber \\= & {} \int _C \frac{\chi {\sqrt{h}}}{2} \bigg [h^{bd}\delta h^{ac}\cdot D_c D_d\delta h_{ab}+h^{ac}\delta h^{bd}\nonumber \\&\cdot D_c D_d\delta h_{ab} +h^{ac}h^{bd}\delta (D_c D_d\delta h_{ab}) \nonumber \\&-\,\delta h^{ab}\cdot D^cD_c\delta h_{ab}-h^{ab}\delta (D^cD_c\delta h_{ab})\nonumber \\&+\,\delta R^{(3)}_{ab}\delta h^{ab}+R^{(3)}_{ab}\delta ^2h^{ab} \bigg ] \nonumber \\&+\,\chi {\sqrt{h}}(h^{ab}h^{cd}-h^{ac}h^{bd}) u_{e}u_{f}\delta \Gamma ^{e}_{ab}\cdot \delta \Gamma ^{f}_{cd}, \end{aligned}$$
(36)

where \(h^{ac}h^{bd}\delta (D_c D_d\delta h_{ab})\) and \(-h^{ab}\delta (D^cD_c\delta h_{ab})\) could be calculated as

$$\begin{aligned}&h^{ac}h^{bd}\delta (D_c D_d\delta h_{ab}) \nonumber \\&\quad = h^{ac}h^{bd}D_c D_d\delta ^2 h_{ab}-\frac{1}{2}\Box \delta h^{ab}\cdot \delta h_{ab}\nonumber \\&\qquad +\,D_aD^b\delta h_{bc}\cdot \delta h^{ac} -\frac{1}{2}h^{ab}D_c D_d\delta h_{ab}\cdot \delta h^{cd} \nonumber \\&\qquad +\,\frac{1}{2}D^{c}\delta h^{ab}\cdot D_c\delta h_{ab}-2h^{ab}D^{c}\delta h_{ac}\nonumber \\&\qquad \cdot D^{d}\delta h_{bd} +h^{ab}D^{c}\delta h_{ab}\cdot D^{d}\delta h_{cd}\nonumber \\&\qquad +\,D^{c}\delta h^{ab}\cdot D_b\delta h_{ac} \end{aligned}$$
(37)

and

$$\begin{aligned}&-h^{ab}\delta (D^cD_c\delta h_{ab}) = -\delta [D^cD_c(h^{ab}\delta h_{ab})]\nonumber \\&\qquad +\, \delta h^{ab}\cdot D^cD_c\delta h_{ab} \nonumber \\&\quad = -\, [D_c\delta (D^{c}(h^{ab}\delta h_{ab}))+\delta C^c_{cd}D^{d}(h^{ab}\delta h_{ab})]\nonumber \\&\qquad +\, \delta h^{ab}\Box \delta h_{ab} \nonumber \\&\quad = -\, h^{ab}D_c\delta h^{cd}\cdot D_d\delta h_{ab}-\Box (\delta h^{ab}\delta h_{ab})-h^{ab}\Box \delta ^2 h_{ab} \nonumber \\&\qquad -\,\frac{1}{2}h^{ab}h^{cd}D^e\delta h_{ab}D_e\delta h_{cd}+\delta h^{ab}\Box \delta h_{ab}. \end{aligned}$$
(38)

So we obtain

$$\begin{aligned} \mathcal {G}_2= & {} \int _C \frac{\chi {\sqrt{h}}}{2} \Big [2\delta h^{ac}D_c D^{b}\delta h_{ab}+2\delta h^{bd}D^{a} D_d\delta h_{ab}\nonumber \\&+\,h^{ac}h^{bd}D_c D_d\delta ^2 h_{ab}-3\delta h_{ab}\Box \delta h^{ab}-h^{ab}D_c D_d\delta h_{ab}\nonumber \\&\cdot \delta h^{cd}-\frac{3}{2}D^{c}\delta h^{ab}\cdot D_c\delta h_{ab}\nonumber \\&-\,2h^{ab}D^{c}\delta h_{ac}\cdot D^{d}\delta h_{bd}+2h^{ab}D^{c}\delta h_{ab}\cdot D^{d}\delta h_{cd} \nonumber \\&+\,D^{c}\delta h^{ab}\cdot D_b\delta h_{ac}-h^{ab}\Box \delta ^2 h_{ab}\nonumber \\&-\,\frac{1}{2}h^{ab}h^{cd}D^e\delta h_{ab} \cdot D_e\delta h_{cd}+R^{(3)}_{ab}\delta ^2h^{ab} \bigg ] \nonumber \\&+\,\chi {\sqrt{h}}(h^{ab}h^{cd}-h^{ac}h^{bd})u^{e}u^{f}(\nabla _a\delta g_{eb}+\nabla _b\delta g_{ae}\nonumber \\&-\,\nabla _{e}\delta g_{ab})\cdot (\nabla _c\delta g_{fd}+\nabla _d\delta g_{cf}-\nabla _{f}\delta g_{cd}). \end{aligned}$$
(39)

Meanwhile the third term of Eq. (5) gives

$$\begin{aligned} \mathcal {G}_3&= \int _C\frac{1}{12}\chi \bigg (R^{(3)}h^{cd}+2R_{ef}h^{ce}h^{df}-Rh^{cd}\bigg )h_{cd}{\sqrt{h}}\delta h^{ab}\cdot \delta h_{ab} \nonumber \\&\quad +\, \frac{1}{4}\chi \bigg (R^{(3)}h^{ab}+2R_{cd}h^{ac}h^{bd}-R h^{ab}\bigg ){\sqrt{h}}\delta ^2 h_{ab} \nonumber \\&\quad +\, \frac{1}{8}\chi \bigg (R^{(3)}h^{ab}+2R_{ef}h^{ae}h^{bf}-R h^{ab}\bigg ){\sqrt{h}}h^{cd}\delta h_{ab}\cdot \delta h_{cd}.\nonumber \\ \end{aligned}$$
(40)

Note that some relations satisfied in a background spacetime could simplify the calculation of the fourth term, such as [5]

$$\begin{aligned} (p+\rho )h^{ab} = R^{(3)}{}^{ab}-A^{a}A^{b}-D^{b}A^{a}+h^{ab}\nabla _c A^{c}, \end{aligned}$$
(41)

which yields \(3(p+\rho )=R^{(3)}+2\nabla _aA^a\). So the fourth term of Eq. (5) becomes

$$\begin{aligned} \mathcal {G}_4= & {} \int _C -\frac{{\sqrt{h}}}{T}\frac{\delta p\delta \rho }{p+\rho } \nonumber \\= & {} \int _C -\frac{3\chi {\sqrt{h}}}{2(R^{(3)}+2\nabla _cA^c)}\delta R^{(3)}\cdot [\delta (p h_{ab})-p\delta h_{ab}]h^{ab}.\nonumber \\ \end{aligned}$$
(42)

The standard calculation yields [26]

$$\begin{aligned} {\delta }{R^{(3)}}= & {} h^{ab}{\delta }{R^{(3)}}_{ab}+{R^{(3)}}_{ab}{\delta } h^{ab} \nonumber \\= & {} D^{a}D^{b}\delta h_{ab}-h^{bc}\Box \delta h_{bc}-{R^{(3)}}{}^{ab}\delta h_{ab}. \end{aligned}$$
(43)

We have

$$\begin{aligned}&h^{ab}[\delta (p h_{ab})-p\delta h_{ab}] = h^{ab}\delta \left( R^{cd}h_{ac}h_{bd}-\frac{1}{2}Rh_{ab}\right) \nonumber \\&\qquad -\, \left( R_{cd}h^{ac}h^{bd}-\frac{1}{2}Rh^{ab}\right) \delta h_{ab} \nonumber \\&\quad = 2R^{ab}\delta h_{ab}-h^{ab}\delta R_{ab}-\frac{3}{2}\delta R-R_{cd}h^{ac}h^{bd}\delta h_{ab} \nonumber \\&\quad = 2R^{ab}\delta h_{ab}-[\delta (R_{ab}h^{ab})-R_{ab}\delta h^{ab}]\nonumber \\&\qquad -\, \frac{3}{2}\delta R-R_{cd}h^{ac}h^{bd}\delta h_{ab} \nonumber \\&\quad = R^{ab}\delta h_{ab}-\frac{1}{2}\delta (R^{(3)}-R)-\frac{5}{2}\delta R-R_{cd}h^{ac}h^{bd}\delta h_{ab} \nonumber \\&\quad = R^{ab}\delta h_{ab}{-}\frac{1}{2}D^{a}D^{b}\delta h_{ab}{+}\frac{1}{2}h^{bc}\Box \delta h_{bc}+\frac{1}{2}R^{(3)}{}^{ab}\delta h_{ab}\nonumber \\&\qquad -\,R_{cd}h^{ac}h^{bd}\delta h_{ab} \nonumber \\&\qquad -\,2\nabla ^{a}\nabla ^{b}\delta g_{ab}+2g^{bc}\nabla ^{a}\nabla _{a}\delta g_{bc}+2R^{ab}\delta g_{ab}, \end{aligned}$$
(44)

where the variation of Eq. (26) has been used. Then the fourth term can be written as

$$\begin{aligned} \mathcal {G}_4= & {} -\frac{3\chi {\sqrt{h}}}{4(R^{(3)}+2\nabla _hA^h)}(D^{a}D^{b}\delta h_{ab}-h^{ab}\Box \delta h_{ab} \nonumber \\&-\,R^{(3)}{}^{ab}\delta h_{ab}) \cdot \Big [2R^{cd}\delta h_{cd}-2R_{ef}h^{ce}h^{df}\delta h_{cd} \nonumber \\&-\,D^{c}D^{d}\delta h_{cd}+h^{cd}\Box \delta h_{cd}+R^{(3)}{}^{cd}\delta h_{cd}\nonumber \\&-\,4\nabla ^{c}\nabla ^{d}\delta g_{cd}+4g^{cd}\nabla ^{e}\nabla _{e}\delta g_{cd}+4R^{cd}\delta g_{cd}\Big ]. \end{aligned}$$
(45)

Note that the state under perturbation deviates only slightly from the equilibrium state; the \(\delta ^2h_{ab}\) terms in the expression for \(\delta ^2S\) should vanish. In fact, denoting the sum of all terms containing \(\delta ^2h_{ab}\) by \(\mathcal {G}_{\delta ^2h_{ab}}\), we have

$$\begin{aligned} \mathcal {G}_{\delta ^2 h_{ab}}= & {} \int _C\frac{\chi {\sqrt{h}}}{2}\bigg [h^{ac}h^{bd}D_c D_d\delta ^2 h_{ab} -h^{ab}\Box \delta ^2 h_{ab} \nonumber \\&+\,R^{(3)}_{ab}\delta ^2h^{ab}\bigg ]+\,\frac{\chi {\sqrt{h}}}{4}\bigg (R^{(3)}h^{ab}\nonumber \\&+\,2R_{cd}h^{ac}h^{bd}-R h^{ab}\bigg )\delta ^2 h_{ab} \nonumber \\= & {} \int _C\frac{{\sqrt{h}}}{2}\bigg \{-D_b D_a\chi +h_{ab}\Box \chi +\chi R^{(3)}_{ab}\nonumber \\&+\,\frac{\chi }{2}\bigg (-R^{(3)}h_{ab}-2R_{cd}h_a^ch_b^d+R h_{ab}\bigg ) \bigg \}\delta ^2h^{ab} \nonumber \\= & {} \int _C\frac{{\sqrt{h}}}{2}\bigg [-D_b(\chi A_a)+h_{ab}D^{c}(\chi A_c)\nonumber \\&+\,\chi R_{aeb}^lu^eu_l-\chi R_{cd}u^cu^d h_{ab}\bigg ]\delta ^2h_{ab}. \end{aligned}$$
(46)

In the last step, we use the relation \(D_a\chi =\chi A_a\) in a static background spacetime. Combining with the fact [5]

$$\begin{aligned} \nabla _au_b=-A_bu_a, \end{aligned}$$
(47)

it could be proved that

$$\begin{aligned} \mathcal {G}_{\delta ^2 h_{ab}}\equiv 0. \end{aligned}$$
(48)

Therefore, substituting Eqs. (35), (39), (40) and (45) into Eq. (5), we have

$$\begin{aligned} \delta ^2S= & {} \int _C\frac{\chi {\sqrt{h}}}{2}\Big [h^{cd}\delta h_{cd}\cdot D^aD^b\delta h_{ab}\nonumber \\&-\,h^{ab}h^{cd}\delta h_{cd}\cdot \Box \delta h_{ab}-h^{cd}R^{(3)}{}^{ab}\delta h_{ab}\delta h_{cd}\nonumber \\&+\,2\delta h^{ac}\cdot D_c D^{b}\delta h_{ab} \nonumber \\&+\,2\delta h^{bd}\cdot D^{a} D_d\delta h_{ab}-3\delta h_{ab}\cdot \Box \delta h^{ab}\nonumber \\&-\,h^{ab}D_c D_d\delta h_{ab}\cdot \delta h^{cd}-\frac{3}{2}D^{c}\delta h^{ab}\cdot D_c\delta h_{ab} \nonumber \\&-\, 2h^{ab}D^{c}\delta h_{ac}\cdot D^{d}\delta h_{bd}+2h^{ab}D^{c}\delta h_{ab}\cdot D^{d}\delta h_{cd}\nonumber \\&+\,D^{c}\delta h^{ab}\cdot D_b\delta h_{ac}-\frac{1}{2}h^{ab}h^{cd}D^e\delta h_{ab}\cdot D_e\delta h_{cd} \Big ] \nonumber \\&+\, \chi {\sqrt{h}}(h^{ab}h^{cd}-h^{ac}h^{bd})u^{e}u^{f}(\nabla _a\delta g_{eb}+\nabla _b\delta g_{ae}\nonumber \\&-\, \nabla _{e}\delta g_{ab})\cdot (\nabla _c\delta g_{fd}+\nabla _d\delta g_{cf}-\nabla _{f}\delta g_{cd}) \nonumber \\&+\, \frac{\chi {\sqrt{h}}}{12}\bigg (R^{(3)}h^{cd}+2R_{ef}h^{ce}h^{df}-Rh^{cd}\bigg )h_{cd} \cdot \delta h^{ab}\delta h_{ab} \nonumber \\&+\, \frac{\chi {\sqrt{h}}}{8}\bigg (R^{(3)}h^{ab}+2R_{ef}h^{ae}h^{bf}-R h^{ab}\bigg )h^{cd}\delta h_{ab}\cdot \delta h_{cd} \nonumber \\&-\, \frac{3\chi {\sqrt{h}}}{4(R^{(3)}+2\nabla _hA^h)}(D^aD^b\delta h_{ab}-h^{ab}\Box \delta h_{ab}\nonumber \\&-\, R^{(3)}{}^{ab}\delta h_{ab}) \cdot \Big [2R^{cd}\delta h_{cd}-2R_{ef}h^{ce}h^{df}\delta h_{cd} \nonumber \\&-\,D^cD^d\delta h_{cd}+h^{cd}\Box \delta h_{cd}+R^{(3)}{}^{cd}\delta h_{cd}\nonumber \\&-\,4\nabla ^{c}\nabla ^{d}\delta g_{cd} +4g^{cd}\nabla ^{e}\nabla _{e}\delta g_{cd}+4R^{cd}\delta g_{cd}\Big ]. \end{aligned}$$
(49)

It should be noted that for a static background spacetime, \(g_{0i}=0\) (here the index \(i=1,2,3\)). However, the system may not remain static after the perturbation, which means that the perturbation fields include \(\delta g_{00}\), \(\delta g_{0i}\) and \(\delta h_{ij}=\delta g_{ij}\). Now we decompose \(\delta g_{ab}\) into \(\delta g_{00}\), \(\delta g_{0i}\) and \(\delta h_{ij}\). The term \(-4g^{ab}\nabla ^e\nabla _e\delta g_{ab}\) in Eq. (49) can be calculated as

$$\begin{aligned} -\,4\nabla ^c\nabla ^d\delta g_{cd}= & {} -4g^{ca}\partial _a(g^{db}\nabla _b\delta g_{cd})\nonumber \\&+\,4g^{ca}\Gamma ^b_{ac}g^{de}\nabla _e\delta g_{bd} \nonumber \\= & {} -\,4g^{ca}\partial _a[g^{db}(\partial _b\delta g_{cd}-\Gamma ^e_{bc}\delta g_{ed}\nonumber \\&-\,\Gamma ^e_{bd}\delta g_{ce})]+\,4g^{ca}g^{de}\Gamma ^b_{ac}(\partial _e\delta g_{bd}\nonumber \\&-\,\Gamma ^f_{eb}\delta g_{fd}-\Gamma ^f_{ed}\delta g_{bf}). \end{aligned}$$
(50)

Then \(\delta g_{ab}\) can be decomposed into \(\delta g_{00}\), \(\delta g_{0i}\) and \(\delta h_{ij}\) (see Appendix A). The term \(4g^{cd}\nabla ^{e}\nabla _{e}\delta g_{cd}\) in Eq. (49) can be written as

$$\begin{aligned} 4g^{cd}\nabla ^{e}\nabla _{e}\delta g_{cd}= & {} 4\nabla ^e\nabla _e(g^{cd}\delta g_{cd})\nonumber \\= & {} 4\nabla ^e\nabla _e(g^{00}\delta g_{00}+h^{ij}\delta h_{ij}). \end{aligned}$$
(51)

It can be proved that \(R^{0i}=0\) in a background static spacetime, so

$$\begin{aligned} 4R^{cd}\delta g_{cd}=4R^{00}\delta g_{00}+4R^{ij}\delta h_{ij}. \end{aligned}$$
(52)

Finally, the expression for \(\delta ^2S\) can be written as

$$\begin{aligned} \delta ^2S= & {} \int _C\frac{\chi {\sqrt{h}}}{2}\Bigg [h^{cd}\delta h_{cd}\cdot D^aD^b\delta h_{ab} -h^{ab}h^{cd}\delta h_{cd}\cdot \nonumber \\&\Box \delta h_{ab}-h^{cd}R^{(3)}{}^{ab}\delta h_{ab}\delta h_{cd}+2\delta h^{ac}\cdot D_c D^{b}\delta h_{ab} \nonumber \\&+\,2\delta h^{bd}\cdot D^{a} D_d\delta h_{ab}-3\delta h_{ab}\cdot \Box \delta h^{ab}\nonumber \\&-\,h^{ab}D_c D_d\delta h_{ab}\cdot \delta h^{cd}-\frac{3}{2}D^{c}\delta h^{ab}\cdot D_c\delta h_{ab} \nonumber \\&-\,2h^{ab}D^{c}\delta h_{ac}\cdot D^{d}\delta h_{bd}+2h^{ab}D^{c}\delta h_{ab}\cdot D^{d}\delta h_{cd}\nonumber \\&+\,D^{c}\delta h^{ab}\cdot D_b\delta h_{ac}-\frac{1}{2}h^{ab}h^{cd}D^e\delta h_{ab}\cdot D_e\delta h_{cd} \Bigg ] \nonumber \\&+\,\frac{{\sqrt{h}}}{\chi } (h^{ab}h^{cd}-h^{ac}h^{bd})(\partial _a\delta g_{0b}+\partial _b\delta g_{a0}-\partial _{0}\delta h_{ab} \nonumber \\&-2\Gamma ^{i}_{ab}\delta g_{0i})(\partial _c\delta g_{0d}+\partial _d\delta g_{c0}-\partial _{0}\delta h_{cd} -2\Gamma ^{i}_{cd}\delta g_{0i}) \nonumber \\&+\,\frac{\chi {\sqrt{h}}}{12}\bigg (R^{(3)}h^{cd}+2R_{ef}h^{ce}h^{df}-Rh^{cd}\bigg )h_{cd} \nonumber \\&\cdot \,\delta h^{ab}\delta h_{ab}+\, \frac{\chi {\sqrt{h}}}{8}\bigg (R^{(3)}h^{ab}+2R_{ef}h^{ae}h^{bf}\nonumber \\&-R h^{ab}\bigg )h^{cd} \cdot \delta h_{ab}\delta h_{cd}\nonumber \\&-\, \frac{3\chi {\sqrt{h}}}{4(R^{(3)}+2\nabla _hA^h)}\nonumber \\&(D^aD^b\delta h_{ab}-h^{ab}\Box \delta h_{ab}-R^{(3)}{}^{ab}\delta h_{ab}) \nonumber \\&\cdot \Big [2R^{cd}\delta h_{cd}-2R_{ef}h^{ce}h^{df}\delta h_{cd}-D^cD^d\delta h_{cd} \nonumber \\&+\, h^{cd}\Box \delta h_{cd}+R^{(3)}{}^{cd}\delta h_{cd}+4\nabla ^{c}\nabla _{c}(g^{00}\delta g_{00})\nonumber \\&+\, 4\nabla ^{c}\nabla _{c}(h^{ij}\delta h_{ij})+4R^{00}\delta g_{00}+4R^{ij}\delta h_{ij}\nonumber \\&-\,4g^{ab}\nabla ^e\nabla _e\delta g_{ab}\Big ], \end{aligned}$$
(53)

where the last term of Eq. (53) can be written as

$$\begin{aligned}&-4g^{ab}\nabla ^e\nabla _e\delta g_{ab} \nonumber \\&\quad = -\,4(g^{00})^2\partial _0^2\delta g_{00}+4h^{ij}\partial _j(g^{00}\Gamma ^0_{0i}\delta g_{00})\nonumber \\&\qquad -\,4(g^{00})^2\Gamma ^i_{00}\Gamma ^0_{0i}\delta g_{00} -4g^{00}h^{jk}\Gamma ^i_{jk}\Gamma ^0_{0i}\delta g_{00} \nonumber \\&\qquad -\, 4g^{00}h^{ij}\partial _0\partial _j\delta g_{0i}-4h^{ij}\partial _j(g^{00}\partial _0\delta g_{0i})\nonumber \\&\qquad +\,4g^{00}h^{ij}\Gamma ^0_{j0}\partial _0\delta g_{0i} \nonumber \\&\qquad +\, 12(g^{00})^2\Gamma ^i_{00}\partial _0\delta g_{0i}+8g^{00}h^{jk}\Gamma ^i_{jk}\partial _0\delta g_{0i} \nonumber \\&\qquad -\, 4h^{ik}\partial _k(h^{jl}\partial _l\delta h_{ij})+4h^{lm}\partial _m(h^{jk}\Gamma ^i_{kl}\delta h_{ij})\nonumber \\&\qquad +\, 4h^{ik}\partial _k(g^{00}\Gamma ^j_{00}\delta h_{ij})+4h^{ik}\partial _k(h^{lm}\Gamma ^j_{lm}\delta h_{ij}) \nonumber \\&\qquad +\, 4g^{00}h^{jk}\Gamma ^i_{00}\partial _k\delta h_{ij}+4h^{lm}h^{jk}\Gamma ^i_{lm}\partial _k\delta h_{ij}\nonumber \\&\qquad -\, 4g^{00}h^{jk}\Gamma ^l_{00}\Gamma ^i_{kl}\delta h_{ij}-4h^{mn}h^{jk}\Gamma ^l_{mn}\Gamma ^i_{kl}\delta h_{ij} \nonumber \\&\qquad -\, 4(g^{00})^2\Gamma ^i_{00}\Gamma ^j_{00}\delta h_{ij}-8g^{00}h^{kl}\Gamma ^i_{00}\Gamma ^j_{kl}\delta h_{ij}\nonumber \\&\qquad -\, 4h^{kl}h^{mn}\Gamma ^i_{kl}\Gamma ^j_{mn}\delta h_{ij}. \end{aligned}$$
(54)

It is shown that whether the system is stable depends on \(\delta ^2S<0\) under the perturbation \(\delta g_{\mu \nu }\). This result gives the criterion for the thermodynamical stability for a perfect fluid star in a static background without spherical symmetry. It is worth noting that this result also gives the criterion for the cases of non-radial perturbations.

4 Summary and discussion

The maximum entropy principle suggests close relations between thermodynamics and gravity. In Ref. [5], we obtained the first variation of the total entropy for a perfect fluid in static spacetime and proved that the Einstein equations could be derived from the extrema of the total entropy and the constraint equation with some boundary conditions. That is to say, the gravitational equations could be replaced by thermodynamical relations and a constraint equation.

In this manuscript, we investigate the thermodynamical stability of an adiabatic, self-gravitating perfect fluid system deviating only slightly from equilibrium state. With thermodynamical relations, we obtain the expression of the second variation of the total entropy and the criterion for the thermodynamical stability. Specific to Einstein’s gravity with spherical symmetry spacetime and a radial perturbation, we give the explicit expression of our criterion and show that it is the same as the one in [19], which was obtained by the dynamical method. For more general cases without spherical symmetry, we transform all variations of the thermodynamical quantities to variations of geometrical quantities. Considering a perfect fluid system in a static background spacetime, our criterion could be used directly to determine whether the system is stable under any specified perturbations. Our result establishes a connection between thermodynamic and gravity under higher order variations.

Using the dynamical method, it is hard to handle the stability problems of general cases without spherical symmetry or under non-radial perturbations. However, in the framework of the thermodynamical method, the stability only depends on the signature of \(\delta ^2S\). Furthermore, if the Lagrangian for diffeomorphism invariant theories is constructed by metric and its symmetrized derivatives, the criterion for the thermodynamical stability of Eq. (5) could also be used in these modified theories, such as f(R) theories. In fact, we also proved that the thermodynamical stability is equivalent to dynamical stability in f(R) theories [27]. We found that using the thermodynamical method to obtain the stability criterion is much more direct than the dynamical method. Note that if the Lagrangian contains other scalar or vector parts, Eq. (1) need to be modified [6], which shows that the criterion for thermodynamical stability also needs to be modified.