\section{Introduction} 
\label{Intro} 

Since Einstein submitted papers on general relativity \cite{Einstein:1916}, classical or quantum field theory on a curved spacetime has extensively been investigated. When spacetime is curved, the physical quantities defined on flat spacetime are required to be modified suitably in accordance with general covariance. 
For example, a conserved current, which exists in the presence of global symmetry in the system \cite{Nother:1918}, is modified to be a covariantly conserved one on a general curved spacetime. 

However there has been no general argument to define a conserved charge from a covariantly conserved current, which inevitably causes a problem to define energy and momentum. 
Einstein originally argued that the conservation law of energy and momentum for matter follows as long as they are combined with those for gravitational field \cite{Einstein:1916}. (See also \cite{Landau:1975}.)
The corresponding energy momentum tensor for the gravitational field, however, is not covariant under general coordinate transformation. As a result an energy defined as in the case of flat spacetime depends on a coordinate system and conserves only in the particular frame. 

One way to circumvent this issue is to define an energy locally on the asymptotic region of spacetime called quasi-local energy. 
Initially the quasi-local energy and momentum were studied on an asymptotically flat spacetime by recasting gravity system into the Hamiltonian dynamics known as the ADM formalism \cite{Arnowitt:1962hi}. (See also \cite{Bondi:1962px}.) They are defined by a surface integral in the asymptotic region, by which the invariance under a class of general coordinate transformations preserving a boundary condition was achieved. 
This result has been further extended for a more general curved spacetime with surface terms suitably incorporated 
\cite{Brown:1992br,Hawking:1995fd,Horowitz:1998ha,Balasubramanian:1999re,Ashtekar:1999jx}. 
A caveat in this extension is that boundary terms accompany with divergence even in the flat spacetime, so that one needs to subtract it by comparing a reference frame or by adding local counter terms.

The authors of the present paper investigated a property of a black hole holographically realized by the flow equation method \cite{Aoki:2020ztd}. 
In the study we encountered a situation to evaluate the energy of the total system with matter spread all over the space, which is required to be evaluated by the volume integral of the energy density.
We  reached a manifestly covariant  definition of a conserved charge from a covariantly conserved current  in a general curved spacetime, which 
extends the one in Ref.~\cite{Fock1959},  and
improves the one given in \cite{Komar:1962xha,Abbott:1981ff} for special backgrounds. 
This allows us to define energy and momentum for matter in a form of the volume integral at an arbitrary time slice of a given curved spacetime.
A virtue is to enable one to evaluate charges of black holes just like an electric charge of the electron in electro-magnetism by the volume integral of 
an integrable singularity such as the delta function.
Furthermore, applying
our definition to the 
energy
of a compact star,
we
discover a correction to the mass formula obtained from the Oppenheimer-Volkoff equation, 
which represents a contribution from the gravitational interaction and becomes
68\% of the mass 
at the maximum
for a constant density.

\section{ Conserved charge from covariantly conserved current}
\label{sec:Proposal}

Consider any classical or quantum field theory on a general curved spacetime.  
Suppose there exists a covariantly conserved current $J^\mu$, $\nabla_\mu J^\mu = 0$, where $\nabla_\mu$ is the covariant derivative for the metric $g_{\mu\nu}$. 
Then we claim that the following quantity is conserved under the given time evolution
\begin{equation}
Q(t) := \int_{M_t} d^{d-1}\vec x\,  \sqrt{\vert g\vert} J^0 (t,\vec x), 
\label{eq:def_q}
\end{equation}
where $M_t$ represents a time slice of the spacetime $M$ at the time $t$, $g$ denotes the determinant of $g_{\mu\nu}$, and $d$ is the dimension of the spacetime $M$. If there exists boundary for $M_t$, we set the boundary condition for
the fields to fall off sufficiently fast at boundary of $M_t$ for all $t$. We emphasize that $g$ is the determinant of the metric in the total spacetime, which contains the time components. 

To show this, we assume the spacetime has the foliation structure for simplicity. 
(The same argument is used in literature. For example, see \cite{Townsend:1997ku}.)
Let us consider the same quantity defined by \eqref{eq:def_q} at another time slice with $t'$ greater than $t$, and take a submanifold $M'$ with the foliation structure whose boundary contains $M_t$ and $M_{t'}$. Such a manifold may be written formally as $M' = \prod_{t\leq s \leq t'} M_{s}$.
Under the boundary condition, the difference between $Q(t^\prime)$ and $Q(t)$ becomes
\beal{ 
Q(t') - Q(t) 
=& \int_{M'} d^{d} x\, \partial_\mu( \sqrt{\vert g\vert} J^\mu (t,\vec x) ) = 0,
}
where we used
$\partial_\mu( \sqrt{\vert g\vert} J^\mu (t,\vec x) ) =\sqrt{\vert g\vert} \nabla_\mu J^\mu (t,\vec x) =0$. 
This proves that $Q(t)$ is independent of $t$. 
The charge $Q$ is a scalar under the assumption, though the generalization is straightforward. Note  that $Q$ is not a scalar if it is defined from a higher rank tensor rather than the vector.

This formula can be applied to the computation of a conserved charge for any gravitational systems with a Killing vector. A covariantly conserved current associated with a Killing vector $\xi$ can be constructed as
\begin{equation} 
J^\mu = T^\mu{}_\nu \xi^\nu, 
\label{currentKilling} 
\end{equation}
where $T^\mu{}_\nu$ is the given energy momentum tensor for matter. It can be easily shown that this is covariantly conserved by using $\nabla_\mu T^\mu{}_\nu =0$ and $\nabla_{\mu} \xi_{\nu}+\nabla_{\nu} \xi_{\mu} = 0$. 
Note that the definition of a conserved  charge using eq.~\eqref{currentKilling} in $d=4$ appears in \cite{Fock1959,Trautman:2002zz}, but the formula has been rarely used in literature as far as we know.
If $\xi^\mu$ is a Killing vector associated with the time translation, the conserved charge 
becomes the total energy of the system
\begin{eqnarray}
E = \int_{M_t} d^{d-1}\vec x \sqrt{\vert g\vert} T^0{}_\mu \xi^\mu,
\end{eqnarray}
which 
agrees with the standard definition of the energy
in the flat 
background with $\xi^\mu=-\delta^\mu_0$.
In the next section we compute charges of several black holes by using this formula.

\section{Conserved charges of black holes} 
\label{sec:Blackhole}

In this section we compute a conserved charge for Schwarzschild and BTZ black holes employing the presented formula. 

\subsection{ Schwarzschild black hole} 


In order to explain the key idea of the calculation we start with the simplest setup. That is, 
we begin with the Einstein equation 
\begin{equation} 
R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R + \Lambda g_{\mu\nu} = 0. 
\label{EOMV}
\end{equation}
This is satisfied by the Schwarzschild black hole solution: 
\beal{
d s^2 =& -f(r) (dx^0)^2  + {1\over f(r) } dr^2  + r^2\tilde g_{ij} dx^i dx^j,
\label{eq:SBH}
}
where $r$ is the radial coordinate and the $d-2$ dimensional manifold fibered over the cone is an Einstein manifold, whose Ricci tensor is given by ${}^{(d-2)} R_{ij} = (d-3) k \tilde g_{ij}$ with a constant $k$, and 
\begin{equation}
f(r) =\frac{-2 \Lambda r^2}{(d-2)(d-1)} + k - \frac{2G_N M }{r^{d-3}}. 
\label{f}
\end{equation}
Note that for a positive or non-positive $k$ the submanifold is compact or non-compact, respectively.

Since this is a static solution, there exists a Killing vector with $\xi^\mu=-\delta^\mu_0$, which corresponds to the time translation. 
Thus the corresponding charge is the energy of the system:
\begin{equation} 
E= \int d^{d-1}\vec x\,  \sqrt{\vert g\vert} (-T^0{}_0), 
\label{mass}
\end{equation}
where the matter energy momentum tensor is given by
\begin{equation} 
T_{\mu\nu}
= {1\over 8\pi G_N} ( 
R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R + \Lambda g_{\mu\nu} ).
\end{equation}
with $G_N$ the Newton constant.
According to the equation of motion \eqref{EOMV} this energy momentum tensor seems to vanish on shell, but it does not. We emphasize that it vanishes except a singularity located at $r=0$. 
This singularity contributes to the charge.

In order to evaluate the contribution, 
we first compute the energy by expanding the energy momentum tensor in terms of the weak field around $r=\infty$. That is,
we write 
$ f=\bar f +\delta f$ with $\delta f = -2G_N M/r^{d-3}$
and expand the stress tensor perturbatively around $\bar f$ 
to extract a pole.
This can be done by separating the metric into the regular part $\bar g_{\mu\nu}$ and the singular part $h_{\mu\nu}$, the latter of which is given by
\begin{eqnarray}
h_{\mu\nu} dx^\mu dx^\nu &=& -\delta f (dx^0)^2 +\left(\frac{1}{f} -\frac{1}{\bar f}\right) dr^2.
\end{eqnarray}  
At the leading order, we  have
\begin{eqnarray}
T^0{}_0
&=&-{1\over2} \left( {1\over \sqrt{|\bar g|}} \partial_\mu( \sqrt{|\bar g|} \bar g^{\mu\nu} \partial_\nu h_{0}^0) - \bar f^{-1}\bar f'^2 h^0_{0} \right) 
+  \bar\nabla^{0} \bar\nabla_\sigma h^\sigma_{0}- h_0^0 -{1\over2} \bar\nabla_{\mu} \bar\nabla_\sigma h^{\sigma\mu} + \cdots ,~~~~~
\end{eqnarray}
where $\bar\nabla_\mu$ is the covariant derivative with respect to the metric $\bar g_{\mu\nu}$,
and the ellipsis represents the higher order terms with respect to $h$. 
This must vanish except at the origin, and indeed this can be written as  
\begin{eqnarray}
T^0{}_0 &=&  \frac{d-2}{16\pi G_N r^{d-2}}\partial_r \left( r^{d-3}\delta f\right) + \cdots , 
\label{eq:T00}
\end{eqnarray}
which has the desired property. 
Plugging this into \eqref{mass} we can compute the charge as 
\begin{eqnarray} 
E
&=&- \int d^{d-1}\vec x\,  \sqrt{\vert \tilde g\vert} \frac{d-2}{16\pi G_N}\partial_r \left( r^{d-3} \delta f \right) 
= \rho V_{d-2},~~~
\label{massS}
\end{eqnarray}
where $V_{d-2} = \int d^{d-2} x \sqrt{\vert \tilde g\vert }$ is the volume of the Einstein manifold
with $\tilde g$ being the determinant of $\tilde g_{ij}$, and
$\rho = (d-2) M/(8\pi)$ is the energy or mass density. To evaluate $r$ integral we employ the Stokes' theorem.
Note that the higher order terms do not contribute to the surface integral. 
Our result reproduces the known result obtained by other methods. (For example, see eq.(2.5) in Ref.~\cite{Horowitz:1998ha}.)

On the other hand, it is also possible to compute the contribution of the singularity by the direct calculation of the Ricci tensor. 
From the direct calculation one finds that the
first term in Eq.~\eqref{eq:T00} is indeed exact. 
More detailed results are as follows.
\begin{eqnarray}
R^0{}_0 &=& -\frac{1}{2r^{d-2}}\partial_r\left( r^{d-2}\partial_r f(r) \right) = R^r{}_r,\\
R^i{}_j &=& \delta^i_j\left[\frac{(d-3)k}{r^2} -\frac{1}{r^{d-2}}\partial_r\left(r^{d-3} f(r)\right)\right], 
\end{eqnarray}
which leads to
a form of $T^\mu{}_\nu$ proportional to a delta function:
\begin{eqnarray}
 T^0{}_0 &=&  \frac{d-2}{16\pi G_N r^{d-2}}\partial_r  \left( r^{d-3} 
 \delta f(r)\theta (r)
 \right)
 = - \rho \frac{\delta(r)}{r^{d-2}} =T^r{}_r, \notag\\ 
 T^i{}_j&=& -\delta^i_j \frac{\rho}{d-2} \frac{r\partial_r\delta(r)}{r^{d-2}}
 = \delta^i_j \frac{\rho}{d-2} \frac{\delta(r)}{r^{d-2}}.
  \label{eq:T00B}
\end{eqnarray}
We here inserted the step function $\theta$ with $\theta(0)=0$ for the singular term in \eqref{f} 
to explicitly extract the singular contribution, 
and used $\delta(r)=\frac{d \theta(r)}{d r} $.
Writing $\delta(r)$ in terms of $\delta^{(d)}(\vec x)$, this agrees with a result at $d=4$ by distributional techniques\cite{Balasin:1993fn}.
Thus the matter energy momentum tensor $T^\mu{}_\nu$ for the Schwarzschild black hole can be understood as a distribution.

Although we do not encounter any mathematical problem to derive this result, one may wonder the physical validity to perform the volume integral of the constant $x^0$ slice, which becomes time-like inside the horizon.\footnote{
We would like to thank M.~Sasaki, T.~Shiromizu and S.~Sugimoto for raising this question and for discussion using the Penrose diagram.
}
To clear up this subtlety, let us consider a simplified situation where the cosmological constant vanishes with $k=1$. Then we can move to the Eddington-Finkelstein coordinates \cite{Eddington:1924aa,Finkelstein:1958zz}
\begin{eqnarray}
ds^2 = -(1+u) d\tau^2 -2u d\tau dr +(1-u) dr^2 +r^2\tilde g_{ij} dx^idx^j,
\end{eqnarray}
where $u =-2G_N M/r^{d-3}$ and we changed the time variable from $x^0$ to $\tau = x^0 -g(r) $ with  $
\frac{d g(r)}{d r} =\frac{u}{1+u}$. Then the unit normal vector to the constant $\tau$ slice defined by $n_\mu= -{1\over \sqrt{ 1 - u}} \delta^\tau_\mu$ is always well-defined and time-like for any nonzero $r$. 
The conserved energy, which takes the same form due to its manifest general covariance, is computed as
\begin{eqnarray}
E &=& -\int_{M_\tau} d^{d-1}\vec x\, \sqrt{\vert g\vert} T^\tau{}_\tau = \rho V_{d-2},
\end{eqnarray}
which agrees with Eq.~\eqref{massS},
where $T^\tau{}_\tau = \frac{(d-2)}{16\pi G_N} \frac{\partial_r( r^{d-3} u)}{r^{d-2}}.$



\subsection{Reissner-Nordstr\"om black hole} 
Below we perform the same computation of a mass of a charged black hole in general $d$ dimensions,
for which the the metric is given in the same form as \eqref{eq:SBH} except that  $f(r)$ is replaced by  $f_q(r) = f (r)+{d-3 \over d-2} 8\pi G_N q^2 r^{-2(d-3)}$, 
together with the gauge potential 
$A_\mu = (-\frac{q}{r^{d-3}} +\frac{q}{r_+^{d-3}}) \delta^0_{\mu}
$, where $q, r_+$ are constants\cite{Chamblin:1999tk}.
This configuration of gravitational and gauge fields satisfies the equations of motion given by
\begin{eqnarray}
G_{\mu\nu} +\Lambda g_{\mu\nu} &=&  
8\pi G_N ( T^G_{\mu\nu} +T^A_{\mu\nu}),
~ \nabla_\mu F^\mu{}_{\nu} =J_\nu, ~~~
\end{eqnarray}
where
$F_{\mu\nu} :=\nabla_\mu A_\nu -\nabla_\nu A_\mu $ and
$T^A_{\mu\nu} := F_{\mu}{}^\alpha F_{\nu\alpha} -\frac{1}{4} g_{\mu\nu} F_{\alpha\beta} F^{\alpha\beta}$.
Here $ T^G_{\mu\nu}$ and $J_\nu $ explicitly represent the singular contributions of the metric and the gauge potential at $r=0$, respectively.
Explicitly $(T^G + T^A)^0{}_0$ is given in \eqref{eq:T00B} by replacing $\delta f$ with
$\delta f_q= \delta f + {d-3 \over d-2} 8\pi G_N q^2 r^{-2(d-3)}$.
This also agrees with the distributional result at $d=4$\cite{Balasin:1993kf}.

Since this metric is also static, the energy defined by \eqref{mass} is conserved. 
However this charge diverges, due to the contribution of the electromagnetic field.
Physically, this divergence can be interpreted as 
that of
the self-energy for the charged point particle.
Indeed it remains even for the flat space-time with $M=0$ and $\Lambda=0$.
Therefore,
classically,  the charged black hole has the infinite energy due to the infinite electromagnetic
energy. Thus the renormalization as well as the quantization of the gauge field on the curved space are needed to fix this problem, as is the case
 on the flat space.

Fortunately, since $\nabla_\mu(T^G)^{\mu}{}_ 0 =0$ (thus  $\nabla_\mu(T^A)^{\mu}{}_ 0 =0$),  we can define an energy from the covariantly conserved $T^G$  alone without electromagnetic energy as
\begin{eqnarray}
(T^G)^0{}_0 &=& -\frac{(d-2)}{16\pi G_d r^{d-2} }\partial_r\left(r^{d-3} \delta f(r)\right),
\end{eqnarray}
where  $\delta f$ is given before.
We thus obtain
\begin{eqnarray}
E_G &=& \int d^{d-2} \vec x\,\int dr \sqrt{\vert g\vert} (T^{G})^0{}_0 \xi^0 = V_{d-2}\rho ,
\end{eqnarray}
which reproduces the result in 
\cite{Chamblin:1999tk} 
as a special case choosing a sphere as the internal manifold.

This system allows another conserved quantity, thanks to the invariance  under the $U(1)$ gauge transformation by $\delta A_\mu =\partial_\mu \theta$, which leads to
\begin{eqnarray}
\partial_\mu j^\mu &=& 0, \quad j^\mu = \nabla_\nu \left(\sqrt{\vert g\vert} F^{\mu\nu}\right)
\end{eqnarray}
without using the Maxwell equation.  
According to our prescription, $Q_c= \int d^{d-2} x\,\int dr\, \sqrt{\vert g\vert} J^0$
with $J^0=j^0/\sqrt{\vert g\vert}$
gives
the conserved electric charge, which is evaluated as $Q_c = V_{d-2} (d-3) q .$
At $d=4$ for $k>0$, $Q_c= 4\pi q$.

\subsection{BTZ black hole} 

As a final example, we compute a charge different from a mass. 
To this end we consider a BTZ black hole and compute its angular momentum \cite{Banados:1992wn}. 
The metric 
\begin{equation} 
ds^2 = -f(r) dt^2 + {1\over f(r)} dr^2 + r^2(d\phi -\omega(r)dt) ^2, 
\end{equation}
where 
\begin{equation}
f(r)= {r^2 \over L^2} - 2G_N M \theta(r) + {G_N^2 J^2 \over 4r^2} , ~~ \omega(r) = {G_N J \over 2 r^2} , 
\end{equation}
with $M,J$ are constants, 
satisfies the Einstein equation in three dimensions. We insert the step function to the constant part to emphasize that this solution is valid except the origin. 

This BTZ black hole has not only a Killing vector with respect to the time translation but also the one which rotates the system, $\xi^\mu = \delta^\mu_\phi$. 
As in the previous cases the first one defines the mass, which can be similarly computed as $E=\frac{M}{4}$. 
On the other hand, the second Killing vector define an angular momentum: 
\beal{
P_\phi=& \int d^{2}x \sqrt{|g|} T^0{}_{\phi}. 
}
$T^0{}_{\phi}$ is computed from the Einstein tensor as 
$T^0{}_{\phi} 
=-\frac{1}{ 16\pi G_N r} \partial_r \left(r^3 \omega'(r) \right) $. Thus we find $P_\phi ={J\over 8} $, which reproduces the known result \cite{Banados:1992wn}. 


\section{Mass of a compact star}

In this section we apply a conserved charge to the computation of the total energy of a static compact star setting $k=1$ with $d\geq3$. 
\subsection{Oppenheimer-Volkoff equation}
Let us consider 
a spherically symmetric system such that the metric is given by
\beal{
d s^2 =& -f(r) (dx^0)^2  + {h(r) } dr^2  + r^2\tilde g_{ij} dx^i dx^j,
\label{star}
}
and the matter energy momentum tensor is described by the perfect fluid as 
\begin{eqnarray}
T^0{}_0 &=& -  \rho(r), \quad 
T^r{}_r =  P(r), \quad
T^i_j = \delta^i_j P(r),
\label{perfectFulid}
\end{eqnarray}
where $\rho(r)$ is the 
density and $P(r)$ is the pressure. 
From the Einstein equation, we can derive useful formulas 
\beal{
P(r)+\rho(r) =& \frac{(d-2) \left(\frac{h'(r)}{h(r)} +\frac{ f'(r)}{f(r)} \right)}{16\pi G_N r h(r)}, ~~~ 
P'(r)   = -\frac{(P(r)+\rho(r))}{2} \frac{f'(r)}{f(r)}. 
\label{useful}
}
From these equations we can derive
the Oppenheimer-Volkoff 
or TOV
equation 
\cite{Oppenheimer:1939ne,Tolman:1939jz}
\begin{eqnarray}
\frac{d P(r)}{d r} &=&- \frac{ G_N \rho(r) M(r)}{ r^{d-2}} 
\Biggl( 1 +\frac{P(r)}{\rho(r)}\Biggr) h(r)
 \left\{ d-3+ \frac{ r^{d-1}}{(d-2)M(r)}
 \Biggl(8\pi P(r)  -\frac{2\Lambda}{(d-1) G_N} \Biggr)\right\},\notag\\
\end{eqnarray}
where 
\begin{eqnarray}
\frac{1}{h(r)} &:=&\frac{-2\Lambda r^2}{(d-2)(d-1)} + 
1 -\frac{2G_N M(r)}{ r^{d-3}} ,
\end{eqnarray}
and
\begin{eqnarray}
M(r) &=& \frac{8\pi}{d-2} \int_0^{r} ds\, s^{d-2} \rho(s).
\label{MS}
\end{eqnarray}
In order for this system to describe a static compact star we impose a boundary condition such that the pressure vanishes and the energy momentum tensor  is covariantly conserved at the surface of the star.
This makes the pressure and the density vanish outside the star, so that the metric becomes the Schwarzschild outside the star
, namely 
$f(r) =1/h(r)$ with $M=M(R)$ in \eqref{eq:SBH}, where $R$ is the radius of the star.
We here remark that this model does not admit the zero radius limit with positive $M(R)$ fixed, since the stress tensor of the Schwarzschild black hole given by eq.~\eqref{eq:T00B} is not the form of the perfect fluid, \eqref{perfectFulid}.
This shows that there needs to be a certain dynamical process for a star described by \eqref{star} and \eqref{perfectFulid} to collapse into the Schwarzschild blackhole. We shall confirm this explicitly in the analysis for constant density.

\subsection{Total energy with an equation of state
} 
We define the total energy of this system 
by the conserved charge corresponding to the Killing vector $\xi^\mu = -\delta^\mu_0$.
\footnote{
There is a traditional argument for the definition of the total energy of a compact star and its
interpretation, which is different from our result presented below. 
Just for clarity and to avoid confusion, we present our argument in main text and make the comparison to a traditional result in appendix.}
Then it can be computed as 
\begin{eqnarray}
E 
&=& V_{d-2} \int_0^R dr \sqrt{f(r)h(r)}r^{d-2}  \rho(r).
\label{TotalEnergy}
\end{eqnarray}
Employing \eqref{useful}, \eqref{MS} ,
we find 
\beal{
E 
=&  \frac{(d-2) V_{d-2}}{8\pi} ( M(R) +\Delta M ), \label{eq:NS_mass}
\\
\Delta M=& -\frac{8\pi G_N}{d-2} \int_0^R dr\, \sqrt{f(r) h(r)^3} r M(r)\{\rho(r) + P(r)\} .
\label{binding}
}
The first leading term gives the gravitational mass of the star also called as the Misner-Sharp mass\cite{Misner:1964je}, which is directly related to an actual observable quantity appearing in the metric around the spacial infinity, 
while the second term is a deviation from the gravitational mass, which corresponds to the gravitational binding energy as we shall see. 

In order to investigate the deviation term in more detail,\footnote{
We would like to thank M.~Shibata, S.~Yamaguchi and S.~Mukohyama to motivate us to study the deviation from the Misner-Sharp mass.
} let us impose the matter consisting of the star to satisfy an equation of state 
\begin{equation} 
P(r) = w\rho(r).
\label{EoS}
\end{equation}
Plugging this into \eqref{useful} leads to 
$f= \rho^{-\frac{2w}{(1+w)} },$ where we fixed the integration constant to satisfy $f\to1$ in the Newtonian limit $w\to0$.\footnote{
This integration constant can be also absorbed by rescaling the time variable.
} 
Then the deviation term becomes  
\begin{equation}
\Delta M  = - \frac{8\pi G_N (1+w)}{d-2} \int_0^R dr r  \rho^{\frac{1}{(1+w)} }{1\over \left({-2\Lambda r^2 \over (d-1)(d-2)} + 1 - {2 G_N M(r) \over r^{d-3}}\right)^{3/2} }   M(r) .
\end{equation}
Let us expand this in terms of the cosmological constant and the Newton constant. 
The leading and next-to leading terms are given by 
\beal{
\Delta M_1  =& - \frac{8\pi G_N (1+w)}{d-2} \int_0^R dr r  \rho^{\frac{1}{(1+w)} }   M(r),   \\
\Delta M_2  =& - \frac{24\pi G_N (1+w)}{d-2} \int_0^R dr  \rho^{\frac{1}{(1+w)} } \left({\Lambda r^3 \over (d-1)(d-2)} + { G_N M(r) \over r^{d-4}}\right)  M(r). 
}
Here $\Delta M_1$ represents the leading order of the gravitational self-interaction energy of the matter inside the star obeying the equation of state \eqref{EoS}.
In particular, at four dimensions in the Newtonian limit $w\to 0$, this can be rewritten as 
\begin{eqnarray}
\Delta M_1 =  -\frac{G_N}{2}  \int 
d^3x\, d^3y  \, \frac{\rho(\vec x) \rho(\vec y)}{\vert \vec x-\vec y\vert} ,
\end{eqnarray}
which is nothing but Newtonian gravitational energy inside the star including the symmetric factor $1/2$.

\subsection{Estimation with constant density}
\label{constantDensity} 
In order to estimate the size of the correction $\Delta M$,
let us consider a case of constant density $\rho(r) = \rho_0$ and vanishing cosmological constant. 
In this case, we can compute $M(r)$ given in \eqref{MS} as $
M(r) = \frac{8\pi r^{d-1}\rho_0}{(d-2)(d-1)}$, which leads to  $h(r)^{-1} = 1 - {r^2 \over r_0^2}$ with $r_0 =\sqrt \frac{(d-2) (d-1)}{16\pi G_N \rho_0 } $. 
The TOV equation with the boundary condition is easily solved as 
\beal{
P(r)=&\frac{\rho_{0} (d-3) \left(\sqrt{r_0^2 - r^2}-\sqrt{r_0^2 - R^2}\right) }{-(d-3) \sqrt{r_0^2 - r^2} +(d-1)   \sqrt{r_0^2 - R^2}}. 
}
For a stable star the pressure has to be finite for all $r$. 
 This leads to an inequality\footnote{
This inequality leads to the lower bound for the radius of the star $R_{\rm min}$ as $
R^{d-3} \ge  \frac{(d-1)^2}{2(d-2)} G_N M_0 =R_{\rm min}^{d-3}$ for 
keeping $M(R)=M_0$ independent on $R$. 
This is consistent with the previous argument that this system does not admit the zero size limit fixing a positive $M(R)=M_0$.
} 
\begin{equation} 
R < \frac{2 \sqrt{d-2} }{d-1 }r_0 =:R_{*}.
\end{equation}
$\Delta M$ can be written as 
\beal{
\Delta M
=-{(d-1) \sqrt{r_0^2 -R^2} \over 4G_N r_0^2 } \int_0^R dr { r^d  \over (r_0^2 -r^2)^{3/2}}.
}
At four dimensions, this is computed as 
\begin{equation} 
F(R):= \frac{\Delta M}{M(R)}
=-\frac{3}{8G_NM(R)} \left(-\frac{R^3}{r_0^2}+3 R-3 \sqrt{r_0^2-R^2} \arctan\left(\frac{R}{\sqrt{r_0^2-R^2}}\right)\right),
\end{equation}
where $F(R)$ turns out to be a monotonically decreasing function 
from $F(0)=0$ to $F(R_{*})\simeq - 0.68$. 
Therefore the total energy $E$ can be about one third of $M(R)$ so that the gravitational binding energy could be considerably large. 
In particular the total energy of a static compact star becomes smaller than the gravitational mass observed at spacial infinity due to the gravitational self-interaction. 



\section{Discussion}
\label{sec:Discussion}
We have proposed a general definition of a conserved charge from any covariantly conserved current,
which requires no specific asymptotic behaviors/approximations for the metric, or no subtraction of boundary contributions, 
as long as a Killing vector exists.
Our definition has reproduced the mass, electric charge and angular momentum of black holes known in the literature.  
Since the presented
formula requires the
matter energy momentum tensor to define the mass,
it is clear that black holes inevitably
have non-zero matter
energy momentum tensor at singularity.
(See also \cite{Buoninfante:2018rlq} in the case of higher derivative gravity.)
One of important consequences in this paper is that any black hole is not a vacuum solution in this treatment.
\footnote{ Traditionally black holes including the ones dealt with in this paper may have been regarded as vacuum solutions in spacetime with their singularity deficient. The method presented in this paper does not return inconsistent results with the traditional treatment, but is applicable to any gravitational object on a general curved spacetime such as a black hole and a compact star.}
This is similar to an electron in electrodynamics: its electric charge distribution is described by a delta function at its position and it is not regarded as a vacuum solution.
We have also demonstrated that the total energy 
of any spherically symmetric compact star defined as the conserved charge cannot be written as a surface term alone, and it contains the gravitational binding energy in addition to the gravitational mass observed at the asymptotic spatial infinity. 

Our definition of charges is formulated as a  generalization of ones used in the flat spacetime to be available in a general curved spacetime. This is achieved to enjoy general covariance manifestly. Therefore our definition has clear physical meaning and is of a generic use with precision compared to other definitions. We discuss difference of a couple of other approaches more in appendix. 

In this paper we focused on a few well-known blackholes such as the Schwarzschild and BTZ ones to compute their charges. In these simple cases the matter stress tensor is
described by a delta function singularity at the origin. It would be interesting to study a more intricate blackhole whose singularity is not point-like any more. In such a case we expect that matter energy momentum tensor is still described by a certain integrable distribution and our definition of charges is valid and useful also as numerical evaluation. 
(Distributional expressions for the Kerr black hole can be found in \cite{Balasin:1993kf}.)

Our proposal is quite generic, so we expect plenty of applications in future. 
As such a potential application, we consider a more  general case where any Killing vectors do not exist. 
We can still consider a charge or a generator associated with a general vector field $\xi^\mu$ as
\begin{equation}
Q[\xi](t) = \int_{M_t} d^{d-1} \vec x\,  \sqrt{\vert g\vert} T^0{}_\nu \xi^\nu .
\label{chargeglobal}
\end{equation}
Using the similar argument before, we obtain
\begin{eqnarray}
\frac{d Q[\xi]}{d t} &=& \int_{M_t} d^{d-1}\vec x  \sqrt{\vert g\vert}
\rho(x),\label{eq:conv} \quad
\rho (x) := \frac{\left( G^{\mu\nu} +\Lambda g^{\mu\nu}\right) }{16\pi G_N}  (\nabla_{\mu} \xi_{\nu}+\nabla_{\nu} \xi_{\mu})
\end{eqnarray}
where  $\rho=0 $ if  $\xi$ is a Killing vector. 
A change of the charge $Q[\xi]$ can be calculated by the volume integral of $\rho$, 
expressed in terms of the gravitational field  through the Einstein equation {with the manifest covariance being kept. 
Eq. \eqref{eq:conv} may give a hint for a {\it  generic conservation equation} 
in general relativity. Note that this does not require any pseudo-tensor. This argument will hold not only for a Lorentzian manifold but also for a more general one.
We will return to this interesting problem in future studies. 

\section*{Acknowledgement}
We would like to thank S.~Mukoyama,  J.~H.~Park,  M.~Sasaki,  M.~Shibata, T.~Shiromizu, S.~Sugimoto and S.~Yamaguchi for useful discussions and 
valuable
comments.
This work is supported in part by the Grant-in-Aid of the Japanese Ministry of Education, Sciences and Technology, Sports and Culture (MEXT) for Scientific Research (Nos.~JP16H03978,  JP18K03620, JP18H05236, JP19K03847). 
S. A. is also supported in part  
by a priority issue (Elucidation of the fundamental laws and evolution of the universe) to be tackled by using Post ``K" Computer, and by Joint Institute for Computational Fundamental Science (JICFuS).
T. O. would like to thank YITP for their kind hospitality during his
stay for the sabbatical leave from his home institute.

