\section{Annexe}

The laser source is a titanium-sapphire mode-locked oscillator delivering $\sim 140 \textrm{fs}$ pulses ($\Delta \lambda \simeq 6 \textrm{nm}$ FWHM) centered at $795 \textrm{nm}$ with a repetition rate of $76 \textrm{MHz}$. Its second harmonic at $397 \textrm{nm}$ serves to synchronously pump an OPO, which consists of a 2mm BIBO crystal contained within a 4m long ring cavity exhibiting a finesse of 27 and escape efficiency of 95\%. In this configuration, the OPO reaches its threshold at a mean pump power of $\sim 75 \textrm{mW}$. Below threshold and in the absence of a seed, the device generates squeezed vacuum with a squeezing level of $\sim -6 \textrm{dB}$ (corrected) when projected onto the unshaped LO as seen in Fig.~\ref{fig-fullsqz}.

\PDFfig{fig-fullsqz}{fullsqz-v2.pdf}{Noise trace observed with an unshaped LO pulse shape. The state exhibits squeezing and anti-squeezing values of $-5.9 \textrm{dB}$ and $7.8 \textrm{dB}$, respectively.}{80mm}

In order to lock the cavity length to the inter-pulse spacing, a near-infrared beam is phase-modulated at $1.7 \textrm{MHz}$ with an electro-optic modulator (EOM) and injected into the cavity in a direction counter-propagating to the pump and seed. Locking of the cavity is then accomplished with a Pound-Drever-Hall strategy.

Frequency correlations of the state are investigated with homodyne detection where the local oscillator pulse form is manipulated with a pulse shaper. The $4 f$-configuration shaper is constructed in a reflective geometry with a programmable 512\,x\,512-element liquid-crystal modulator in the Fourier plane.  Application of a periodic spatial grating to the spatial light modulator induces diffraction of the spectrally-dispersed light. The amplitude and phase of the diffracted spectrum are independently controlled by the groove depth and position of the spatial grating, respectively \cite{nelson2005}.

The OPO output is projected by the homodyne detection technique onto the spectrally filtered LO pulse shape to assess the quadrature noise content of the corresponding spectral region. Light detection is performed with silicon photodiodes ($\sim 90\%$ detection efficiency, 100MHz detection bandwidth), and the noise level of the squeezed vacuum is examined at 1MHz. The homodyne visibility is $92 \%$. The cumulative loss of the system is taken to be $\sim 25 \%$ and the measured signals are corrected accordingly.

\PDFfig{fig-pquad}{figure-pquad.pdf}{Noise correlation matrix defined as $C^{p}_{i,j} = \langle p_{i}p_{j} \rangle / \sqrt{\langle p_{i}^2 \rangle \langle p_{j}^2 \rangle} -  \delta_{i,j} \langle p_{i}^2 \rangle_{\textrm{vac}} / \langle p_{i}^2 \rangle$ for the $\hat{p}$ quadrature. The diagonal displays excess noise while the off-diagonal elements demonstrate quantum correlations among disparate spectral regions.}{75mm}

The optical spectrum of the local oscillator is divided into ten frequency bands of equal energy as seen in Fig.~1 of the main text, and spectral holes between the individual regions are imposed to ensure the absence of any spectral overlap. Importantly, the supplemental loss incurred from including these gaps is not accounted for when correcting the noise levels.
The $\hat{x}$ quadrature is defined as the field quadrature of lowest noise for the unshaped LO pulse.
The noise dependence on LO phase is then examined for all 55 possible pulse shapes and compared to that of the unshaped LO reference. Each individual trace is observed to follow the phase dependence of the reference, which confirms that the lowest noise level for every spectral combination is present in the $\hat{x}$ quadrature, i.e., there is no rotation of the squeezing ellipse between successive measurements.

The noise content of both the $\hat{x}$- and $\hat{p}$-quadratures for each spectral region and all possible pairs of regions are measured as shown in Fig.~\ref{fig-sqzTraces}. The individual covariance elements are constructed according to the following equation:

\begin{eqnarray}
\langle x_{i} x_{j} \rangle &=&
\left[ \langle (x_{i} + x_{j})^2 \rangle - \frac{P_{i}}{P_{i}+P_{j}} \langle x_{i}^2 \rangle - \frac{P_{j}}{P_{i}+P_{j}} \langle x_{j}^2 \rangle \right] \nonumber \\ && \times \, \frac{P_{i}+P_{j}}{2 \sqrt{P_{i} P_{j}}},
\end{eqnarray}
where $P_{i}$ and $P_{j}$ are the optical powers of frequency bands $i$ and $j$, respectively, which are measured with the homodyne photodiodes.

\textbf{Data Analysis}. Noise fluctuations and correlations departing from the vacuum level for the $\hat{p}$-quadrature are represented in Fig.~\ref{fig-pquad}. Similar to the $\hat{x}$ quadrature, the vast majority of both the excess noise (diagonal peaks) and correlation (off-diagonal peaks) are localized in the spectral wings. As opposed to the $\hat{x}$ quadrature, however, the correlations within the $\hat{p}$ quadrature are positive.

The purity $\mathcal{P}$, which is an intrinsic property of the state, is accessible from the covariance matrix with the relation $\mathcal{P} = 1 / \sqrt{\textrm{det}(\mathcal{C})}$, where $\mathcal{C}$ is the full covariance matrix of the state \cite{Adesso2007}. This relationship reveals that $\mathcal{P} \simeq 0.45 \pm 0.09$.
The full covariance matrix is then eigendecomposed, where it is observed that although the individual $\hat{x}$ and $\hat{p}$ block eigenvectors are very similar, they are not identical.
In general, two matrices are simultaneously diagonalizeable only if they commute. For instance, the commutator of the $\hat{x}$ and $\hat{p}$ covariance matrices for two-mode squeezed vacuum is zero even in the presence of spectrally-uniform loss, and a single eigenvector set perfectly decorrelates both quadratures. The introduction of a spectrally-dependent loss, however, spoils the commutation of the $\hat{x}$ and $\hat{p}$ blocks, and the magnitude of a commutator element is directly proportional to the spectral asymmetry of the loss. The output coupler of the SPOPO does not have a perfectly flat transmission window over the bandwidth of the pulse, and the magnitude of the commutator elements is largest in the spectral wings. Consequently, this asymmetry explains the fact that it is not possible to perfectly remove both the $\hat{x}$ and $\hat{p}$ correlation with a single basis change.

\begin{table}
\begin{center}
\caption{Squeezing values and uncertainties (linear scale) for the modes retrieved from the twenty dimensional covariance matrix. Eight modes are clearly non-classical while a ninth mode is not included in a conservative estimate. \label{tab-sqz}}
\medskip
{
\begin{tabular}{c|c|c}
\hline
{Mode} & Squeezing Level & Anti-Squeezing Level \\
\hline
1 & $0.38 \pm 0.07$ & $3.86 \pm 0.12$ \\
2 & $0.48 \pm 0.06$ & $3.62 \pm 0.07$ \\
3 & $0.58 \pm 0.07$ & $2.74 \pm 0.10$ \\
4 & $0.74 \pm 0.06$ & $1.96 \pm 0.06$ \\
5 & $0.83 \pm 0.05$ & $1.41 \pm 0.06$ \\
6 & $0.90 \pm 0.03$ & $1.17 \pm 0.04$ \\
7 & $0.93 \pm 0.03$ & $1.11 \pm 0.03$ \\
8 & $0.95 \pm 0.02$ & $1.06 \pm 0.03$ \\
9 & $0.97 \pm 0.02$ & $1.03 \pm 0.03$ \\
10 & $0.98 \pm 0.02$ & $1.00 \pm 0.02$ \\
\end{tabular}
}
\end{center}
\end{table}

Yet, the commutator elements have a small magnitude.
Therefore, as a means of generating a mode set that efficiently removes the correlation of the original frequency bands, the ten most noise-robust eigenmodes (see statistical sampling procedure detailed below) are selected from the total eigenvector set. These vectors, which are drawn from both the $\hat{x}$ and $\hat{p}$ blocks, are orthogonalized with a Gram-Schmidt procedure, and the covariance matrix is subsequently re-expressed in this basis. The resulting matrix is nearly diagonal and contains the squeezing values for each orthogonalized mode on its diagonal.

The error level associated with each of these squeezing values is subsequently estimated with a stochastic method.
The data underlying every trace of Fig.~\ref{fig-sqzTraces} is collected for twenty seconds in order to accrue statistics about the individual squeezing and anti-squeezing values. The noise value for a given trace is then drawn from a normal distribution with a mean provided by the red (green) line, which represents the average of all identified peaks (valleys), and a variance specified by the variance of the peak (valley) amplitudes. After randomly sampling from these 110 normal distributions, the corresponding covariance matrix is constructed in the appropriate manner. A collection of $10^{4}$ individual matrices is assembled, and the Gram-Schmidt orthogonalization procedure is repeated for every matrix, which provides both a squeezing spectrum and orthonormal mode set. The mean squeezing spectrum from this procedure is shown in Fig.~\ref{fig-eigentrace} and the individual noise levels are detailed in Table~\ref{tab-sqz}. Likewise, the average mode structures corresponding to the six highest squeezing values are shown in Fig.~4 of the main text.

\PDFfig{fig-eigentrace}{eigentrace-v2.pdf}{Mean noise levels and uncertainties (dB) for each of the orthogonalized Gramm-Schmidt modes. The eight modes that are clearly distinguished from the vacuum level are indicated in blue, and modes ascribed to vacuum fluctuations are indicated in red.}{85mm}

\PDFfigDouble{fig-sqzTraces}{sqz-Traces-v2.pdf}{Squeezing traces for determination of the ten frequency band covariance matrix. The traces are sequentially arranged according to the index scheme $\{i,j\}$, which indicates that the pulse shaper selects frequency bands $i$ and $j$ (e.g., \{1,1\}, \{1,2\},...,\{1,10\},\{2,2\},\{2,3\},...). In each trace, the black line represents the vacuum level while the green and red lines indicate the average squeezing and anti-squeezing levels, respectively. The LO phase is depicted with arbitrary units.}{180mm}

The mean spectrum is generally noise-robust and reveals that eight modes are statistically distinct from the vacuum level. Although a ninth mode may also be considered to exhibit non-classicality, its uncertainty is at the vacuum level, and is therefore not included in a conservative estimate.

\

The research is supported by the ANR project \emph{Qualitime} and
the ERC starting grant \emph{Frecquam}. C. Fabre is a member of
the Institut Universitaire de France.

\bibliographystyle{apsrev}


\section{Introduction}

People have been concerned with how to define the energy-momentum and angular momentum of gravitating systems since even before general relativity was born. Different pseudotensors were proposed to describe the local gravitational energy-momentum density. This choice between the different available expressions is one  ambiguity in characterizing the energy-momentum of any local region of space-time. Moreover there is an additional ambiguity: since a pseudotensor is not a covariant object, it leads to values which inherently depend on the reference frame. This feature is actually an inevitable consequence of the fundamental geometric nature of gravity. Mathematically it is due to diffeomorphism invariance (as shown by Noether in her famous paper~\cite{noether}); physically it has been said that the equivalence principle forbids the localization of gravitational energy~\cite{MTW73}. The more modern view is that energy-momentum should be regarded as being \emph{quasi-local} (associated with a closed 2-surface) rather than local. The various quasi-local ideas have been comprehensively reviewed~\cite{Sza09}. There have been many quasi-local proposals, but no consensus.  In fact the quasi-local approach has ambiguities analogous to the two already mentioned.  However, it has been argued that the Hamiltonian approach, in particular the covariant Hamiltonian formalism~\cite{N91,N95,N99,N00,Nester04,N05,N08}, can tame both of these ambiguities, clarifying their geometric and physical nature.    The quasi-local values are then determined by the Hamiltonian boundary term.  A preferred boundary term has been identified; this gives a resolution to the first type of ambiguity. Recent work has finally led to identifying a good procedure for addressing the second type of ambiguity: how to best select the necessary reference values on the boundary~\cite{ae100}. We now test this program on axisymmetric spacetimes, showing that it works well; in particular it yields sensible energy and angular momentum values for the Kerr metric.

Early quasi-local investigations (an outstanding one is the seminal work of Brown and York~\cite{BY93}) naturally explicitly worked out quasi-local energy values for spherically symmetric spacetimes.  As far as we know Matrinez~\cite{Mart94} was the first to consider an axisymmetrical spacetime.  We are concerned with the covariant Hamiltonian formalism; within this approach the quasi-local energy for spherically symmetric spacetime was examined in depth by Wu {\it et al}.~\cite{Wu10,Wu11,Liu11,Wu12}  Here we will now go further: using the results of the covariant Hamiltonian formalism we consider an axisymmetric spacetime, use a new approach to find the best matched reference, and determine both the optimal quasi-local energy and angular momentum.

\section{Covariant Hamiltonian Formalism}

Detailed discussions of the covariant Hamiltonian formalism as developed by our research group can be found in the references~\cite{N91,N95,N99,N00,Nester04,N05,N08}. Here we give a brief summary of some key features.

A \emph{first order Lagrangian 4-form} for a $k$-form field $\phi$ has the form $\mathcal{L}=d\phi\wedge p-\Lambda(\phi,p)$. A vector field $N$ can be regarded as an ``infinitesimal translation''. It generates a local infinitesimal diffeomorphism under which the action associated with $\mathcal{L}$ should be invariant. Using the well-known Lie derivative formula $\pounds_N \equiv di_N + i_Nd$, the associated conserved\footnote{Conserved in the first Noether theorem sense: $d{\cal H}(N) \propto $ field equations.} Noether current 3-form---which is moreover the Hamiltonian---
\begin{equation}
\mathcal{H}(N) := \pounds_{N} \phi \wedge p - i_{N} \mathcal{L} \equiv N^{\mu} \mathcal{H}_{\mu} + d\mathcal{B}(N),
\end{equation}
can be split as indicated into two pieces with distinct roles. The part algebraic in the vector field generates the Hamiltonian equations; from Noether's second theorem concerning local invariance with respect to the choice of $N$ one finds that it is proportional to field equations and thus vanishes ``on shell''. Then the \emph{value} of the Hamiltonian,
\begin{eqnarray} \label{eq1}
E(N, \Sigma) = H(N, \Sigma) = \int_{\Sigma} \mathcal{H}(N) = \oint_{\partial\Sigma} \mathcal{B}(N),
\end{eqnarray}
is {\em quasi-local\/}; it is determined entirely by quantities on the closed 2-surface $S := \partial\Sigma$.

For \emph{any} choice of $N$ this expression defines a conserved quasi-local quantity.  To get a physically meaningful result one needs to have a good way of selecting the specific vectors $N$ that correspond to the quasi-local symmetries that yield the quasi-local energy-momentum and angular momentum/center-of-mass.

It must be noted that the Hamiltonian boundary 2-form ${\cal B}(N)$ can be modified ``by hand''; this would still preserve the conservation property, but change the conserved value. It was shown that in this way the Hamiltonian approach can include the various pseudotensors~\cite{Chang99}. But this ``Noether current'' ambiguity is tamed by the Hamiltonian \emph{boundary variation principle\/}: the boundary term in the variation of the Hamiltonian must be required to vanish.  This determines the related boundary conditions: different boundary terms are associated with different boundary conditions.

Here we consider specifically only Einstein's theory, general relativity (GR).
Our preferred\footnote{At spatial infinity it gives the accepted energy, momentum, angular momentum, and center-of-mass. At null infinity it gives the Bondi energy and energy flux, for small spheres it is a positive multiple of the Bel-Robinson tensor, it has a positivity property, and for spherically symmetric solutions it has the hoop property~\cite{Omur10}.}
boundary term for GR (which corresponds to holding the metric fixed on the boundary) is
\begin{eqnarray}\label{eq2}
16 \pi \mathcal{B}(N) = \Delta\Gamma^\alpha{}_\beta \wedge i_{N} \eta_{\alpha}{}^{\beta} + \overline{D}_{\beta} N^{\alpha} \Delta\eta_{\alpha}^{~\beta}, \label{B}
\end{eqnarray}
where $\eta^{\alpha\beta} := *(dx^\alpha\wedge dx^\beta)$ is the dual form basis, $\Gamma^\alpha{}_\beta = \Gamma^\alpha{}_{\beta\mu}dx^\mu$ is the connection one form, and $\Delta\Gamma := \Gamma - \bar\Gamma$, $\Delta\eta := \eta - \bar\eta$, with the bar refering to non-dynamical reference values
(which are needed only in a neighborhood of the boundary $S$).
The component form of our expression\footnote{An equivalent expression was proposed by Katz, Bi\v{c}\'ak, and Lynden-Bel~\cite{LBKB95,Katz:1996nr}.} is
\begin{eqnarray}
16 \pi \mathcal{B}(N) \label{eq3}
&=& \frac{1}{4} [ N^{\sigma} \sqrt{-g} g^{\gamma\beta} \Delta\Gamma^{\alpha}{}_{\beta\lambda} \delta^{\lambda\tau\rho}_{\alpha\gamma\sigma}
\nonumber\\
&& \quad + \overline{D}_{\beta} N^{\alpha} \Delta(\sqrt{-g} g^{\gamma\beta}) \delta^{\tau\rho}_{\alpha\gamma} ] \epsilon_{\tau\rho\mu\nu} dx^{\mu} \wedge dx^{\nu}.
\end{eqnarray}

The reference values are non-dynamic, they represent the ``vacuum'' or ground state; how to best select these values is our concern here.
The natural reference for dynamic geometry is Minkowski space.  Thus we need to select a specific Minkowski space in the neighborhood of the two boundary.  Any such choice can be specified by selecting some \emph{quasi-Minkowski coordinates}; in fact any four independent functions $y^i,\ i = 0 \dots 3$ defined in a neighborhood of the the 2-boundary $S$ will define a Minkowski metric by $d\bar s^2 = \bar g_{ij} dy^i dy^j, \ \bar g_{ij} = \hbox{diag}(-1,1,1,1)$.
A Killing field of this reference has the form $N^i = N^i_0 + \lambda_0^i{}_j y^j$, where the translation parameters $N_0^i$ and the boost-rotation parameters $\lambda_0^{ij} = \lambda_0^{[ij]}$ are constants.
Using this reference the 2-surface integral  of the Hamiltonian boundary term then gives a value of the form
\begin{equation}
\oint_S {\mathcal B}(N) = - N_0^i p_i(S) + \frac12 \lambda^0_{ij} J^{ij}(S)\,,
\end{equation}
which yields not only a quasi-local \emph{energy-momentum} but also a quasi-local \emph{angular momentum/center-of-mass}. These quasi-local values of course depend on the chosen reference.  As long as the reference approaches at an appropriate rate the flat Minkowski space at spatial infinity, the integrals $p_i(S),\ J^{ij}(S)$ in the spatial asymptotic limit will agree with accepted expressions for these quantities~\cite{MTW73,Regge:1974zd,Beig:1987zz,Szabados:2003yn,Sza06}.

\section{Optimal choice for reference}

To explicitly determine the specific values of the quasi-local quantities one needs some good way to choose the reference.  Minkowski spacetime is the natural choice, especially for asymptotically flat spacetimes~\cite{Kij97}.  However, as noted above, almost any four functions will determine some reference. With such freedom one can still get almost any value for the quasi-local quantities.  This freedom is the quasi-local version of the second type of ambiguity mentioned in the introduction in connection with the reference frame dependence of energy-momentum pseudotensor expressions.

Recently we have proposed a procedure~\cite{ae100} to fix the ``best'' choice of reference.  It has two features: 4D isometric matching\footnote{This was already proposed by Szabados back in 2000 (at a workshop in Hsinchu, Taiwan); he has since extensively explored this idea~\cite{Sza05}.} and energy optimization.  Here we shall show that our procedure works well for a certain class of axisymmetric spacetimes.

\subsection{4D isometric matching}

Previously the program of 4-dimensional isometric matching with energy extremization was explicitly worked out for the most important special case: spherically symmetric spaces (both static and dynamic~\cite{Wu11,Wu12}). Here we turn to what is probably the next most important special case: {\em axisymmetry}. More specifically we will consider Kerr-like metrics, with just one non-vanishing off diagonal term. (The most general axisymmetric metric and the general case without any symmetry are also being investigated; the results of these investigations will be presented in due course.)
Specifically, in terms of suitable spherical coordinates
$x^\mu = \{ x^{0}, x^{1}, x^{2}, x^{3} \} = \{ t, r, \theta, \varphi\}$, we consider metrics of the form
\begin{eqnarray}\label{eq20}
ds^2 = g_{00} dt^2 + g_{11} dr^2 + g_{22} d\theta^2 + g_{33} d\varphi^2 + 2 g_{03} dt d\varphi, \label{metric}
\end{eqnarray}
where the metric components are independent of $\varphi$.
Although it is not the most general metric of an axisymmetric spacetime, it is enough to unveil the property of a non-static spacetime with angular momentum so we can test our ideas.

We assume that the spherical-like coordinates can be chosen such that the boundary of 2-sphere $S = \partial\Sigma$ is given by $t=t_{0}, r=r_{0}$. In general the reference Minkowski metric is determined by four quasi-Minkowski coordinates defined in a neighborhood of $S$ by four suitable functions, $y^i = (T, X, Y, Z)$ of $(t, r, \theta, \varphi)$; from the $y^i$ we find $dy^i = y^i{}_\mu dx^\mu$.

Since the dynamic spacetime metric is not diagonal when expanded in terms of $x^\mu$ the reference metric must likewise allow for off diagonal terms.
In the axisymmetric case one convenient way to achieve this is to use a certain cylindrical form. Let
\begin{eqnarray}\label{eq21}
X = \rho \cos(\varphi + \Phi), \qquad Y = \rho \sin(\varphi + \Phi),
\end{eqnarray}
and assume that $(T, \rho, \Phi, Z)$ are axisymmetric, i.e., they do not depend on $\varphi$.
Then the reference Minkowski metric in a neighborhood of $S$ can be taken to have the form
\begin{eqnarray}\label{eq22}
d \bar{s}^2 &=& \bar g_{ij} dy^i dy^j
\nonumber\\
&=& - dT^2 + dX^2 + dY^2 + dZ^2 = - dT^2 + d\rho^2 + \rho^2 (d\varphi + d\Phi)^2 + dZ^2
\nonumber\\
&=& - ( T_{t} dt + T_{r} dr + T_{\theta} d\theta )^2 + ( \rho_{t} dt + \rho_{r} dr + \rho_{\theta} d\theta )^2
\nonumber\\
&& + \rho^2 ( d\varphi + \Phi_{t} dt + \Phi_{r} dr + \Phi_{\theta} d\theta )^2 + ( Z_{t} dt + Z_{r} dr + Z_{\theta} d\theta )^2
\\
&=& \bar g_{ij} y^i{}_\mu y^j{}_\nu dx^\mu dx^\nu = \bar g_{\mu\nu} dx^\mu dx^\nu. \nonumber
\end{eqnarray}
On the 2-sphere $S$ determined by $t = t_0, r = r_0$ the reference metric is a function of 12 functions of $\theta$, namely\footnote{In this work, we will just deal with a special type of metric with only one cross component. However, (\ref{eq21}), (\ref{eq22}),  and (\ref{isometric})--(\ref{g33}) are quite general, so we can straightforwardly extend our procedure to cover the most general axisymmetric metric.}
\begin{eqnarray}
\{ T_{t}, T_{r}, T_{\theta}, \rho_{t}, \rho_{r}, \rho, \Phi_{t}, \Phi_{r}, \Phi_{\theta}, Z_{t}, Z_{r}, Z_{\theta} \}
\end{eqnarray}
(note that on the 2-sphere $\rho$ determines $\rho_\theta := \partial_\theta\rho$).  We take $T_{r}$ and $T_{\theta}$ as a suitable set of controlling functions on the 2-sphere; then the other 10 can be found algebraically from the 4D isometric matching conditions, $ds^2 \doteq d\bar s^2$, i.e.,
\begin{equation}
g_{\mu\nu} \doteq \bar g_{\mu\nu} := \bar g_{ij} y^i{}_\mu y^j{}_\nu. \label{isometric}
\end{equation}
(Here the notation $\doteq$ means that the relation holds just on the 2-surface $S$.)
For a general metric these relations have the detailed form
\begin{eqnarray}
g_{00} &\doteq& -T_{t}^2 + \rho_{t}^2 + \rho^2 \Phi_{t}^2 + Z_{t}^2, \label{g00} \\
g_{01} &\doteq& -T_{t} T_{r} + \rho_{t} \rho_{r} + \rho^2 \Phi_{t} \Phi_{r} + Z_{t} Z_{r}, \label{g01} \\
g_{11} &\doteq& -T_{r}^2 + \rho_{r}^2 + \rho^2 \Phi_{r}^2 + Z_{r}^2, \label{g11} \\
g_{02} &\doteq& -T_{t} T_{\theta} + \rho_{t} \rho_{\theta} + \rho^2 \Phi_{t} \Phi_{\theta} + Z_{t} Z_{\theta}, \label{g02} \\
g_{03} &\doteq& \rho^2 \Phi_{t}, \label{g03} \\
g_{12} &\doteq& -T_{r} T_{\theta} + \rho_{r} \rho_{\theta} + \rho^2 \Phi_{r} \Phi_{\theta} + Z_{r} Z_{\theta}, \label{g12} \\
g_{13} &\doteq& \rho^2 \Phi_{r}, \label{g13} \\
g_{22} &\doteq& -T_{\theta}^{2} + \rho_{\theta}^{2} + \rho^{2} \Phi_{\theta}^{2} + Z_{\theta}^{2}, \label{g22} \\
g_{23} &\doteq& \rho^2 \Phi_{\theta}, \label{g23} \\
g_{33} &\doteq& \rho^2. \label{g33}
\end{eqnarray}

For our special axisymmetric metric with only one cross term it is easy to see from (\ref{g03}, \ref{g13}, \ref{g23}, \ref{g33}) that
\begin{equation}
\rho \doteq \sqrt{g_{33}}, \quad \Phi_{t} \doteq {g_{03}}/{g_{33}}, \quad \Phi_{r} \doteq g_{13}/g_{33} = 0, \quad \Phi_{\theta} \doteq g_{23}/g_{33} = 0.
\end{equation}
Recalling that $\rho_\theta$ follows from $\rho$, from (\ref{g22}) one can now get\footnote{Consider the spherical coordinate system in $\mathbf{R}^3$.
The relation between the Cartesian coordinates $X, Y, Z$ and the spherical coordinates $r, \theta, \varphi$ is $X = r \sin\theta \cos\varphi$, $Y = r \sin\theta \sin\varphi$, $Z = r \cos\theta$, which is such that $Z_{\theta} = -r \sin\theta$.  Hence the negative root is the reasonable one here, for we wish our procedure to give the usual result for metrics which are flat.}
\begin{equation}
\nonumber Z_{\theta} \doteq - \sqrt{g_{22} + T_{\theta}^2 - \rho_{\theta}^2}.
\end{equation}
This leaves the remaining unknowns $\{ T_{t}, \rho_{t}, Z_{t}, \rho_{r}, Z_{r} \}$ to be determined.
Consider next the pair (\ref{g11}, \ref{g12}):
\begin{eqnarray}
\rho_{r} \rho_{\theta} + Z_{r} Z_{\theta} &\doteq& T_{r} T_{\theta} + g_{12} = T_{r} T_{\theta},
\\
\rho_{r}^2 + Z_{r}^2 &\doteq& g_{11} + T_{r}^2.
\end{eqnarray}
One can first get $Z_{r}$,
\begin{eqnarray}
Z_{r} \doteq \frac{T_{r} T_{\theta} - \rho_{r} \rho_{\theta}}{Z_{\theta}}, \label{Zr}
\end{eqnarray}
and then a quadratic equation for $\rho_{r}$, namely
\begin{eqnarray}
\rho_{r}^2 \left( 1 + \frac{\rho_{\theta}^2}{Z_{\theta}^2} \right) - 2 \frac{\rho_{r} \rho_{\theta} T_{r} T_{\theta}}{Z_{\theta}^2} + \frac{T_{r}^2 T_{\theta}^2}{Z_{\theta}^2} - T_{r}^2 - g_{11} \doteq 0\,.
\end{eqnarray}
Generally $\rho_{r}$ has two solutions,
but one should select the sign that gives a positive value---at least when the dynamical geometry is near flat space.
Once $\rho_{r}$ is determined, $Z_{r}$ can be found from (\ref{Zr}).

As to $T_{t}, \rho_{t}, Z_{t}$, they can now be found through the remaining relations (\ref{g00}, \ref{g01}, \ref{g02}) which are effectively one quadratic and two linear equations:
\begin{eqnarray}
-T_{t}^2 + \rho_{t}^2 + Z_{t}^2 &\doteq& g_{00} - \rho^2 \Phi_{t}^2, \label{00}
\\
-T_{t} T_{r} + \rho_{t} \rho_{r} + Z_{t} Z_{r} &\doteq& g_{01} = 0,
\\
-T_{t} T_{\theta} + \rho_{t} \rho_{\theta} + Z_{t} Z_{\theta} &\doteq& g_{02} = 0.
\end{eqnarray}
From the two linear equations we find
\begin{eqnarray}
\rho_{t} &\doteq& T_{t} \, \frac{T_{r} Z_{\theta} - T_{\theta} Z_{r}}{\rho_{r} Z_{\theta} - \rho_{\theta} Z_{r}}, \label{rhot}
\\
Z_{t} &\doteq& T_t\, \frac{\rho_r T_{\theta} - T_r \rho_{\theta}}{\rho_{r} Z_{\theta} - \rho_{\theta} Z_{r}}. \label{Zt}
\end{eqnarray}
Substituting this into (\ref{00}) gives
a new quadratic equation for $T_{t}$:
\begin{eqnarray}
T_{t}^2 \left[ -1 + \frac{(T_{r} Z_{\theta} - T_{\theta} Z_{r})^2}{(\rho_{r} Z_{\theta} - \rho_{\theta} Z_{r})^2} + \frac{(\rho_{r} T_{\theta} - \rho_{\theta} T_{r})^2}{(\rho_{r} Z_{\theta} - \rho_{\theta} Z_{r})^2} \right] \doteq g_{00} - \rho^2 \Phi_{t}^2.
\end{eqnarray}
This has two solutions, but one should choose the positive value.
Once $T_{t}$ is found, $\rho_{t}$ and $Z_{t}$ can be found from (\ref{rhot}, \ref{Zt}).

Algebraically solving these equations and collecting our isometric matching results together, we have in addition to vanishing $\Phi_r,\ \Phi_\theta$,
\begin{eqnarray}
\rho &\doteq& \sqrt{g_{33}}, \\
\Phi_{t} &\doteq& g_{03}/g_{33}, \\
Z_{\theta} &\doteq& - L_{1}, \\
T_{t} &\doteq& L_{2} L_{3}, \\
\rho_{r} &\doteq& (\rho_{\theta} T_{r} T_{\theta} + {L_{1} L_{3}})/{L_{4}}, \\
Z_{r} &\doteq& ({\rho_{\theta} L_3} - {L_{1} T_{r} T_{\theta}})/{L_{4}}, \\
\rho_{t} &\doteq& L_2 ({L_{1} T_{r} g_{22}} + {L_3 T_{\theta} \rho_{\theta}})/{L_{4}}, \\
Z_{t} &\doteq& L_2 ({\rho_{\theta} T_{r} g_{22}} - {L_{1} L_3 T_{\theta}})/{L_{4}},
\end{eqnarray}
where
\begin{eqnarray}
L_{1} &:=& \sqrt{g_{22} - \rho_{\theta}^2 + T_{\theta}^2}, \\
L_{2} &:=& \sqrt{\frac{g_{03}^2 - g_{00} g_{33}}{g_{11} g_{22} g_{33}}}, \\
L_{3} &:=& \sqrt{T_{r}^2 g_{22} + g_{11} g_{22} + g_{11} T_{\theta}^2}, \\
L_{4} &:=& g_{22} + T_{\theta}^2.
\end{eqnarray}

With 4D isometric matching on the 2-boundary $S$ the second term in our quasi-local expression (\ref{B}) vanishes, so it simplifies considerably to
\begin{equation}
16 \pi \mathcal{B}(N) = (\Gamma^\alpha{}_\beta - \bar\Gamma^\alpha{}_\beta) \wedge i_{N} \eta_{\alpha}{}^{\beta}. \label{B2}
\end{equation}
We will use this expression to determine our quasi-local energy and angular momentum.

\subsection{Energy extremization}
Even with 4D isometric matching the quasi-local values still have a lot of freedom. The value of the energy-momentum can be regarded as a measure of the difference between the dynamical geometry and the reference geometry.  With this in mind it was proposed to look to the critical points of $m^2 := - \bar g^{ij} p_i p_j$~\cite{ae100}.  Here we shall use a simpler expression which should lead to the same answer.
To properly select  $T_{r}$ and $T_{\theta}$ we propose extremizing the quasi-local energy to find the best matching for these embedding functions.

The quasi-local energy is determined by selecting the reference geometry unit timelike Killing vector\footnote{The subscript in $N_{\rm E}$ is used
to denote the vector field for energy}: $N_{\rm E} = \partial_{T}$.
To calculate the quasi-local energy one needs  the components of $N_{\rm E}$ and $\bar\Gamma^\alpha{}_\beta$ in the dynamic spacetime.
The reference connection is determined by pullback from the flat space reference (which vanishes in Minkowski coordinates):
\begin{equation}
\bar \Gamma^\alpha{}_\beta = x^\alpha{}_i (\bar\Gamma^i_{~j} y^j{}_\beta + dy^i{}_\beta) = x^\alpha{}_i dy^i{}_\beta.
\end{equation}
Here $x^\alpha{}_i$ is the inverse of $y^i{}_\alpha$, and the latter is determined from $dy^i = y^i{}_\alpha dx^\alpha$, i.e.,
\begin{eqnarray} \label{reference choice}
\left(
  \begin{array}{c}
    dT \\
    dX \\
    dY \\
    dZ \\
  \end{array}
\right)=\left(
          \begin{array}{cccc}
            T_t & T_r & T_\theta & 0 \\
            \rho_t \cos\xi - \rho \sin\xi \Phi_t\ \ & \rho_r \cos\xi\ \ & \rho_\theta \cos\xi\ \ & -\rho \sin\xi \\
            \rho_t \sin\xi + \rho \cos\xi \Phi_t\ \ & \rho_r \sin\xi\ \ & \rho_\theta \sin\xi\ \ & \rho \cos\xi \\
            Z_t & Z_r & Z_\theta & 0 \\
          \end{array}
        \right)\left(
                 \begin{array}{c}
                   dt \\
                   dr \\
                   d\theta \\
                   d\varphi \\
                 \end{array}
               \right),
\end{eqnarray}
where $\xi = \varphi + \Phi$.

Using these results, the determined values for the components of the vector are
\begin{eqnarray}\label{eq46}
&& N^{t}_{\rm E} \doteq \frac{L_3}{g_{11} g_{22} L_2}, \qquad N^{r}_{\rm E} \doteq - \frac{T_{r}}{g_{11}}, \qquad N^{\theta}_{\rm E} \doteq -\frac{T_{\theta}}{g_{22}},
\nonumber\\
&& N^{\varphi}_{\rm E} \doteq - \Phi_t N^t_{\rm E} = - \frac{g_{03} L_3}{g_{11} g_{22} g_{33} L_2}. \label{N}
\end{eqnarray}
With $N = N_{\rm E}$ the quasi-local boundary term (\ref{B2}) works out to be
\begin{eqnarray}
16 \pi \mathcal{B}(N_{\rm E}) \!&\!\!=\!\!&\! \sqrt{g_{33}} \biggl\{ \frac{2 T_{\theta}}{L_{3}} T_{r\theta} + \frac{2}{L_4} \left[ -\frac{T_{\theta}^2 T_{r}}{L_{3}} + \frac{T_{\theta} \rho_{\theta}}{L_{1}} \right] T_{\theta\theta}
- \frac{g_{11\theta} T_r T_\theta}{L_3 g_{11}}
\nonumber\\
&& \quad + \frac{g_{22\theta}}{L_4} \left[ \frac{T_r T_\theta^3}{g_{22} L_3} + \frac{\rho_\theta}{L_1} \right]
- \frac{2 \rho_{\theta\theta}}{L_{1}} + \frac{2 L_{1}}{\sqrt{g_{33}}} + L_{3} W \biggr\} d\theta \wedge d\varphi, \qquad \label{BE}
\end{eqnarray}
with
\begin{equation}
W := -\frac{1}{g_{11}} \left( \frac{g_{22r}}{g_{22}} + \frac{g_{33r}}{g_{33}} \right).
\end{equation}
Here $g_{22\theta}$ denotes its derivative by $\theta$ and so on.  Also it should be noted that in this expression and in all of the following, to convey our ideas without too many complications, we have additionally assumed that the  metric components in (\ref{metric}) are all time independent.

The quasi-local energy value is the integral of the 2-form $\mathcal{B}(\partial_T)$ (\ref{BE}) over $S$.  The value is a functional of the two control functions $T_\theta, T_r$.  The distinguished values of this energy functional are the ``critical points'', which can be found as the solutions of  the respective variational derivatives, the two Euler-Lagrange equations obtained by variation with respect to the control variables $T_\theta,T_r$.  The respective variational equations are
\begin{eqnarray}\label{eq50}
0 &=& \frac{\partial \mathcal{B}(\partial_T)}{\partial T_{\theta}} - \frac{d}{d\theta} \frac{\partial \mathcal{B}(\partial_T)}{\partial T_{\theta\theta}}
\nonumber\\
&=& \frac{\sqrt{g_{33}}}{L_3} \left[ 2 T_{r\theta} + \left( -\frac{g_{11\theta}}{g_{11}} + \frac{g_{33\theta} T_{\theta}^2}{g_{33} L_{4}}\right) T_{r} + \frac{2 L_{1} L_3 T_{\theta}}{\sqrt{g_{33}} L_{4}} + g_{11} W T_{\theta} \right], \label{ELtheta}
\\
\label{eq51}
0 &=& \frac{\partial \mathcal{B}(\partial_T)}{\partial T_{r}} - \frac{d}{d\theta} \frac{\partial \mathcal{B}(\partial_T)}{\partial T_{r\theta}}
\nonumber\\
&=& \frac{\sqrt{g_{33}}}{L_3} \left[ -2 T_{\theta\theta} + \left( \frac{g_{22\theta}}{g_{22}} - \frac{g_{33\theta}}{g_{33}} \right) T_{\theta} + g_{22} W T_{r} \right]. \label{ELr}
\end{eqnarray}
The second relation can be used to eliminate $T_r$ from the first relation.  The result is a quasi-linear (but highly nonlinear) second order equation for $T_\theta$.  For such equations there is generally little hope of finding any analytic solution.  However, this case is still rather special: there is one quite obvious simple solution that applies to all of these Kerr-like metrics, namely the trivial solution $T_r = 0 = T_\theta$.

For the simple solution $T_r=0=T_\theta$ the vector field takes the simpler form
\begin{eqnarray}\label{eq58}
N^{t}_{\rm E} = \frac{1}{T_{t}} = \sqrt{\frac{g_{33}}{g_{03}^2 - g_{00}g_{33}}}, \quad N^{\varphi}_{\rm E} = -\frac{g_{03}}{g_{33}} N^t_{\rm E}, \quad N^{r}_{\rm E} = 0, \quad N^{\theta}_{\rm E} = 0.\label{vector}
\end{eqnarray}
This vector is orthogonal to the constant $t$ hypersurface; physically it can be identified with the \emph{static\/} observer on the boundary.

From (\ref{BE}) with $T_\theta=0=T_r$, on the $t,r$ constant 2-surface of any axisymmetric time independent metric of the form (\ref{eq20}), our optimal boundary term for the quasi-local energy is
\begin{eqnarray}\label{eq4}
\mathcal{B}_{\rm E} &=& \frac{\sqrt{g_{33}}}{16\pi} \Bigl\{ -\Bigl( \frac{\partial_{r} g_{22}}{g_{22}} + \frac{\partial_{r} g_{33}}{g_{33}} \Bigr) \frac{\sqrt{g_{22}}}{\sqrt{g_{11}}}
\nonumber\\
&+& \frac{(\partial_{\theta} g_{22})(\partial_{\theta} g_{33}) - 2 g_{22} (\partial_{\theta\theta} g_{33}) + 4 g_{22}^2}{g_{22} \sqrt{4 g_{33} g_{22} - (\partial_{\theta} g_{33})^2}} ~\Bigr\} d\theta \wedge d\varphi. \label{optimalenergy}
\end{eqnarray}

Turning to the quasi-local angular momentum, note that it can have only one non-vanishing component.  To select this component from our quasi-local boundary term (\ref{B2}) the appropriate vector field is the \emph{reference rotational Killing field}, $N_{\rm AM} = X \partial_{Y} - Y \partial_{X}$, around the axis. In the dynamic spacetime, this vector field works out to have a very simple form, namely
\begin{eqnarray}\label{eq54}
N_{\rm AM} = \partial_{\varphi},
\end{eqnarray}
and it is in fact the \emph{axisymmetric Killing field} for the metrics considered.
For the choice of $N_{\rm AM} = \partial_{\varphi}$, which implies that the only non-vanishing vector component is $N^{\varphi} = 1$, the boundary expression (\ref{B2}) takes the specific form
\begin{eqnarray}\label{AM N3=1}
16 \pi \mathcal{B}_{\rm AM} = \sqrt{-g} (g^{1\beta} \Delta\Gamma^{0}{}_{\beta3} - g^{0\beta} \Delta\Gamma^{1}{}_{\beta 3} ) d\theta \wedge d\varphi.
\end{eqnarray}
Using the reference choice (\ref{reference choice}), a computation shows that the reference connection contribution vanishes, so the boundary term depends only on the dynamical part:
\begin{eqnarray}\label{AM N3=1 dyn}
16 \pi \mathcal{B}_{\rm AM} &=& \sqrt{-g} (g^{11} \Gamma^{0}{}_{13} - g^{00} \Gamma^{1}{}_{03} - g^{03} \Gamma^{1}{}_{33}) d\theta \wedge d\varphi
\nonumber\\
&=& \sqrt{-g} g^{11} (g^{00} \partial_{r} g_{03} + g^{03} \partial_{r} g_{33}) d\theta \wedge d\varphi. \label{J}
\end{eqnarray}

\section{Detailed results}

To see whether this simple choice is really a good choice---maybe even the best---for determining the reference we examined some physically meaningful special cases.

\subsection{Minkowski metric}
For the Minkowski metric it is easy to see that there is actually a simple non-trivial solution to the two embedding equations (\ref{eq50}) and (\ref{eq51}), namely
\begin{equation}
T_{r} = -\lambda \cos\theta, \qquad T_{\theta} = \lambda r \sin\theta , \label{boostedref}
\end{equation}
with a constant $\lambda$. With this choice the vector field (\ref{N}) corresponds to that of a constant boost along the axis:
\begin{eqnarray}\label{eq53}
N^{t}_{\rm E} = \sqrt{\lambda^2 + 1}, \quad N^{r}_{\rm E} = \lambda \cos\theta, \quad N^{\theta}_{\rm E} = -\frac{\lambda\sin\theta}{r}, \quad N^{\varphi}_{\rm E} = 0.
\end{eqnarray}
For all $\lambda$ the value of the optimal boundary expression (\ref{BE}) and the associated energy vanish:
\begin{eqnarray}\label{eq6}
\mathcal{B}(N_{\rm E}) = 0, \qquad {E} = \oint \mathcal{B}(N_{\rm E}) = 0.
\end{eqnarray}
This is sensible: a sphere in Minkowski space has vanishing energy, and for a boosted observer it still has vanishing energy.  It should be noted that, because we constrained everything to be axisymmetric, we did not find general boost freedom---only boost freedom along the axis.
Because the Minkowski geometry has this symmetry, naturally our procedure cannot find a unique reference.  Furthermore, because the symmetries of Minkowski are  known we need not look further for any other solutions to our equations to find a good reference.

It is easy to verify that, just as it should, for any $\lambda$ our expression also gives vanishing angular momentum for Minkowski space.

\subsection{Schwarzschild metric}
We next consider the Schwarzschild metric.  To identify the desired physical solution, consider first the limit for small $m$ and/or large $r$; then the physical solution should approach the one just found for Minkowski space.  But expressions of the form (\ref{boostedref}) can only solve the equations (\ref{ELtheta},\ref{ELr}) for the Schwarzschild metric for the trivial value $\lambda=0$.  This is the physical solution we want to use for our preferred reference, both for this limit and for more general values of $m$ and $r$.  With vanishing $T_r,T_\theta$, from (\ref{vector}) for the Schwarzschild metric is that of a static observer:
\begin{eqnarray}\label{eq57}
N^{t}_{\rm E} = 1/\sqrt{-g_{00}},
\quad N^{\varphi}_{\rm E} = 0, \quad N^{r}_{\rm E} = 0, \quad N^{\theta}_{\rm E} = 0,
\end{eqnarray}
which is only meaningful for $-g_{00} = 1 - 2m/r > 0$, i.e., outside the horizon.  This strategy cannot be used to find a reference for the dynamic region inside the horizon (alternatives that can cope with the interior region have been discussed elsewhere~\cite{Wu10,Wu11,Wu12}).
For this case the energy boundary term (\ref{optimalenergy}) and the associated energy values are
\begin{eqnarray}
&& \mathcal{B}_{\rm E} = \frac{4r}{16\pi} (1 - \sqrt{1 - {2m}/{r}}) \sin\theta d\theta \wedge d\varphi,
\nonumber\\
\label{eq8}
&& {E} = \oint \mathcal{B}_{\rm E} = r (1 - \sqrt{1 - {2m}/{r}}) = \frac{2m}{1 + \sqrt{1-{2m}/{r}}}.
\end{eqnarray}
This is the well known quasi-local energy first found by Brown and York~\cite{BY93}.

From (\ref{J}) it is easy to see that, as it should, our quasi-local angular momentum vanishes for the Schwarzschild metric.

\subsection{Kerr metric}
Let us now consider specifically the Kerr metric. In terms of the convenient Boyer-Linquist coordinates, the Kerr solution is
\begin{eqnarray}
ds^2 &=& - \frac{(\Delta - a^2 \sin^2\theta)}{\Sigma} dt^2 + \frac{\Sigma}{\Delta} dr^2 + \Sigma d\theta^2
\nonumber\\
&&-\frac{4 m a r \sin^2\theta}{\Sigma} dt d\varphi + \left[ \frac{(r^2 + a^2)^2 - \Delta a^2 \sin^2\theta}{\Sigma} \right] \sin^2\theta d\varphi^2, \label{Kerr}
\end{eqnarray}
where $\Sigma(r, \theta) := r^2 + a^2 \cos^2\theta$, $\Delta(r) := r^2 - 2 m r + a^2$.
For the Kerr metric the special relations $g_{22}/g_{11} = \Delta$ and $g_{03}^2 - g_{00} g_{33} = \Delta\sin^2\theta$ hold.

By the same argument used for the Schwarzschild metric, we take as our desired reference the trivial solution $T_\theta = 0 = T_r$.  The associated displacement vector field from (\ref{vector},\ref{Kerr}) has a well-known form.
For the Kerr metric the expression (\ref{optimalenergy}) is too complicated to find an analytical result.  So we turn to other procedures.  We first find an approximate quasi-local energy for a slowly rotating Kerr black hole, then we present a numerical evaluation.

\subsubsection{When $|a| \ll m$}
With the Kerr metric values the slow rotation approximation is
\begin{eqnarray}\label{eq12}
\mathcal{B}_{\rm E} &\approx& \frac{1}{16\pi} \Biggl\{ 4 r \Bigl( 1 - \sqrt{1 - \frac{2m}{r} + \frac{a^2}{r^2}} \Bigr) \sin\theta
\nonumber\\
&& + \Bigl[ \Bigl( \frac{2 \sin\theta \cos^2\theta}{r} + \frac{2 m \sin^3\theta}{r^2} \Bigr) \sqrt{1 - \frac{2m}{r} + \frac{a^2}{r^2}}
\nonumber\\
\label{eq13}
&& + \frac{2 \sin^3\theta}{r} - \frac{20 m \sin\theta \cos^2\theta}{r^2} + \frac{8 m \sin\theta}{r^2} \Bigr] a^2 \Biggr\} d\theta \wedge d\varphi,
\nonumber\\
\label{eq14}
{E} &=& \oint \mathcal{B}_{\rm E}
\nonumber\\
&\approx& \!\! r \Bigl( 1 - \sqrt{1 - \frac{2m}{r} + \frac{a^2}{r^2}} \Bigr) + \frac{(1 + \frac{2m}{r}) \sqrt{1 - \frac{2m}{r} + \frac{a^2}{r^2}} + 2 + \frac{2m}{r}}{6r} a^2. \qquad
\end{eqnarray}
Our quasi-local energy value for the Kerr metric using the slow rotation approximation turns out to be identical with that found by Martinez from the Brown-York quasi-local expression~\cite{Mart94}.

\subsubsection{General cases: $0\leq|a|\leq m$}
On the $t,r$ constant closed 2-surface we calculated numerically the value of our quasi-local energy with
our best matched reference for the Kerr metric.
In Figure~1 we show how the  quasi-local energy in units of $m$ (i.e., $E/m$) depends on $a/m$ and $r/m$.  Some of the numerical data is presented in Table 1.
\begin{figure}[pb]
\begin{center}
  \includegraphics[trim=20 20 20 20, clip,width=3in,angle=270]{3d}\\
  \vspace*{8pt}
  \caption{A plot of our Kerr quasi-local energy for $r\geq2$, with all the quantities given in units of $m$.  No matter what the value of $a\le 1$, as $r$ approaches infinity, $E$  approaches $1$, as is well known. The highest point corresponds to the Schwarzschild case, with $E=2$ at $r=2$. For any $r$, the larger $|a|$ is, the smaller $E$ is. }
  \end{center}
\end{figure}

\begin{table}
\caption{The quasi-local energy ${E}$ dependance on $|a|$ and  $r$.}
{\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
  ${E}/m$&$r/m=2$&$r/m=3$&$r/m=10$&$r/m=20$&$r/m=50$\\\hline
  $|a|/m=0.0$ & 2.00000000 & 1.26794919 & 1.05572809 & 1.0263340 & 1.0102052  \\\hline
  $|a|/m=0.1$ & 1.90258182 & 1.26708057 & 1.05571462 & 1.0263325 & 1.0102050  \\\hline
  $|a|/m=0.2$ & 1.81063676 & 1.26449931 & 1.05567426 & 1.0263278 & 1.0102048  \\\hline
  $|a|/m=0.3$ & 1.72457134 & 1.26027757 & 1.05560706 & 1.0263200 & 1.0102042  \\\hline
  $|a|/m=0.4$ & 1.64468939 & 1.25453045 & 1.05551316 & 1.0263091 & 1.0102036  \\\hline
  $|a|/m=0.5$ & 1.57117096 & 1.24740880 & 1.05539271 & 1.0262951 & 1.0102028  \\\hline
  $|a|/m=0.6$ & 1.50406598 & 1.23909055 & 1.05524594 & 1.0262779 & 1.0102017  \\\hline
  $|a|/m=0.7$ & 1.44330154 & 1.22977111 & 1.05507309 & 1.0262577 & 1.0102004  \\\hline
  $|a|/m=0.8$ & 1.38869944 & 1.21965415 & 1.05487447 & 1.0262344 & 1.0101990  \\\hline
  $|a|/m=0.9$ & 1.33999980 & 1.20894300 & 1.05465044 & 1.0262080 & 1.0101974  \\\hline
  $|a|/m=1.0$ & 1.29688677 & 1.19783368 & 1.05440138 & 1.0261786 & 1.0101956  \\\hline
\end{tabular}}
\end{table}

\subsubsection{Kerr angular momentum}

Evaluating (\ref{J}) explicitly for the Kerr metric (\ref{Kerr}) we find that the exact boundary expression for our quasi-local angular momentum on the $t,r$  constant 2-surface works out to
\begin{eqnarray}
16 \pi \mathcal{B}_{\rm AM} &=& \frac{2 m a \sin^{3}\theta}{\Sigma} [ (1 - 2 r^{2}/\Sigma) (r^{2} + a^{2}) - 2 r^{2} ]\, d\theta \wedge d\varphi
\nonumber\\
&=& - \frac{2 m a [ 3 r^4 + r^2 a^2 (1 + \cos^2\theta) - a^4 \cos^2\theta ] {\sin^3\theta}}{(r^2 + a^2 \cos^2\theta)^2} ~d\theta \wedge d\varphi. \qquad \label{KerrBAM}
\end{eqnarray}
Hence the one non-vanishing component of our quasi-local angular momentum for the Kerr metric is
\begin{eqnarray}
\label{eq11} {J_\varphi} = \oint \mathcal{B}_{\rm AM} = \frac{m a \cos\theta (3 r^2 - r^2 \cos^2\theta + a^2 + a^2 \cos^2\theta)}{4 (r^2 + a^2 \cos^2\theta)} \Big|^{\pi}_{0} = - m a.
\end{eqnarray}

It is amazing that our quasi-local angular momentum turns out to have a simple natural $r$-independent value, even though the boundary term ${\cal B}(\partial_\varphi)$ (\ref{KerrBAM}) is a not so simple function of $r$.

\section{Conclusion}

Our approach to identifying the quasi-local energy and angular momentum for gravitating systems is via the covariant Hamiltonian formalism.
How to choose the reference geometry and the evolution vector that appear in the Hamiltonian boundary term have been outstanding issues in the quasi-local program. Minkowski geometry and its associated Killing fields are the natural choice, but one actually needs a program for selecting a specific Minkowski geometry on the boundary.
Here for axisymmetric metrics we show that this can be accomplished using an appealing program: \emph{four}-dimensional isometric matching on the boundary 2-surface and energy optimization.

Exact 4D isometric matching includes 2-surface isometric matching.  In general this is the problem of embedding a closed 2-surface into Minkowski space.  Wang and Yau have carefully considered this problem in connection with their quasi-local energy investigations~\cite{WY09PRL,WY09CMP}.  The problem generally involves the solution of a nonlinear partial differential equation.  Although they proved an existence and uniqueness theorem, it does not seem possible to find a general analytic formula for the embedding functions.  However, for the axisymmetric case  considered here the situation is quite different.
The 2-surface isometric matching, and indeed the whole 4D isometric matching is in this case essentially a simple \emph{algebraic} problem. With 4D isometric matching the choice of Minkowski reference on the boundary still has two degrees of freedom.

To determine the remaining two degrees of freedom we have proposed energy optimization.
Our energy optimization procedure in the axisymmetric case leads to a pair of first order quasi-linear \emph{ordinary\/} differential equations.
A certain class of Kerr-like metrics fortunately have a simple analytic solution (the trivial solution) to these non-linear energy optimization equations.  One cannot expect to find such simple analytic solutions in more general cases, such as a dynamical metric with more non-vanishing cross terms, a time dependent metric, or for some other 2-surface $S$ described by a formula like $R = R(r,\theta) = \hbox{constant}$.  That we can have our solution of trivial embedding variables can be attributable to the fact that the Boyer-Lindquist coordinate $r$ has a special significance for these Kerr-like geometries, and the constant $r$ surfaces are thus quite special and meaningful.

For the considered metrics it is noteworthy that our optimal timelike vector field turned out to be hypersurface orthogonal, corresponding to a static observer. Our quasi-local expression with  best matched reference gives for Minkowski space, as expected, vanishing energy and angular momentum.  For the Schwarzschild metric it gives the well-known energy value first found by Brown and York along with vanishing angular momentum.

We found the exact analytic formula for our optimal quasi-local energy boundary term for these Kerr-like metrics on the surface $r=\hbox{constant}$.
For the Kerr metric in the slow rotation approximation our quasi-local energy agrees with that  found by Martinez using the Brown-York quasi-local expression.
We plotted our Kerr quasi-local energy and note that it is a decreasing function of $|a|$.  At present we have no physical understanding as to why this is the case.

For these Kerr-like metrics we found our optimal quasi-local angular momentum boundary term.  Our quasi-local angular momentum for the Kerr metric worked out to be an unexpectedly simple but satisfying constant value ($-ma$).
This may again be attributed to the chosen 2-surface $S$, which is given by a constant value of the Boyer-Lindquist coordinate $r$.  It is another indication of the special significance of this coordinate for these metrics.  Our quasi-local angular momentum would not be so simple for closed surfaces of the form $R=R(r,\theta)=\hbox{constant}$.

Encouraged by these results, the application of these ideas is being extended to both general axisymmetric metrics and completely general metrics.

\section*{Acknowledgments}
This work was supported by the National Science Council of the R.O.C. under the grants NSC 101-2112-M-008-006 and 99-2112-M-008-005-MY3, the Institute of Physics, Academia Sinica, and the National Center of Theoretical Sciences (NCTS).


\section{Introduction}
Traditionally, power distribution systems have been engineered as radial networks that carry power from the utility feeder to the end consumers of electricity. These consumers are seen as a set of passive loads, and generally, the utility has no control of their electricity demand, and does not collect any real-time information from them. In the future, under the \emph{SmartGrid} initiative, distribution systems are envisioned to undergo dramatic transformation in structure and functionality. Instead of passive demands, buses in the network could include clusters of households with rooftop solar installations, electric vehicles, and storage; these generation and storage devices are commonly referred to as distributed energy resources(DERs). It has been shown that massive penetration of rooftop solar can cause voltages to rise in the distribution system, while large-scale deployment of electrical vehicles can cause voltages to drop \cite{Carvalho08,Guille09}. The objective of this paper is to address this voltage regulation problem.

Voltage regulation in distribution networks is traditionally achieved via tap-changing under-load transformers and switched capacitors; however, while ensuring that voltage requirements are met, the setting of these devices might not optimal, in the sense that they may not minimize the thermal loss or some other objectives. On the other hand, by relying on a cyber-infrastructure for sensing, communication, and control, it is also possible to utilize DERs with capability to provide reactive power for optimal voltage regulation (see, e.g., \cite{Lam12} and the references therein). Unfortunately, the current speed of installation of DERs far outpaces the deployment of the cyber-infrastructure.\footnote{The advanced metering infrastructure (e.g. smart meters) transmit their information on a daily bases and are not used for real-time communication.} For example, most distribution networks in California have not implemented real-time two-way communication between DERs and the utilities, but the state already has 2900 MW of installed capacity \cite{SEIA13}, enough for about 600,000 households. Therefore we are faced with a situation where DERs are already in the distribution network with no communication network to coordinate them. Thus, an urgent question arises: in order to utilize DERs for voltage regulation in distribution networks, how should they be managed with no communication or very limited communication between the buses in the network?

In our setting, we assume each bus can measure its own voltage magnitude, active and reactive power injections, but cannot communicate with other buses in the network. The available control action at a bus is the adjustment of its active and reactive power injections. Altering the active power injection at certain buses corresponds to demand response (DR) actions. Because of economical and regulatory considerations, dynamic DR is difficult to implement without a communication network; therefore, in this paper we only consider reactive power injections at each bus as the main control knob for DER-based voltage regulation. We assume that reactive power can be provided by, e.g., the power electronics grid interfaces of rooftop solar installations \cite{JoOo:00,LoKr:01,solar_bridge}.

 In this paper, we propose a method for voltage regulation that relies on the local control actions\footnote{We call the control actions \emph{local control actions} to emphasize that only local information is available at each bus.} of reactive-power-capable DERs. We provide sufficient conditions such that local control is guaranteed to be successful in maintaining any specified voltage profile. We also provide conditions on the network topology and active power injection regions that shows when local control cannot maintain voltage stability in the system. We arrive at our results by casting the voltage regulation problem as an optimization problem, and investigate the conditions under which optimality can be achieved without communication between the buses.

The problem of voltage regulation in distribution networks has received considerable attention in recent years. Centralized algorithms have been proposed in \cite{Baran07,Turitsyn10,Villacci06} and decentralized algorithms have been proposed in \cite{Farivar12,Lam12,Bolognani12,Robbins12}. For the latter, communication is generally assumed to be possible between subsets of the buses (e.g., \cite{Lam12} provides a distributed optimization algorithm that adheres to neighbor-to-neighbor communications, and \cite{Bolognani12} can choose different communication topology using an approximate model for power flow). In contrast, our proposed algorithm does not rely on a communication infrastructure. A related questions in context of voltage stability has been studied in \cite{Ilic:1986}.

 This paper is organized as follows.  Section \ref{sec:model} introduces the notations and the power system distribution model adopted throughout the paper. Section \ref{sec:result} gives the main algorithm. Section \ref{sec:sr} provides convergence analysis of the algorithm. Section \ref{sec:sim} provides case studies,  and Section \ref{sec:con} concludes the paper. The proofs of the main results are provided in the Appendix.

\section{Model and Problem Formulation} \label{sec:model}
This section introduces the power distribution system models and also provides the formulation of the problem. Consider a distribution network with $n$ buses. Today, most power distribution systems operate as radial networks connected to a single feeder\cite{Kersting06}; therefore we model the network as a connected tree rooted at the feeder. By convention, the feeder is labelled as bus $0$. The distance from bus $i$ to the root (bus 0) is the number of vertices on the unique path from bus $i$ to the root. The depth of the tree is the maximum of the bus distances to the root.

Let $V_i=|V_i| e (j \theta_i)$ denote bus $i$ complex voltage, and define  vector $\mathbf v =[ V_0 \; V_1 \dots V_{n-1}]^T \in \mathbb{C}^n$. Let $I_i$, $P_i$ and $Q_i$ denote, the complex current, active and reactive power injections at bus $i$, respectively; and define $\mathbf i=[I_0 \; \dots \; I_{n-1}] \in \mathbb{C}^n$, $\mathbf p=[P_0 \; \dots \; P_{n-1}] \in \mathbb{R}^n$ ,and $\mathbf q=[Q_0 \; \dots \; Q_{n-1}] \in \mathbb{R}^n$. Let $\mathbf s=\mathbf p + j \mathbf q$ denote the complex power. The active powers in a distribution network are typically constrained by upper and lower capacity limits, i.e.,  of $\underline{P}_i \leq P_i \leq \overline{P}_i$, where the upper and lower bound are determined by by the types of devices at bus $i$. Let $\mathcal P =\{\mathbf p \in \mathbb{R}^{n-1}: \underline{P}_i \leq P_i \leq \overline{P}_i \}$ be the feasible region formed by the bus injections (see \cite{Zhang13,LTZ12} for more detailed discussions).

 Given two connected buses $i$ and $k$, let $y_{ik}=g_{ik}-jb_{ik}$ denote the admittance of the transmission line between them, and let $\mathbf Y$ be the bus admittance matrix \cite{Overbye04}. Then, the complex power injections are related to the voltage by $\mathbf s = \diag (\mathbf v \mathbf v^H \mathbf Y^H)$,
 where $(\cdot)^H$ denotes the Hermitian transpose, and the $\diag(\cdot)$ operator returns the diagonal if the argument is a square matrix or a diagonal matrix if the argument is a vector. We are interested in regulating the voltage in the distribution network, so let $V_0^{ref}, V_1^{ref},\dots,V_{n-1}^{ref}$ be the set of reference voltage magnitudes that we are interested in tracking. We assume that the feeder behaves as an infinite bus, thus $V_0=V_0^{ref}$. We assume that bus $i$ can control its reactive power, whereas the active power is set at a prescribed value $P_i^*$. We formulate the voltage regulation problem as the following feasibility problem:
 \begin{subequations} \label{eqn:feasible}
 \begin{align}
 \mbox{find } & \mathbf q \\
 \mbox{s.t. } & |V_i|=V_i^{ref},~ i=0,\dots,n-1 \\
 &  P_i=P_i^*, \; i=1,\dots,n-1 \label{eqn:p_f} \\
  & \mathbf p + j \mathbf q =\diag(\mathbf v \mathbf v^H \mathbf Y^H), \label{eqn:s_f}
\end{align}
\end{subequations}
 where \eqref{eqn:p_f} corresponds to the fact that active powers are fixed, and the equality constraints in (1d) correspond to the power flow equations.
That is, we seek to find a reactive power injection vector $\mathbf q$ such that the voltages are at their desired levels. The feasibility of \eqref{eqn:feasible} can be decided by the following optimization problem:
 \begin{subequations} \label{eqn:main_v}
 \begin{align}
 \min_{\mathbf q }\; & \sum_{i=1}^{n-1} (|V_i|-V_i^{ref})^2 \label{eqn:obj_v} \\
 \mbox{s.t. } & P_i=P_i^*, \; i=1,\dots,n-1 \label{eqn:p_v} \\
 & V_0=V_0^{ref} \label{eqn:v0_v} \\
 & \mathbf p + j \mathbf q =\diag(\mathbf v \mathbf v^H \mathbf Y^H), \label{eqn:s_v}
 \end{align}
 \end{subequations}
 where the feeder bus (bus 0) always holds its voltage to be at $V_0^{ref}$. The objective function means that the system is trying to minimize the deviation of the voltage magnitudes to their reference values. The feasibility problem in \eqref{eqn:feasible} has a solution if and only if the optimal value of the optimization problem \eqref{eqn:main_v} is 0.

 It turns out the objective function in \eqref{eqn:obj_v} is not easy to analyze. We rewrite \eqref{eqn:main_v} with a new objective function as:
  \begin{subequations} \label{eqn:main}
  \begin{align}
  \min_{\mathbf q }\; & \sum_{i=1}^{n-1} \left(|V_i|^2-(V_i^{ref})^2 \right)^2 \label{eqn:obj} \\
  \mbox{s.t. } & P_i=P_i^*, \; i=1,\dots,n-1 \label{eqn:p} \\
  & V_0=V_0^{ref} \label{eqn:v0} \\
  & \mathbf p + j \mathbf q =\diag(\mathbf v \mathbf v^H \mathbf Y^H),  \label{eqn:s}
  \end{align}
  \end{subequations}
 Intuitively, the reason that \eqref{eqn:obj} is much easier to handle than \eqref{eqn:obj_v} is that $|V_i|$ is not differentiable with respect to $\mathbf q$ but $|V_i|^2$ is. The rest of the paper presents a local control algorithm to solve \eqref{eqn:main}.

\textbf{Remark:} In this paper we do not consider the constraints on bus reactive power injections $\mathbf q$. In an earlier work, where limited communication is present, all constraints were considered \cite{Lam12}. The next step of this work is to incorporate the constraints on reactive power injections. \hfill $\Box$

\subsection{Uniqueness of Power Flows} \label{sec:2A}
Given a power system with $n$ buses, the state of the system is described by $2n-1$ real numbers: the $n$ voltage magnitude and $n-1$ angles (taking bus $0$ to be the reference bus). This suggests that there are $2n-1$ degree of freedom in the power system. However, for transmission networks, specifying the voltage magnitudes and active powers $P_1,\dots,P_{n-1}$ is not enough to determine the state of the system since there could be multiple angles that result in those active powers \cite{BeVi:00}. In contrast, due to the tree structure of a distribution network \cite{LTZ12}, specifying voltage magnitudes and active powers at $n-1$ buses is enough to fully specify the state of the system. Therefore we often give the power system state in terms of voltage magnitude and active power; the corresponding angles and reactive powers are uniquely specified.

\section{Local Control Algorithm} \label{sec:result}
Consider a distribution system that is operating at its reference voltage, i.e., with $|V_i|=V_i^{ref}$ for $i=0 \dots n-1$. Denote the corresponding active and reactive injection vectors by $\overline{ \mathbf p}$ and $\overline{\mathbf q}$, respectively. Due to changes in renewables and loads, the active power injections in the system may change from $\overline{\mathbf p}$ to $\mathbf p'$. If the reactive powers stayed constant at $\overline{\mathbf q}$, the voltage magnitudes are no longer at the reference values. This section proposes an algorithm to bring the voltage magnitudes back to the reference values. By adjusting the reactive power, the algorithm relies  purely on local voltage measurements and does not involve communication among nodes in the network.

To correct the voltage magnitudes, there are $n-1$ control variables: the reactive power injections at all the buses except the feeder (bus $0$). We assume that the problem in \eqref{eqn:feasible} has a solution, or equivalently the optimal value of the problem in \eqref{eqn:main} is 0, and we design a local control algorithm that finds the optimal reactive power injection vector $\mathbf q$. Since bus $i$ can only measure local quantities, we propose the following iterative update algorithm:

\textbf{Local Control Algorithm:}
\begin{itemize}
\item Initialize the system at $t=0$.
\item At time $t$, bus $i$ measures $|V_i[t]|$ and makes the following update:
\begin{equation} \label{eqn:update_i}
Q_i[t+1]=Q_i[t]-d_i \left(|V_i[t]|^2-(V_i^{ref})^2 \right)
\end{equation}
for $i=1,\dots,n-1$. Intuitively, a bus should inject more reactive power to raise its voltage, so $d_i$ is assume to be positive (see \cite{Robbins12} for a deeper discussion).
\item $t \leftarrow t+1$ and repeat
\end{itemize}
By a slight abuse of notation, let $\mathbf q=[Q_1 \; \dots Q_{n-1}]^T$ since we do not explicitly control the reactive injection at the feeder. The matrix version of the update in \eqref{eqn:update_i} can be written as
\begin{equation} \label{eqn:update}
\mathbf q[t+1]= \mathbf q[t]- \mathbf D \left(|\mathbf v[t]|^2 - (\mathbf v^{ref})^2 \right),
\end{equation}
where $\mathbf D=\diag(d_1,\dots,d_{n-1})$ is a diagonal matrix and the squares are taken component-wise. Note both $|\mathbf v[t]|$ and $\mathbf v^{ref}$ are both vectors. We envision the local control algorithm would run in a continuous fashion to supplement the actions of devices such as tap-changing under-load transformers (TCULs) and switched capacitors. This local control algorithm is similar to the first-stage algorithm in \cite{Robbins12} where we replaced $|V_i|$ by $|V_i|^2$.

\section{Convergence of the Local Control Algorithm}
The convergence properties of the local control algorithm are investigated in this section. Since the active bus injections are the parameters in the voltage regulation problem, the convergence of the algorithm also depends on the prescribed active power injections. Recall from Section \ref{sec:model} that $\mathcal P$ is the feasible injection region determined by the active bus power constraints. In this section, we answer the following two questions:
\begin{enumerate}
\item Under what condition can we ensure that the algorithm converges for every $\mathbf p \in \mathcal P$?
\item How does network topology influence the convergence of the algorithm? That is, for the two networks in Fig. \ref{fig:star_line} and same active injection vector, does the algorithm converge for one of the networks but not the other one?
\begin{figure}[ht]
\centering
\subfigure[Star Topology]{
\includegraphics[scale=0.7]{star.eps}}
\subfigure[Line Topology]{
\includegraphics[scale=0.7]{line_s.eps}}
\caption{We are interested in how the topology of the network affect the convergence property of the local control algorithm. We show that it is easier for the local algorithm to converge on the star network than the tree network. Therefore the depth of the network limits the convergence of the algorithm.}
\label{fig:star_line}
\vspace{-0cm}
\end{figure}
\end{enumerate}

Since power system planners generally plan for the worst-case scenario, answering the first question would ensure that the algorithm would work for every possible $\mathbf p$ and communication is not needed to regulate the voltage. Theorem \ref{thm:sr} below states that it is in fact sufficient to check the convergence of the algorithm at a single $\mathbf p$ to ensure the converge for the entire feasible injection region $\mathcal{P}$.

In section IV-C, we show that the local control algorithm converges for the star topology under a much broader set of conditions than for the line topology. Since typical distribution networks are more like lines than stars, they are in some sense inherently more difficult to control only using local actions. Before stating the convergence theorems for the entire $\mathcal P$, we first state the result for a single power injection vector.
\subsection{Convergence for a Power Injection Vector }
The local control algorithm shows that we can view the $n-1$ voltages on the non-feeder buses as a function of the $n-1$ reactive power injections at these buses. The active powers and the feeder voltage can be interpreted as parameters in this function. From this point of view, we make the following definition:
Let $\mathbf M$ be the $n-1$ by $n-1$ Jacobian matrix where the $(i,k)$-th entry is
\begin{equation} \label{eqn:M}
M_{ik}= \frac{\partial |V_i|^2}{\partial Q_k};
\end{equation}
the partial derivative is evaluated at some active power and voltage magnitudes.

Note that $\mathbf M$ depends on both values of both the active powers and voltage magnitudes; however,  the variation of $\mathbf M$ with respect to voltage magnitude is relatively small over a wide range of operating conditions \cite{BeVi:00}. In this regard, we will always evaluate $\mathbf M$ at some $\mathbf p$ and the reference voltage magnitudes. To emphasize the dependence of $\mathbf M$ on $\mathbf p$, we write it as $\mathbf M (\mathbf p)$.

The next lemma gives the sufficient condition when the local control algorithm converges.
\begin{lem} \label{lem:1}
Suppose the optimization problem in \eqref{eqn:main} has optimal value $0$ (that is, the problem in \eqref{eqn:feasible} is feasible). Let $\mathbf p$ be the active power injection vector. The update in \eqref{eqn:update_i} converges if $\mathbf D \mathbf M(\mathbf p) + \mathbf M(\mathbf p)^T \mathbf D \succ 0 $, where $ \succ 0 $ denotes positive definiteness.
\end{lem}
The proof is an application of Prop. 2.1 in \cite{Bertsekas97} and is given in Appendix \ref{app:1}.

\subsection{Stability Region} \label{sec:sr}
The matrix inequality $\mathbf D \mathbf M(\mathbf p) + \mathbf M(\mathbf p)^T \mathbf D \succ 0 $ is in the form of a Lyapunov equation, and this motivates the next definition
\begin{defn}
An active power injection vector is called stable if there exists a diagonal matrix $\mathbf D$ such that $\mathbf D \mathbf M(\mathbf p) + \mathbf M(\mathbf p)^T \mathbf D \succ 0$. The set of all stable active power injections is called the stability region and is denoted by $\mathcal S_p$.
\end{defn}

Deciding the existence of a stabilizing $\mathbf D$ for a given $\mathbf p$ is a convex problem, and can be solved using semidefinite programming techniques \cite{Boyd04}. Theorem \ref{thm:sr} below shows how the stability for the entire region $\mathcal P$ can be checked.

Given the feasible injection region $\mathcal {P}$, let $\mathbf p_{\min}$ be the minimum injection vector in $\mathcal P$. That is, $\mathbf p_{\min}=[\underline{P}_1 \; \dots \; \underline{P}_{n-1}]$ (recall $\underline{P_i}<0$ is the lower bound on the active power injection, which correspond to the maximum demand). The next theorem gives a sufficient condition for the stability of $\mathcal P$ in terms of the stability of $\mathbf p_{\min}$.
\begin{thm} \label{thm:sr}
Let $\mathbf p_{\min}$ be the minimum injection vector.  Let $\bm \theta_{\min}$ be the corresponding angles when the voltage magnitudes are at their reference values (recall Section \ref{sec:2A}). Suppose $|\theta_{\min,ik}| \leq \tan^{-1}(b_{ik}/g_{ik})$. Then if $\mathbf M(\mathbf p_{\min})$ is stable, $\mathbf M(\mathbf p)$ is stable for all $\mathbf p \in \mathcal P$. That is, $\mathcal P \subseteq \mathcal S_p$.
\end{thm}

The proof of Theorem \ref{thm:sr} is given in the Appendix. The significance of the theorem is that it suffices to check the stability of a single point, and this would ensure the global stability of the entire $\mathcal P$. Thus, in practice it is easy to check whether $\mathcal P \subseteq \mathcal S_p$ and how much demand need to be reduced to achieve stability under local controls.

The condition on $\theta_{\min,ik}$ is discussed in detail in \cite{LTZ12}, and is equivalent to the thermal constraints on the transmission lines. They are expected to hold in almost all practical situations.

\subsection{Effects of Topology on the Stability Region} \label{sec:long}
Theorem \ref{thm:sr} states that the stability of $\mathcal P$ is determined by the minimum injection vector. In this section, we show that the stability region $\mathcal S_p$ shrinks as the depth of the network increases. This shows that long networks are inherently more difficult to control, or equivalently, a line network is inherently more unstable than a star network. Due to space constraints, rather than stating a general theorem, we state a result pertaining to line networks.
\begin{lem} \label{lem:long}
Consider a $n$-bus homogeneous line network where every line has admittance $1-j$. Then the stability region of $\mathcal S_p$ approaches the single point $\mathbf 0$ as $n$ goes to infinity.
\end{lem}

This lemma states that as the depth of the network increases, communication becomes critical in maintaining the voltage stability of the network.  Practical networks often have a large depth, e.g., the 123-bus system has a depth of 23; therefore, some communication would likely be necessary for large distribution networks.

\section{Case Studies} \label{sec:sim}
We analyze the stability of the proposed local control algorithm for three of the IEEE test feeders: the 8-bus, 34-bus and the 123-bus network (see \cite{testfeeders} for system data). Since only one demand point is given in the data, we compute $\mathbf M$ at that point. By Theorem \ref{thm:sr}, all demand vectors less than the given demand is stable. Table \ref{tab:stable} shows whether the networks is stable under its current demand.
\begin{table}[ht]
\centering
\begin{tabular}{|c|c|c|c|}
Networks & 8-bus & 34-bus & 123-bus \\
\hline
Depth & 6 & 19 & 23 \\
Stable & yes & yes & no
\end{tabular}
\caption{Stability of the three network at the demand given in \cite{testfeeders}.}
\label{tab:stable}
\vspace{-0.5cm}
\end{table}

\begin{figure}[!t]
\centering
\includegraphics[scale=0.35]{34sim_a.eps}
\caption{The behavior of voltage magnitudes under the local control algorithm. Each curve is the normalized voltage, i.e., $|V_i[t]|/V_i^{ref}$. All of the curves converges to $1$, which means that all voltages are at their reference.}
\label{fig:34sim_a}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[scale=0.35]{34sim_b.eps}
\caption{The behavior of voltage magnitudes under the local control algorithm. Each curve is the normalized voltage, i.e., $|V_i[t]|/V_i^{ref}$. The demand in the system is 3.5 times the demand in Fig. \ref{fig:34sim_a}. Due the heavier loading, the network becomes unstable under the local control algorithm and experiences large voltage swings.}
\label{fig:34sim_b}
\vspace{-0.7cm}
\end{figure}

We focus on the 34-bus network and demonstrate the action of the local control algorithm. Figures \ref{fig:34sim_a} and \ref{fig:34sim_b} respectively show the voltage magnitudes in the network when the local action is able to stabilize the voltage and when it fails to do so. In Fig. \ref{fig:34sim_a} the demand is set to be the current demand in \cite{testfeeders}, and initially, we vary the voltage randomly by 10 \% and apply the local control algorithm in Section \ref{sec:result}. The voltage magnitudes converge to the reference voltages.

The local control algorithm fails to converge once the demand in the network is increased to 3 times its current value. In Fig. \ref{fig:34sim_b}, we increase the demand to $3.5$ times the current demand, and again vary the voltage randomly by 10 \% from its reference value; the figure shows the unstable behavior of the voltage magnitudes.

\section{Conclusion} \label{sec:con}
In this paper we investigated the utilization of reactive-power-capable DERs for voltage regulation in distribution networks without explicit communications among the buses. We proposed an algorithm where each bus controls its reactive power injections sensing only its own voltage magnitude. We provided sufficient conditions when the algorithm converges for the entire feasible active power injection region of a network, and necessary conditions that shows longer and heavier loaded networks are harder to control using localized algorithms. We validated our claims in case studies using the IEEE 34-bus test feeder.

\bibliographystyle{IEEEtran}