# Extended thermodynamics of the bumblebee black holes

Zhan-Feng Mai<EMAIL_ADDRESS>Kavli Institute for Astronomy and
Astrophysics, Peking University, Beijing 100871, China  Rui Xu
<EMAIL_ADDRESS>Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Department of Astronomy, Tsinghua
University, Beijing, 100084, China  Dicong Liang  Kavli Institute for
Astronomy and Astrophysics, Peking University, Beijing 100871, China  Lijing
Shao  Kavli Institute for Astronomy and Astrophysics, Peking University,
Beijing 100871, China  National Astronomical Observatories, Chinese Academy of
Sciences, Beijing 100012, China

###### Abstract

As a vector-tensor theory including nonminimal coupling between the Ricci
tensor and a vector field, the bumblebee gravity is a potential theory to test
Lorentz symmetry violation. Recently, a new class of numerical spherical black
holes in the bumblebee theory was constructed. In this paper, we investigate
the associated local thermodynamic properties. By introducing a pair of
conjugated thermodynamic quantities $X$ and $Y$, which can be interpreted as
an extension of electric potential and charge of the Reissner-Nordström black
holes, we numerically construct a new first law of thermodynamics for
bumblebee black holes. We then study the constant-$Y$ processes in the
entropy-charge parameter space. For the constant-$Y$ processes, we also
calculate the heat capacity to study the local thermodynamic stability of the
bumblebee black holes. For a negative nonminimal coupling coefficient $\xi$,
we find both divergent and smooth phase transitions. For a positive but small
$\xi$, only a divergent phase transition is found. It turns out that there is
a critical value $0.4\kappa<\xi_{c}<0.5\kappa$ such that when
$\xi_{c}<\xi<2\kappa$, even the divergent phase transition disappears and the
bumblebee black holes thus become locally thermodynamically unstable
regardless of the bumblebee charge. As for $\xi>2\kappa$, the smooth phase
transition arises again but there no longer exists any discontinuous phase
transition for the bumblebee black holes.

##  I Introduction

The bumblebee theory is a vector-tensor theory as well as an extension of the
Einstein-Maxwell theory. In the bumblebee theory, the vector field $B_{\mu}$,
which is also called the bumblebee field, is nonminimally coupled with the
Ricci tensor quadratically and has a nonzero background value leading to
spontaneous Lorentz symmetry violation by minimizing its potential term $V$.
The action associated with the bumblebee theory is Kostelecký (2004)

$\displaystyle I$ $\displaystyle=$
$\displaystyle\int\mathrm{d}^{4}x\sqrt{-g}\Big{(}\frac{1}{2\kappa}R+\frac{\xi}{2\kappa}B^{\mu}B^{\nu}R_{\mu\nu}-\frac{1}{4}B^{\mu\nu}B_{\mu\nu}\,$
(1) (3) $\displaystyle-V(B^{\mu}B_{\mu}\pm b^{2})\Big{)}+S_{m}\,,$

where $\kappa\equiv 8\pi$, $B_{\mu}$ denotes the bumblebee field and its field
strength is defined as
$B_{\mu\nu}\equiv\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$. The coupling
constant $\xi$ indicates the strength of the quadratic interaction between the
bumblebee field and the Ricci tensor $R_{\mu\nu}$. The last term $S_{m}$
denotes the action given by normal matter. Furthermore, $V(B^{\mu}B_{\mu}\pm
b^{2})$ refers to a cosmological potential. In general, for obtaining a stable
vacuum of spacetime, the potential $V$ will be minimized if
$B^{\mu}B_{\mu}=\mp b^{2}$. In other words, in the background of
$B^{\mu}B_{\mu}=\mp b^{2}$, there exists a preferred direction of the
spacetime so that the Lorentz symmetry is broken.

If considering $V$ with a nonzero minimal value, this potential provides a
cosmological constant. However, for an unknown expression of the potential
$V$, a practical way to probe the Lorentz breaking is considering a constant
bumblebee background configuration in an asymptotically flat spacetime in the
case of vanishing minimal potential value. In the weak field limit, the
constant bumblebee field background is related to the coefficients associated
with Lorentz violation in the Standard Model extension (SME) framework
Kostelecký and Mewes (2009, 2012). In this paper, we consider the nonzero
bumblebee field background configuration, $B^{\mu}=b^{\mu}$ satisfying
$b^{\mu}b_{\mu}=\pm b^{2}$, which minimizes the potential to zero, namely

$\displaystyle\left.V(B^{\mu}B_{\mu}\pm b^{2})\right|_{B^{\mu}=b^{\mu}}$
$\displaystyle=0\,,$ (4) $\displaystyle\left.V^{\prime}(B^{\mu}B_{\mu}\pm
b^{2})\right|_{B^{\mu}=b^{\mu}}$ $\displaystyle=0\,,$ (5)

where $V^{\prime}(x)\equiv\mathrm{d}V(x)/\mathrm{d}x$.

In addition to effects in the weak-field region, effects in the strong field
region of Lorentz symmetry breaking is an intriguing topic in the bumblebee
theory as well. Black holes, as the most compact objects in our universe, have
many important and intriguing properties in both theoretical study and
astrophysical observation. Recently, several analytical solutions, including
spherical black holes as well as those with slow rotation, and their
associated properties in the bumblebee gravity have been investigated
Bertolami and Paramos (2005); Oliveira  _et al._ (2021); Güllü and Övgün
(2022); Maluf and Neves (2021); Izmailov and Nandi (2022); Liu  _et al._
(2023). Bertolami and Paramos (2005) constructed a class of Schwarzschild-like
black hole solutions analytically where the bumblebee field has only a nonzero
radial component. Furthermore, another class of spherical bumblebee black
holes, involving the bumblebee hair only with time component, has also been
numerically constructed by Xu  _et al._ (2023a), and its extreme mass ratio
inspirals (EMRIs) for probing the bumblebee hair Liang  _et al._ (2023) as
well as the possibility for Event Horizon Telescope (EHT) observation has been
studied Xu  _et al._ (2023b). However, the thermodynamics of this class of
numerical black holes, as a very important theoretical aspect of this topic,
still needs to be addressed.

Black hole thermodynamics provides an enlightening way to study various models
of gravity Hawking (1976); Bekenstein (1973); Hawking (1975); Bardeen  _et
al._ (1973); Davies (1977). Interpreting black holes as thermodynamic systems,
one can investigate black holes by only a few thermodynamic quantities.
Additionally, the thermodynamic quantities of black holes, such as energy,
temperature, entropy, and other quantities depend highly on the gravity model.
Black hole thermodynamics is also closely related to quantum gravity, and
quantum effects of gravity might become very important in the region near the
horizon where gravity is strong. Therefore, the thermodynamics of the
bumblebee black holes might provide a new angle to probe Lorentz symmetry
violation in the strong field region. In addition, the local thermodynamic
stability of the bumblebee black holes is an intriguing topic as well. It is
usually characterized by heat capacity which describes the ability of a
thermodynamic system to resist perturbations caused by exchanging a small
amount of heat with its environment.

Practically, the bumblebee theory can be interpreted as an extension of
Einstein-Maxwell theory by including a nonminimal coupling between the vector
field and gravity. From this point of view, various solutions with vector
hair, such as spherical black holes, wormholes as well as solitons in
asymptotically flat or AdS spacetimes have been constructed as well Geng and
Lu (2016); Heisenberg  _et al._ (2017a, b); Babichev  _et al._ (2017).
Moreover, the thermodynamics of these vectorized black holes in various
vector-tensor theories have been well investigated using the Wald formalism
Liu  _et al._ (2014); Fan and Lu (2015); Fan (2016, 2018). The Wald formalism,
which was developed by Wald Wald (1993); Iyer and Wald (1994), gives a
systematic procedure to analytically derive the first law of thermodynamics
for black holes in various gravity theories. Especially, the first law of
thermodynamics for the black holes in the vector-tensor theory only containing
the nonminimal coupling term $B^{\mu}B^{\mu}R_{\mu\nu}$ has been analytically
constructed by Fan (2018). By introducing the nonminimal coupling term, the
first law of thermodynamics is different from that of Reissner-Nordström (RN)
black holes. In addition to the Hawking temperature $T$ and the black hole
entropy $S$, two pairs of additional conjugated thermodynamic quantities were
introduced, even though such a vectorized black hole has only two free
parameters. The thermodynamics and associated local stability of bumblebee
black holes are thus still interesting topics.

In this paper, based on the extended numerical black hole solutions
constructed by Xu  _et al._ (2023a), we investigate the associated
thermodynamics numerically. Especially, we construct a numerical version of
the first law of thermodynamics for the bumblebee black holes. In addition to
the Hawking temperature and entropy of the bumblebee black holes, we introduce
only one pair of conjugated thermodynamic quantities, denoted as $X$ and $Y$,
where $Y$ can be interpreted as an extension of the electric charge $Q$ of RN
black holes. We further study thermodynamic processes with $Y$ being constant
and calculate the associated heat capacity $C_{Y}$.

This paper is organized as follows. In Sec. II, we briefly overview the
construction of the spherical bumblebee black holes. In Sec. III, after
calculating the Hawking temperature as well as the Wald entropy, we construct
a new version of the first law of thermodynamics for the bumblebee black
holes. As an example, we analyze thermodynamic processes where $Y$ is constant
in the entropy-charge parameter space. In Sec. IV, we study the heat capacity
$C_{Y}$ as a potential indicator for thermodynamic stability of the bumblebee
black holes and further analyze deviations between bumblebee black holes and
RN black holes. In Sec. V, we conclude this paper and give some further
comments. In this paper, we adopt the Planck-Gauss natural units;
$G=c=\hbar=k_{\rm B}=\epsilon_{0}=1$.

##  II Bumblebee black holes

Before investigating the thermodynamics of the spherical bumblebee black
holes, we first briefly overview how to construct them following the scheme
given by Xu  _et al._ (2023a). As we previously mentioned, we consider the
vacuum bumblebee theory with a nonzero bumblebee vector background
$B^{\mu}=b^{\mu}$, which minimizes the unspecified potential. With this in
mind, by performing a variation on the action in Eq. (1) with respect to the
metric $g^{\mu\nu}$ and the background bumblebee field $b_{\mu}$, we obtain
the covariant expressions for the equations of motion,

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\kappa T^{b}_{\mu\nu}\,,\quad
D_{\mu}b^{\mu\nu}=\frac{\xi}{\kappa}R^{\mu\nu}b_{\nu}\,,$ (6)

where the energy-momentum tensor contributed by the bumblebee field $b^{\mu}$
is Casana  _et al._ (2018)

$\displaystyle T^{b}_{\mu\nu}=$
$\displaystyle\frac{\xi}{2\kappa}\bigg{(}g_{\mu\nu}b^{\alpha}b^{\beta}R_{\alpha\beta}-2b_{\mu}b_{\lambda}R_{\nu}^{\phantom{\nu}\lambda}-2b_{\nu}b_{\lambda}R_{\mu}^{\phantom{\mu}\lambda}+\Box_{g}(b_{\mu}b_{\nu})-g_{{\mu\nu}}D_{\alpha}D_{\beta}(b^{\alpha}b^{\beta})+D_{\kappa}D_{\mu}\left(b^{\kappa}b_{\nu}\right)$
(7) (9)
$\displaystyle+D_{\kappa}D_{\nu}(b_{\mu}b^{\kappa})+b_{\mu\lambda}b_{\nu}^{\phantom{\nu}\lambda}-\frac{1}{4}g_{\mu\nu}b^{\alpha\beta}b_{\alpha\beta}\bigg{)}+b_{\mu\lambda}b_{\nu}{}^{\lambda}-\frac{1}{4}g_{\mu\nu}b^{\mu\nu}b_{\mu\nu}\,,$

with $b_{\mu\nu}=\partial_{\mu}b_{\nu}-\partial_{\nu}b_{\mu}$. To construct a
spherical black hole, we consider the following spherical ansatz for
$g_{\mu\nu}$ and $b_{\mu}$ in $(t,r,\theta,\varphi)$ coordinates,

$\displaystyle\mathrm{d}s^{2}=-h(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{f(r)}+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\varphi^{2}\right)\,,$
(10) (12) $\displaystyle b_{\mu}=(b_{t}(r),0,0,0).$

Here we consider the case where the bumblebee field only has the $t$ component
and it only depends on $r$ (see Ref. Xu  _et al._ (2023a) for extended
solutions with a nonzero radial component). Plugging the spherical ansatz in
Eq. (10) into the equations of motion, we have

$\displaystyle 0$ $\displaystyle=$
$\displaystyle-4h^{2}(h(-6+6f+rf^{\prime})+rf(3r\kappa
b_{t}^{2}+5f^{\prime}))+2\xi rh(2rfhb_{t}^{\prime
2}+b_{t}(b_{t}^{\prime}(rhf^{\prime}+f(4h+rh^{\prime}))+2rfh+b_{t}^{\prime\prime})$
(14) $\displaystyle-
b_{t}^{2}((4f+rf^{\prime})h^{\prime}+2rfh^{\prime\prime}))\,,$ (15)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
2h^{2}(-2+2f+rf^{\prime})-r^{2}fh^{\prime
2}+rh((6f+rf^{\prime})h^{\prime}+2rfh^{\prime\prime})\,,$ (16) (17)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle 2h(2(f-1)h^{2}+r^{2}\kappa
b_{t}fb_{t}^{\prime}h^{\prime}+rh(f(r\kappa b_{t}^{\prime
2}+2h^{\prime})-\kappa
b_{t}(b_{t}^{\prime}(4f+rf^{\prime})+2rfb_{t}^{\prime\prime})))$ (18)
$\displaystyle+r\xi
b_{t}h(-2rfb_{t}^{\prime}h^{\prime}+b_{t}((4f+rf^{\prime})h^{\prime}+2rfh^{\prime\prime}))\,.$

Here the primes denote derivatives with respect to $r$. Denoting $r_{h}$ as
the location of the horizon, in the region near the horizon, we consider the
following asymptotic expansions for $h(r),\,f(r)$ and $b_{t}(r)$,

$\displaystyle f(r)=f_{1}(r-r_{h})+f_{2}(r-r_{h})^{2}+\cdots\,,$ (19) (21)
$\displaystyle h(r)=h_{1}(r-r_{h})+h_{2}(r-r_{h})^{2}+\cdots\,,$ (23)
$\displaystyle b_{t}(r)=b_{t0}+b_{t1}(r-r_{h})+b_{t2}(r-r_{h})^{2}+\cdots.$

Then the asymptotic expansions of Eq. (14) give

$\displaystyle
f_{1}=\frac{4h_{1}}{r_{h}(4h_{1}+b_{t1}^{2}r_{h}(2\kappa-\xi))}\,,\quad
f_{2}=\frac{-16h_{1}^{2}\kappa+4b_{0}^{2}h_{1}r_{h}\kappa(2\kappa-\xi)-3b_{t1}^{2}r_{h}^{2}\xi(2\kappa+\xi)^{2}}{4h_{1}r_{h}^{2}\kappa(4h_{1}+b_{t1}^{2}r_{h}(2\kappa-\xi))}\,,\quad
b_{t0}=0\,,$ (24) (25) $\displaystyle
h_{2}=-\frac{h_{1}}{r_{h}}+\frac{1}{4}b_{t1}^{2}(2\kappa-\xi)+\frac{b_{t1}^{4}r_{h}\xi(2\kappa-\xi)^{2}}{16h_{1}\kappa}\,,\quad
b_{t2}=\frac{-16b_{t1}h_{1}^{2}\kappa+2b_{t1}^{3}h_{1}r_{h}(2\kappa-\xi)\xi+b_{t1}^{5}r_{h}^{2}\xi(-2\kappa+\xi)^{2}}{16h_{1}^{2}r_{h}\kappa}\,.$

Three coefficients $(h_{1},b_{t1},r_{h})$ are not determined by the equations
of motion. Selecting appropriate values for these coefficients and performing
numerical integration using the equations of motion to a large enough $r$, we
obtain a class of numerical black hole solutions. In addition, for an
asymptotically flat bumblebee black hole, we expect that in the asymptotic
infinity, $h(r),\,f(r)$ and $b_{t}(r)$ have the following behavior,

$\displaystyle\left.h(r)\right|_{r\to\infty}=1-\frac{2M}{r}+\frac{\tilde{h}_{2}}{r^{2}}+\cdots\,,$
(26) (28)
$\displaystyle\left.f(r)\right|_{r\to\infty}=1-\frac{2M}{r}+\frac{\tilde{f}_{2}}{r^{2}}+\cdots\,,$
$\displaystyle\left.b_{t}(r)\right|_{r\to\infty}=\mu_{\infty}-\sqrt{\frac{2}{\kappa}}\frac{Q}{r}+\frac{\tilde{b}_{2}}{r^{2}}+\cdots\,,$

where $\tilde{h}_{2},\tilde{f}_{2}$ and $\tilde{b}_{2}$ are determined by the
equations of motion. For given $(h_{1},b_{t1},r_{h})$, one may find a black
hole solution with $\left.h(r)\right|_{r\to\infty}=\tilde{h}_{0}\not=1$.
However, the following rescaling transformation

$h(r)\to\frac{h(r)}{\tilde{h}_{0}}\,,\quad
b_{t}(r)\to\frac{b_{t}(r)}{\sqrt{\tilde{h}_{0}}}\,,\quad
t\to\sqrt{\tilde{h}_{0}}\,t\,,$ (30)

keeps the equations of motion invariant, and can generate the desired black
hole solution with $\left.h(r)\right|_{r\to\infty}=1$. Furthermore, it
indicates that the bumblebee black hole has only two parameters, the same as
RN black holes. From Eq. (26), the mass $M$ and the vector charge $Q$ of the
bumblebee black hole read

$M=-\frac{1}{2}\lim_{r\to\infty}r^{2}f^{\prime}(r)\,,\quad
Q=-\sqrt{\frac{\kappa}{2}}\lim_{r\to\infty}r^{2}b_{t}^{\prime}\,.$ (31)

Additionally, $\mu_{\infty}$ is the potential at infinity and it depends on
$M$ and $Q$. When the coupling constant $\xi=0$, namely that the bumblebee
theory reduces to the Einstein-Maxwell theory, the bumblebee black holes
naturally reduce to RN black holes with

$\displaystyle f(r)=h(r)=1-\frac{2M_{\rm RN}}{r}+\frac{Q_{\rm
RN}^{2}}{r^{2}}\,,$ $\displaystyle b_{t}(r)=\mu_{\rm
RN}-\sqrt{\frac{2}{\kappa}}\frac{Q_{\rm RN}}{r}\,.$ (32)

Let us point out some essential differences between RN black holes and
bumblebee black holes. In the case of RN black holes, $\mu_{\rm
RN}=\sqrt{{2}/{\kappa}}\,{Q_{\rm RN}}/{r_{h}}$. In the Einstein-Maxwell
theory, there exists a $\text{U}(1)$ gauge symmetry. Therefore for RN black
holes, one can choose arbitrarily the zero point for the electric potential.
Usually, it is set at infinity so that $\mu_{\rm RN}=0$. The choice of
$\mu_{\rm RN}=\sqrt{{2}/{\kappa}}\,{Q_{\rm RN}}/{r_{h}}$ corresponds to
setting the zero point of the potential at the horizon. Furthermore, the
electric charge $Q_{\rm RN}$ is a conserved Noether charge associated with the
$\text{U}(1)$ symmetry in the Einstein-Maxwell theory. However, in the
bumblebee theory, due to the nonminimal coupling $\xi
b^{\mu}b^{\nu}R_{\mu\nu}$, the $\text{U}(1)$ gauge symmetry is broken.
Therefore the nonvanishing $\mu_{\infty}$ is not a gauge choice any longer,
but a source term for bumblebee black holes in addition to the fact that the
bumblebee vector charge $Q$ does not conserve when the nonminimal coupling
constant $\xi$ is nonzero.

Besides $\xi=0$, $\xi=2\kappa$ is also a special value for the coupling
constant because the solution is a Schwarzschild black hole with a nontrivial
bumblebee vector. Such kind of Schwarzschild black holes with a nontrivial
vector field has been found in other vector-tensor theories too Bertolami and
Paramos (2005); Oliveira  _et al._ (2021); Güllü and Övgün (2022). In the
bumblebee theory, they are simply

$f(r)=h(r)=1-\frac{2M_{\rm Sch}}{r}\,,\quad b_{t}(r)=\mu_{\rm
Sch}\left(1-\frac{2M_{\rm Sch}}{r}\right)\,,$ (33)

where $\mu_{\rm Sch}$ is an undetermined constant. In the Schwarzschild black
hole configuration, the bumblebee field can carry an arbitrarily large vector
charge, with the horizon of the black hole always being $r_{h}=2M_{\rm Sch}$.

For a general $\xi$, analytical black hole solutions are difficult to
construct, thus we investigate the bumblebee black holes numerically. In the
upper panel of Fig. 1, we show the relationship between the vector charge $Q$
and the radius of the event horizon for bumblebee black holes. We present
curves for seven values of $\xi$, namely
$\xi=-\kappa,\,0,\,0.1\kappa,\,\kappa,\,1.5\kappa,\,2\kappa$, and $4\kappa$ as
examples. For the case of $\xi=0$, namely the RN case, we have

$\frac{r_{h_{\rm RN}}}{M_{\rm RN}}=1+\sqrt{1-\frac{Q_{\rm RN}}{M_{\rm
RN}}}\,,$ (34)

which matches the black dashed line. For the cases of $0<\xi\lesssim
1.8\kappa$, we find that the event horizon $r_{h}/M$ is a multivalued function
with respect to $Q/M$ in certain regions, namely that there exist two black
hole solutions for the same vector charge $Q$ and mass $M$. It is clearly
indicated by the case of $\xi=\kappa$ shown as the red line in Fig. 1.
Additionally, when $\xi<2\kappa$, there exist maximal values for the vector
charge of the bumblebee black holes, namely maximal values for the charge to
mass ratio. In the case of RN black holes, the maximal charge to mass ratio is
unity, corresponding to the extremal RN black hole. For the cases of $\xi<0$,
the maximal charge to mass ratios of bumblebee black holes are less than
unity, while for those of $0<\xi<2\kappa$, the maximal charge to mass ratios
are larger than unity. For the case of $\xi=2\kappa$, the horizon is at
$r_{h}=2M$ regardless of the value of the charge, so there is no maximal
charge to mass ratio for those bumblebee black holes. As for $\xi>2\kappa$, we
find that our numerical method can generate solutions having very large charge
to mass ratios but with unacceptable errors when checking against the field
equations. So it is still unclear if the solutions end at finite charge to
mass ratios or approach all the way to infinity.

Figure 1: Upper panel: $r_{h}$ versus $Q$ for the bumblebee black holes. Lower
panel: $T$ versus $Q$ for the bumblebee black holes. When $\xi<2\kappa$,
solutions with maximal charge to mass ratios are found and indicated by solid
dots in the plots.

##  III An extended first law of thermodynamics

Black hole thermodynamics is an essential and rich topic in black hole
physics, connecting classical physics with quantum physics. In the Einstein-
Maxwell theory, an RN black hole obeys the following four laws of
thermodynamics Bardeen  _et al._ (1973):

  * •

_The zeroth law_ : The Hawking temperature $T$ is a constant on the event
horizon;

  * •

_The first law_ : $\delta M_{\rm RN}=T_{\rm RN}\delta S_{\rm RN}+\frac{Q_{\rm
RN}}{\sqrt{4\pi}r_{+}}\delta Q_{\rm RN}$;

  * •

_The second law_ : The entropy of the black hole and its environment does not
decrease with time;

  * •

_The third law_ : The Hawking temperature $T$ can not be reduced to zero by
finite physical processes.

Taking RN black holes as an example, the Hawking temperature and the
Bekenstein entropy are Hawking (1975),

$T_{\rm RN}\equiv\frac{\kappa_{+_{\rm RN}}}{2\pi}\,,\quad S_{\rm
RN}\equiv\frac{A_{\rm RN}}{4}\,.$ (35)

where $\kappa_{+_{\rm RN}}=({r_{+}-r_{-}})/{2r_{+}^{2}}$ denotes the surface
gravity for RN black holes. Additionally, the area of an RN black hole is
$A_{\rm RN}=4\pi r_{+}^{2}$, where $r_{\pm}$ denotes the outer and inner
horizon for RN black holes,

$r_{\pm}=M_{\rm RN}\pm\sqrt{M_{\rm RN}^{2}-Q_{\rm RN}^{2}}\,.$ (36)

When $M_{\rm RN}=Q_{\rm RN}$, the Hawking temperature of a RN black hole
vanishes, and we called it an extremal black hole.

In the case of bumblebee black holes found by Xu  _et al._ (2023a), the zeroth
law is naturally satisfied since they are static and spherical. We then
investigate the first law of thermodynamics. Before achieving this, we first
give the Hawking temperature $T$ and the entropy $S$ of the bumblebee black
holes. Since the bumblebee theory involves a nonminimal interaction between
the bumblebee vector and the Ricci tensor, the entropy of the bumblebee black
hole is generalized according to Wald entropy as follows Wald (1993); Iyer and
Wald (1994),

$S_{W}\equiv-2\pi\int_{r=r_{h}}\mathrm{d}^{2}x\sqrt{h}~{}\frac{\partial{\cal
L}}{\partial R^{\mu\nu\rho\sigma}}\epsilon^{\mu\nu}\epsilon^{\rho\sigma}\,,$
(37)

where the integral is on the two-sphere of the horizon, $h$ is the determinant
of the metric of the two-sphere, and $\epsilon^{\mu\nu}$ is the antisymmetric
binormal to the horizon. Applying Eq. (37) to the spherical bumblebee black
holes, because the bumblebee field $b_{t}$ vanishes on the horizon according
to the expansion coefficients in Eq. (24), one find that $b_{t}$ has no
contribution to the Wald entropy $S_{W}$. Therefore the entropy of a bumblebee
black hole is still a quarter of its area, namely

$S=\pi r_{h}^{2}\,.$ (38)

Additionally, the Hawking temperature of the bumblebee black hole is

$T\equiv\frac{\kappa_{+}}{2\pi}=\frac{1}{2\sqrt{2}\,\pi}\left.\sqrt{(\nabla^{\mu}\eta^{\nu})(\nabla_{\mu}\eta_{\nu})}\right|_{r\to
r_{h}}\,,$ (39)

where $\eta^{\mu}\equiv\left(\partial/\partial t\right)^{\mu}$ represents the
timelike Killing vector. After applying the expansions for $h(r)$ and $f(r)$
in Eq. (19), the Hawking temperature of the bumblebee black hole reads

$T=\frac{\sqrt{f_{1}h_{1}}}{4\pi}\,.$ (40)

In the lower panel of Fig. 1, we show the relation between the temperature $T$
and the vector charge $Q$ for the bumblebee black holes, for several given
values of $\xi$. The RN black holes, represented by the black dashed line
($\xi=0$), has the temperature decreasing to zero when the vector charge
approaches the maximal value, indicating the existence of the extremal RN
black hole. For the other special case, namely the Schwarzschild black holes
with $\xi=2\kappa$, the Hawking temperature is the same as that of the
Schwarzschild black hole and does not change with the vector charge.

For other values of $\xi$, we find that (i) the temperature never decreases to
zero as the vector charge grows if $\xi<\xi_{c}$ where
$0.4\kappa<\xi_{c}<0.5\kappa$, (ii) the temperature does not decrease at all
as the vector charge increases if $\xi>\xi_{c}\,$. We mentioned that for
$\xi<2\kappa$, black holes with maximal charge to mass ratios are found. But
we see that these maximally charged black holes do not have a vanishing
temperature (except for the maximally charged RN black hole), so we do not
consider them as extremal black holes. The fact that there are no extremal
black holes defined by a zero temperature for $\xi<2\kappa$ $(\xi\neq 0)$ is
unexpected, but might still be acceptable as the temperature at least is
finite. For $\xi>2\kappa$, we find that the temperature (or the factor $f_{1}$
precisely) goes to infinity when

$|b_{t1}|\to b_{c}\equiv\sqrt{\frac{4h_{1}}{r_{h}(\xi-2\kappa)}}\,,$ (41)

as shown in Eq. (24). In fact, our numerical solutions seems to indicate that
as $|b_{t1}|\to b_{c}$, the mass and the charge of the black hole go to
infinity as well. But with the increasing numerical error as $|b_{t1}|$
approaches the critical value $b_{c}$, we cannot tell for sure, let alone the
limit for the charge to mass ratio as $|b_{t1}|\to b_{c}$.111This is why we
are unable to mark a maximally charged black hole solution in Fig. 1 for
$\xi=4\kappa$.

With the temperature and the entropy found for the bumblebee black holes, we
are ready to construct the first law of thermodynamics for them. A naive
conjecture is

$\delta M=T\delta S+\mu_{\infty}\delta Q\,,$ (42)

which is a direct use of the first law for RN black holes. However, Eq. (42)
turns out to be incorrect. As we previously mentioned, the spherical bumblebee
black holes have two free parameters, which can be chosen from
$(h_{1},\,b_{t1},\,r_{h})$, or more suitably here as the entropy $S$ and the
charge $Q$. Then the mass of the black hole can be treated as a function of
the entropy $S$ and the vector charge $Q$, namely $M=M(S,Q)$, and

$\delta M=\left(\frac{\partial M}{\partial S}\right)_{Q}\delta
S+\left(\frac{\partial M}{\partial Q}\right)_{S}\delta Q\,.$ (43)

In the case of $\xi=0$ where the bumblebee black holes reduce to RN black
holes, we have

$\displaystyle M_{\rm RN}=\sqrt{\frac{S_{\rm
RN}}{4\pi}}+\sqrt{\frac{\pi}{4S_{\rm RN}}}\,Q_{\rm RN}^{2},$ $\displaystyle
T_{\rm RN}=\frac{1}{4S_{\rm RN}}\left(\sqrt{\frac{S_{\rm
RN}}{\pi}}-\sqrt{\frac{\pi}{S_{\rm RN}}}\,Q_{\rm RN}^{2}\right),$ (44)

so that one can easily check

$T_{\rm RN}=\left(\frac{\partial M_{\rm RN}}{\partial S_{\rm
RN}}\right)_{Q}\,,\quad\mu_{\rm RN}=\left(\frac{\partial M_{\rm RN}}{\partial
Q_{\rm RN}}\right)_{S}\,.$ (45)

However, we numerically find that for a general $\xi\neq 0$,

$T\not=\left(\frac{\partial M}{\partial S}\right)_{Q}\,.$ (46)

Hence, Eq. (42) does not hold in general.

An intuitive explanation for the failure of Eq. (42) is that the $\text{U}(1)$
symmetry is violated due to the nonminimal coupling term, so the vector charge
$Q$ is no longer naturally a good thermodynamic quantity to use. Following the
Wald formalism, the first law of thermodynamics for the bumblebee black holes
is actually Fan (2018)

$\displaystyle\delta M$ $\displaystyle=$ $\displaystyle T\delta
S+\mu_{\infty}\delta Q$ (47) (49)
$\displaystyle-2\xi\left(\delta(Q\mu_{\infty})-\frac{3}{2}M\delta(\mu_{\infty})^{2}-4\delta
M\mu_{\infty}^{2}\right)\,,$

where one has to introduce two pairs of conjugated thermodynamic quantities,
$(\mu_{\infty},Q)$ and $(M,\mu_{\infty}^{2})$. Of course, by treating
$\mu_{\infty}$ as a function of the two free parameters $S$ and $Q$, one is
able to put Eq. (47) solely in terms of $\delta S$ and $\delta Q$ in
principle, but the coefficients before $\delta S$ and $\delta Q$ are no longer
simply $T$ and $\mu_{\infty}$.

Figure 2: Contours of the black hole mass $M$ as well as the family of curves
of constant-$Y$ processes on the $S$-$Q$ plane. Here we take
$\xi=-\kappa,\,0,\,0.1\kappa,\,\kappa,\,1.5\kappa$, and $4\kappa$ as examples.
The yellow, orange, red and dark red lines correspond to constant-$Y$ curves
passing through points $(S_{0},Q_{0})=(2,0.2),\,(4,0.4),\,(6,0.6)$ and
$(8,0.8)$ respectively.

Our numerical solutions permit us to carry out the chain rule numerically to
find out $\delta M$ in terms of $\delta S$ and $\delta Q$. But we prefer to
find another quantity $Y$, which is a function of $S$ and $Q$ so that the
first law is in the familiar form

$\delta M=T\delta S+X\delta Y\,.$ (50)

In this sense, the quantities $X$ and $Y$ are the generalized vector potential
on the horizon and the generalized charge. For RN black holes, they reduce to

$X=\frac{Q_{\rm RN}}{\sqrt{4\pi}r_{+}}\,,\quad Y=Q_{\rm RN}\,.$ (51)

In the following of this section, we describe how $X$ and $Y$ might be
calculated and then present our numerical results of contour plots for
constant-$Y$ processes in the $S$-$Q$ parameter space.

As $Y(S,Q)$ is a function of $S$ and $Q$, we have

$\delta Y=X_{1}\delta S+X_{2}\delta Q\,,$ (52)

so that Eq. (50) is

$\delta M=\left(T+XX_{1}\right)\delta S+XX_{2}\delta Q\,.$ (53)

The unknowns $X,\,X_{1},\,X_{2}$ and $Y$, all treating as functions of $S$ and
$Q$, therefore satisfy

$\displaystyle\left(\frac{\partial M}{\partial S}\right)_{Q}$ $\displaystyle=$
$\displaystyle T+XX_{1}\,,$ (54) (55) $\displaystyle\left(\frac{\partial
M}{\partial Q}\right)_{S}$ $\displaystyle=$ $\displaystyle XX_{2}\,,$ (56)
(57) $\displaystyle\frac{\partial^{2}M}{\partial Q\partial S}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
Q}\bigg{|}_{S}\left(T+XX_{1}\right)=\frac{\partial}{\partial
S}\bigg{|}_{Q}\left(XX_{2}\right),$ (58) (59)
$\displaystyle\frac{\partial^{2}Y}{\partial Q\partial S}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial X_{1}}{\partial
Q}\right)_{S}=\left(\frac{\partial X_{2}}{\partial S}\right)_{Q}.$ (60)

Eliminating $X_{1}$ and $X_{2}$, we find a partial differential equation (PDE)
for $X(S,Q)$,

$\left(\frac{\partial M}{\partial
Q}\right)_{S}X_{S}+\left(T-\left(\frac{\partial M}{\partial
S}\right)_{Q}\right)X_{Q}=\left(\frac{\partial T}{\partial Q}\right)_{S}X\,,$
(61)

where $X_{S}=\left(\partial X/\partial S\right)_{Q}$ and $X_{Q}=\left(\partial
X/\partial Q\right)_{S}$. Our numerical solutions give the functions
$M(S,Q),\,T(S,Q)$ and their partial derivatives. Together with a reasonable
boundary condition, $X$ can be solved from Eq. (61). Once $X=X(S,Q)$ is known,
we have

$\displaystyle\left(\frac{\partial Y}{\partial
S}\right)_{Q}=X_{1}=\frac{1}{X}\left(\left(\frac{\partial M}{\partial
S}\right)_{Q}-T\right),$ $\displaystyle\left(\frac{\partial Y}{\partial
Q}\right)_{S}=X_{2}=\frac{1}{X}\left(\frac{\partial M}{\partial
Q}\right)_{S}.$ (62)

So $Y$ as a function of $S$ and $Q$ can be obtained by integrating Eq. (52)
given an initial value $Y(S_{0},Q_{0})$.

Let us address that the quantities $X$ and $Y$ thus calculated are not unique.
It is clear in Eq. (50) that if one pair of $X$ and $Y$ have been found, then
$\tilde{Y}(Y)$ which is any reasonable function of $Y$ together with
$\tilde{X}=XdY/d\tilde{Y}$ accomplishes the first law too, namely

$\displaystyle\delta M=T\delta S+\tilde{X}\delta\tilde{Y}.$ (63)

This freedom of transforming $Y$ to $\tilde{Y}(Y)$ is equivalent to the
arbitrariness in specifying the boundary condition for $X$ in solving Eq.
(61), leading to numerous choices for $X$ and $Y$. Among all the choices, one
might be able to fix $X$ and $Y$ up to a constant factor if the Smarr relation
Smarr (1973)

$\displaystyle M=2TS+XY$ (64)

is required. However, the Smarr relation is not guaranteed to be consistent
with Eq. (61). It turns out to require the following condition (see Appendix A
for the derivation)

$\displaystyle 0$ $\displaystyle=$ $\displaystyle 2S\left(\left(\frac{\partial
M}{\partial Q}\right)_{S}\left(\frac{\partial T}{\partial
S}\right)_{Q}-\left(\frac{\partial M}{\partial
S}\right)_{Q}\left(\frac{\partial T}{\partial Q}\right)_{S}\right)$ (65)
$\displaystyle+T\left(\frac{\partial M}{\partial
Q}\right)_{S}+M\left(\frac{\partial T}{\partial Q}\right)_{S}.$

One can check that RN black holes satisfy the condition using the analytical
relation in Eq. (44). But for numerical bumblebee black hole solutions, we
cannot say for sure if the condition is satisfied as small residuals are found
in our numerical calculation.

Due to the lack of a unique choice for $X$ and $Y$, we quit calculating them
but focus on investigating processes with $Y$ (thus also any $\tilde{Y}(Y)$)
being constant. These processes are curves on the $S$-$Q$ plane defined by

$\left(\frac{\mathrm{d}Q}{\mathrm{d}S}\right)_{Y}=-\frac{X_{1}}{X_{2}}=\frac{T-\left(\frac{\partial
M}{\partial S}\right)_{Q}}{\left(\frac{\partial M}{\partial Q}\right)_{S}}\,.$
(66)

In Fig. 2, the colored regions show contours of the black hole mass $M$, while
the curves are solutions to Eq. (66), namely isolines of the generalized
charge $Y$ in the $S$-$Q$ plane. Note that black hole solutions do not exist
in the blank zones in the plots. To numerically obtain the curves representing
the constant-$Y$ processes in Fig. 2, we start with an initial point
$(S_{0},Q_{0})$ on the plane and proceed to the next point by the following
recursive relation

$Q_{i+1}=Q_{i}+\frac{T_{i}-\left(\frac{\bar{M}_{i}-M_{i}}{\bar{S}_{i}-S_{i}}\right)_{Q_{i}}}{\left(\frac{\tilde{M}_{i}-M_{i}}{\tilde{Q}_{i}-Q_{i}}\right)_{S_{i}}}(S_{i+1}-S_{i})\,.$
(67)

Here $\bar{M}_{i}$ and $\bar{S}_{i}$ denote the mass and entropy of black
holes sharing the same vector charge $Q_{i}$ with those of mass $M_{i}$ and
entropy $S_{i}$. Similarly, $\tilde{M}_{i}$ and $\tilde{Q}_{i}$ denote the
mass and vector charge of those black holes sharing the same entropy $S_{i}$
with those of mass $M_{i}$ and charge $Q_{i}$. To ensure numerical accuracy,
we choose the step $S_{i+1}-S_{i}\approx 10^{-2}$ (in the Planck-Gauss natural
units). Following this scheme step-by-step, we obtain a series of points
revealing a constant-$Y$ process on the $S$-$Q$ plane. Varying the starting
point, we can obtain a family of constant-$Y$ processes for different values
of $Y$. As shown in Fig. 2, we choose
$(S_{0},Q_{0})=(2,0.2),(4,0.4),(6,0.6),(8,0.8)$ as four initial points, and
calculated four typical constant-$Y$ curves in the $S$-$Q$ plane for several
values of $\xi$.

We notice that in the case of $\xi=0$ (RN black hole), the curves of
constant-$Y$ processes are horizontal straight lines, consistent with the fact
$Y=Q_{\rm RN}$ for RN black holes. For other values of $\xi$, we find that the
lines of constant-$Y$ processes are not parallel anymore. The larger $|\xi|$
is, the more curved the lines are. For example, when $\xi=4\kappa$, the
constant-$Y$ lines are so bent that for a given entropy $S$, there are
generally two values of $Q$ corresponding to the same value of $Y$, as shown
in the third plot in the bottom panel of Fig. 2.

##  IV local thermodynamic stability

After revealing constant-$Y$ processes in the $S$-$Q$ plane, now we focus on
the local thermodynamic stability which is indicated by the heat capacity of a
thermodynamic system in its thermodynamic equilibrium Avramov  _et al._
(2023). In general, heat capacity is defined as the heat transferred per unit
change of temperature,

$C\equiv\frac{\delta{\cal Q}}{\delta T}\,.$ (68)

Heat capacity depends on the specific thermodynamic process under
consideration. Different thermodynamic processes give different heat
capacities. For example, two commonly used heat capacities for an ideal gas
system are under constant volume (isochoric) and under constant pressure
(isobaric) correspondingly. A positive heat capacity, $C>0$, indicates that
the system in equilibrium is locally thermodynamically stable, while a
negative heat capacity, $C<0$, indicates that it is locally thermodynamically
unstable. In general, most matters, such as ideal gas, liquid, and solid have
positive heat capacities. Supposing that a thermodynamic system composed of
ideal gas is in thermodynamic equilibrium with its environment, when the
temperature of the environment raises, the system of ideal gas will absorb
heat from the environment. Since the heat capacity of ideal gas is positive,
through absorbing heat from the environment, the temperature of the ideal gas
will increase and finally reach a new equilibrium with the environment. Thus
the ideal gas system is locally thermodynamically stable.

However, self-gravitating systems are different. In black hole thermodynamics,
black holes can have negative heat capacities, implying that they are locally
thermodynamically unstable. For example, the heat capacity of the
Schwarzschild black hole is

$C_{\rm Sch}=\frac{\delta{\cal Q}}{\delta
T}=\frac{\mathrm{d}M}{\mathrm{d}T}=-8\pi M^{2}\,.$ (69)

The negative heat capacity means that a Schwarzschild black hole is locally
thermodynamically unstable, reflecting the fact that a Schwarzschild black
hole will continuously release energy through Hawking radiation, the so-called
black hole evaporation Traschen (1999). As for the RN black holes, considering
constant-$Q$ processes which are analogs of constant-volume processes for
ideal gas, the associated heat capacity is Altamirano  _et al._ (2014)

$\displaystyle C_{Q}$ $\displaystyle=$
$\displaystyle\left(\frac{\mathrm{d}M}{\mathrm{d}T}\right)_{Q}$
$\displaystyle=$
$\displaystyle-\frac{2M^{2}\pi(1+\sqrt{1-\lambda^{2}})^{2}(1-\lambda^{2}+\sqrt{1-\lambda^{2}})}{1-2\lambda^{2}+\sqrt{1-\lambda^{2}}}\,,$

where $\lambda\equiv|Q_{\rm RN}|/M_{\rm RN}$. It is easy to observe that there
is a discontinuous point when $1-2\lambda^{2}+\sqrt{1-\lambda^{2}}=0$, giving
$\lambda=\sqrt{3}/2$. When $\lambda<\sqrt{3}/2$, the heat capacity of an RN
black hole is negative, giving a locally thermodynamically unstable charged
black hole, while when $\lambda>\sqrt{3}/2$, the heat capacity turns to be
positive through the discontinuous point, giving a locally thermodynamically
stable charged black hole. Additionally, such a change of the sign for heat
capacity due to a divergence indicates a discontinuous phase transition
Avramov  _et al._ (2023).

Figure 3: Relation between the heat capacity associated with constant-$Y$
processes and the charge for the bumblebee black holes.

Based on the constant-$Y$ processes we discussed in the previous section, Eq.
(IV) can be extended for the bumblebee black holes as $Y$ is the generalized
charge appearing in the first law of thermodynamics. The heat capacity for
constant-$Y$ processes can be calculated as

$C_{Y}=\left(\frac{\mathrm{d}M}{\mathrm{d}T}\right)_{Y}\,.$ (72)

Since $C_{Y}$ has no known analytic expression, we numerically calculate
$C_{Y}$ for various $\xi$ and present the results as functions of the charge
in Fig. 3.

First, Fig. 3 indicates that for each $\xi$, when $|Q|\to 0$, the heat
capacity of the bumblebee black holes becomes
$C_{Y}/M^{2}\to-8\pi\approx-25.12$, which is consistent with the fact that the
bumblebee black holes reduce to the Schwarzschild black holes when the vector
charge vanishes. We are not able to show in Fig. 3 the value of $C_{Y}/M^{2}$
at $Q=0$ as the numerical error prevents our calculation when $|Q|/M$ is very
small. But the trend for $C_{Y}/M^{2}\to-8\pi\approx-25.12$ is clear.

Second, when $\xi<0$ (taking $\xi=-\kappa$ as an example), we find that
similar to RN black holes, there exist a divergent point where the sign of the
heat capacity changes. A closer inspection finds that the charge to mass ratio
at the discontinuity is smaller than that of RN black holes ($Q/M\approx 0.52$
for $\xi=-\kappa$). In addition to the discontinuous phase transition, when
$Q/M$ is large enough, we find that there exists a continuous phase transition
point where the heat capacity changes the sign smoothly. Therefore, for
$\xi<0$, a bumblebee black hole begins as a locally thermodynamically unstable
black hole, turns to stable through a discontinuous phase transition, and then
turns to unstable again through a smooth phase transition as the vector charge
grows.

Third, when $\xi>0$, we find that the divergence of the heat capacity, if
exists, appears on the right of the divergence point of the RN black holes. As
$\xi$ increases, we find that the divergence point disappears at
$\xi\to\xi_{c}$, where $0.4\kappa<\xi_{c}<0.5\kappa$ is the same critical
value found in Sec. III. When $\xi_{c}<\xi<2\kappa$ (for example, $\xi=\kappa$
and $1.5\kappa$ in Fig. 3), all bumblebee black holes are locally
thermodynamically unstable and do not have any phase transition point.

As for the case of $\xi=2\kappa$ where the bumblebee black hole reduces to the
Schwarzschild black hole, Eq. (33) the generalized charge $Y$ is undefined.
There is only one definition of heat capacity, which is Eq. (69).

Finally, when $\xi>2\kappa$ (taking $\xi=4\kappa$ as an example), we find that
there exists only a continuous phase transition point ($Q/M\approx 0.82$ for
$\xi=4\kappa$). The bumblebee black hole is locally thermodynamically unstable
when the charge to mass ratio is small and is stable when the charge to mass
ratio is large.

##  V conclusion and discussion

The bumblebee theory is not only one of the modified gravity theories with
Lorentz symmetry violation but also a potential extension of the Einstein-
Maxwell theory. Similar to RN black holes, a spherical black hole with a
nontrivial bumblebee field can also be described by its mass $M$ and charge
$Q$. We have numerically constructed the first law of thermodynamics for the
bumblebee black holes by introducing a pair of new thermodynamic quantities
$X$ and $Y$, which can be interpreted as generalizations of the electric
potential $\mu_{\rm RN}$ and charge $Q_{\rm RN}$ of RN black holes
respectively. Algorithms to calculate $X$ and $Y$ in terms of the entropy $S$
and the charge $Q$ have been presented and constant-$Y$ processes in the
$S$-$Q$ plane have been constructed.

To investigate the local thermodynamic stability of the bumblebee black holes,
we have further calculated the heat capacity associated with the constant-$Y$
processes. When $0<\xi\lesssim 0.4\kappa$, we find that there exists a
divergence with the changing of the sign for the heat capacity, indicating a
discontinuous phase transition from locally thermodynamically unstable black
holes to stable black holes. When $\xi<0$, in addition to the discontinuous
phase transition, we also find that when the vector charge is sufficiently
large, there exists of a continuous phase transition where the bumblebee black
hole becomes locally thermodynamically stable again. As for the case of
$0.4\kappa\lesssim\xi<2\kappa$, the divergence of heat capacity disappears and
the bumblebee black holes are locally thermodynamically unstable for all
charge to mass ratios. When $\xi>2\kappa$, there is no divergence but a smooth
phase transition point that separates locally thermodynamically stable and
unstable bumblebee black holes.

In modified gravity, the Wald formalism is usually used to analytically
investigate black hole thermodynamics. However, the first law of
thermodynamics given by the Wald formalism needs to introduce additional
conjugated thermodynamic quantities, usually more than the number of the free
parameters in the black hole solutions. Therefore, we try another path that
numerically constructs the first law of thermodynamics, using the bumblebee
theory as an example. Its correspondence to the Wald formalism still needs
further investigation.

In addition, we remark that we here focus on the local thermodynamic
properties of the bumblebee black holes, but the global stability is still
unclear. The global and local thermodynamic properties are not equivalent in
black hole thermodynamics since black holes generally have negative heat
capacity. In principle, recall that $\left.b_{t}\right|_{r\to\infty}\not=0$ in
Eq. (26), and the bumblebee vector field can be treated as sourced at
infinity, namely that our bumblebee black holes are globally grand canonical
ensembles. To investigate the global thermodynamic stability of the bumblebee
black holes, one needs to calculate the free energy by fixing $T$ and
$\mu_{\infty}$ Mai and Yang (2021). Furthermore, as also mentioned in Ref. Mai
and Yang (2021), black holes in astrophysical environments, usually with
accretion disk and interaction with other matters, are globally grand
canonical ensembles. The global thermodynamic stability should be considered.
The astrophysical properties of the bumblebee black holes have been studied in
Ref. Xu  _et al._ (2023b). In this paper, we mainly focus on the theoretical
properties of the bumblebee black hole thermodynamics. The application of the
thermodynamics studied in the current work to astrophysical black holes is
unclear at the current stage. We thus leave the application aspect of our work
to future studies.

###### Acknowledgements.

We thank Run-Qiu Yang and Hong Lü for useful discussions. This work was
supported by the National Natural Science Foundation of China (Grants No.
12247128, N0. 11991053, No. 11975027), the China Postdoctoral Science
Foundation (No. 2021TQ0018), the National SKA Program of China (No.
2020SKA0120300), the Max Planck Partner Group Program funded by the Max Planck
Society, and the High-Performance Computing Platform of Peking University. R.
X. is supported by the Boya postdoctoral fellowship at Peking University.

##  Appendix A A necessary condition for the Smarr relation

Here we derive Eq. (65) under the assumption of the Smarr relation in Eq. (64)
with $X$ and $Y$ satisfying the first law in Eq. (50). Using the Smarr
relation, the derivative of $X$ with respect to $S$ along the constant-$Y$
processes defined by Eq. (66) is

$\displaystyle Y\left(\frac{\mathrm{d}X}{\mathrm{d}S}\right)_{Y}$
$\displaystyle=$
$\displaystyle\left(\frac{\mathrm{d}}{\mathrm{d}S}\left(M-2TS\right)\right)_{Y}$
(73) $\displaystyle=$ $\displaystyle\left(\frac{\partial M}{\partial
S}\right)_{Q}-2S\left(\frac{\partial T}{\partial S}\right)_{Q}-2T$
$\displaystyle+\left(\left(\frac{\partial M}{\partial
Q}\right)_{S}-2S\left(\frac{\partial T}{\partial
Q}\right)_{S}\right)\left(\frac{\mathrm{d}Q}{\mathrm{d}S}\right)_{Y}\,.$

In the meanwhile, $X$ satisfies Eq. (61), giving

$\displaystyle\left(\frac{\mathrm{d}X}{\mathrm{d}S}\right)_{Y}$
$\displaystyle=$ $\displaystyle\left(\frac{\partial X}{\partial
S}\right)_{Q}+\left(\frac{\partial X}{\partial
Q}\right)_{S}\left(\frac{\mathrm{d}Q}{\mathrm{d}S}\right)_{Y}$ (74)
$\displaystyle=$ $\displaystyle\frac{X\left(\frac{\partial T}{\partial
Q}\right)_{S}}{\left(\frac{\partial M}{\partial Q}\right)_{S}}\,,$

where Eq. (66) has been used. For Eq. (73) and Eq. (74) to be consistent with
each other, we find

$\displaystyle\left(\frac{\partial M}{\partial
Q}\right)_{S}\left(\left(\frac{\partial M}{\partial
S}\right)_{Q}-2S\left(\frac{\partial T}{\partial S}\right)_{Q}-2T\right)$ (75)
$\displaystyle+\left(\left(\frac{\partial M}{\partial
Q}\right)_{S}-2S\left(\frac{\partial T}{\partial
Q}\right)_{S}\right)\left(T-\left(\frac{\partial M}{\partial
S}\right)_{Q}\right)$ $\displaystyle=$
$\displaystyle\left(M-2TS\right)\left(\frac{\partial T}{\partial
Q}\right)_{S}\,,$

which simplifies to Eq. (65).

## References

  * Kostelecký (2004) V. A. Kostelecký, Phys. Rev. D 69, 105009 (2004), arXiv:hep-th/0312310 .
  * Kostelecký and Mewes (2009) V. A. Kostelecký and M. Mewes, Phys. Rev. D 80, 015020 (2009), arXiv:0905.0031 [hep-ph] .
  * Kostelecký and Mewes (2012) A. Kostelecký and M. Mewes, Phys. Rev. D 85, 096005 (2012), arXiv:1112.6395 [hep-ph] .
  * Bertolami and Paramos (2005) O. Bertolami and J. Paramos, Phys. Rev. D 72, 044001 (2005), arXiv:hep-th/0504215 .
  * Oliveira  _et al._ (2021) R. Oliveira, D. M. Dantas, and C. A. S. Almeida, EPL 135, 10003 (2021), arXiv:2105.07956 [gr-qc] .
  * Güllü and Övgün (2022) I. Güllü and A. Övgün, Annals Phys. 436, 168721 (2022), arXiv:2012.02611 [gr-qc] .
  * Maluf and Neves (2021) R. V. Maluf and J. C. S. Neves, Phys. Rev. D 103, 044002 (2021), arXiv:2011.12841 [gr-qc] .
  * Izmailov and Nandi (2022) R. N. Izmailov and K. K. Nandi, Class. Quant. Grav. 39, 215006 (2022).
  * Liu  _et al._ (2023) W. Liu, X. Fang, J. Jing, and J. Wang, Eur. Phys. J. C 83, 83 (2023), arXiv:2211.03156 [gr-qc] .
  * Xu  _et al._ (2023a) R. Xu, D. Liang, and L. Shao, Phys. Rev. D 107, 024011 (2023a), arXiv:2209.02209 [gr-qc] .
  * Liang  _et al._ (2023) D. Liang, R. Xu, Z.-F. Mai, and L. Shao, Phys. Rev. D 107, 044053 (2023), arXiv:2212.09346 [gr-qc] .
  * Xu  _et al._ (2023b) R. Xu, D. Liang, and L. Shao, Astrophys. J. 945, 148 (2023b), arXiv:2302.05671 [gr-qc] .
  * Hawking (1976) S. W. Hawking, Phys. Rev. D 13, 191 (1976).
  * Bekenstein (1973) J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * Hawking (1975) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Bardeen  _et al._ (1973) J. M. Bardeen, B. Carter, and S. W. Hawking, Commun. Math. Phys. 31, 161 (1973).
  * Davies (1977) P. C. W. Davies, Proc. Roy. Soc. Lond. A 353, 499 (1977).
  * Geng and Lu (2016) W.-J. Geng and H. Lu, Phys. Rev. D 93, 044035 (2016), arXiv:1511.03681 [hep-th] .
  * Heisenberg  _et al._ (2017a) L. Heisenberg, R. Kase, M. Minamitsuji, and S. Tsujikawa, JCAP 08, 024 (2017a), arXiv:1706.05115 [gr-qc] .
  * Heisenberg  _et al._ (2017b) L. Heisenberg, R. Kase, M. Minamitsuji, and S. Tsujikawa, Phys. Rev. D 96, 084049 (2017b), arXiv:1705.09662 [gr-qc] .
  * Babichev  _et al._ (2017) E. Babichev, C. Charmousis, and A. Lehébel, JCAP 04, 027 (2017), arXiv:1702.01938 [gr-qc] .
  * Liu  _et al._ (2014) H.-S. Liu, H. Lü, and C. N. Pope, JHEP 06, 109 (2014), arXiv:1402.5153 [hep-th] .
  * Fan and Lu (2015) Z.-Y. Fan and H. Lu, Phys. Rev. D 91, 064009 (2015), arXiv:1501.00006 [hep-th] .
  * Fan (2016) Z.-Y. Fan, JHEP 09, 039 (2016), arXiv:1606.00684 [hep-th] .
  * Fan (2018) Z.-Y. Fan, Eur. Phys. J. C 78, 65 (2018), arXiv:1709.04392 [hep-th] .
  * Wald (1993) R. M. Wald, Phys. Rev. D 48, R3427 (1993), arXiv:gr-qc/9307038 .
  * Iyer and Wald (1994) V. Iyer and R. M. Wald, Phys. Rev. D 50, 846 (1994), arXiv:gr-qc/9403028 .
  * Casana  _et al._ (2018) R. Casana, A. Cavalcante, F. P. Poulis, and E. B. Santos, Phys. Rev. D 97, 104001 (2018), arXiv:1711.02273 [gr-qc] .
  * Smarr (1973) L. Smarr, Phys. Rev. Lett. 30, 71 (1973), [Erratum: Phys.Rev.Lett. 30, 521–521 (1973)].
  * Avramov  _et al._ (2023) V. Avramov, H. Dimov, M. Radomirov, R. C. Rashkov, and T. Vetsov, (2023), arXiv:2302.11998 [gr-qc] .
  * Traschen (1999) J. H. Traschen, in  _1999 Londrona Winter School on Mathematical Methods in Physics_ (1999) arXiv:gr-qc/0010055 .
  * Altamirano  _et al._ (2014) N. Altamirano, D. Kubiznak, R. B. Mann, and Z. Sherkatghanad, Galaxies 2, 89 (2014), arXiv:1401.2586 [hep-th] .
  * Mai and Yang (2021) Z.-F. Mai and R.-Q. Yang, Phys. Rev. D 104, 044008 (2021), arXiv:2101.00026 [gr-qc] .