\section{Introduction}

In many situations occurring in applied mathematics (for example mathematical models or  numerical simulation), when working with an input-output system with uncertain (random)  inputs, it is crucial to understand the global influence of one of the inputs on the output. This problem is generally called global sensitivity analysis (or in short sensitivity analysis). We refer, for example to
\cite{tarantola2008} and  \cite{saltelli2000} for an overview on practical aspects of sensitivity analysis.  Sensitivity analysis aims to give some quantitative indicator allowing to score the global influence of each input variable of the model. A very popular index well tailored for the case of a real valued output  is the so-called Sobol index early proposed in \cite{sobol1993}.   It is based on second order moments of the distributions related to the input-output system.
Different strategies have been implemented for its statistical estimation (see for instance \cite{gamboa2013_s} or \cite{rugama2018}).  This kind of second order indices may also be considered in the more general  frame where the output is a  real vector or a real function (see for instance \cite{klein2014} or \cite{marrel2011}).
Nevertheless, the variance-based indices have some drawbacks. Indeed, they only study the impact on the variance of the output and so provide  only restricted summary of the output distribution (see \cite{daveiga2015}). Hence,
considering higher order indices involving the whole distributions and not only their second order properties may give more accurate information on the system.
In \cite{pianosi2015},  a sensitivity measure based on the Kolmogorov-Smirnov statistic is proposed and studied while entropy-based sensitivity measures are considered for instance
in \cite{liu2006r}  and \cite{borgonovo2007}.
For a general overview we also refer to  \cite{borgonovo2017}.
When the output $Z=F(X)$ ($X=(X_i)$ is the input vector), is valued in a more general space, the authors of  \cite{borgonovo2016} have proposed the following general sensitivity index for the input $X_i$:
\begin{displaymath}
S_i \mathrel{\mathop:}= \mbox{$\mathbb{E}$}_{X_i} \left[ d \left(P_Z, P_{Z \vert X_i} \right) \right].
\end{displaymath}
Here,  $d(\cdot, \cdot)$ is a given dissimilarity measure between  probability measures, $\mbox{$\mathbb{E}$}_{X_i}(\cdot)$ denotes the expectation with respect to $X_i$ and $P_Z$  (resp. $P_{Z \vert X_i}$) is the unconditional (resp. conditional) probability. Notice that in general, the study of an accurate statistical estimation of such index is not obvious.

Up to our knowledge, the case where the output takes its values on a Riemannian manifold has not yet been studied.
Beside, in \cite{gamboa2015} a  sensitivity measure in the case of a scalar output and based on the Cram\'er von Mises distance is proposed an studied. This paper extends this last approach to the more general frame of an output valued on a Riemannian manifold.

In the last decades starting from the pioneer 1945' work of Rao \cite{rao1945}, the statistical theory for data valued on a Riemannian manifold have received a lot of interest and contributions. One of the reason of the development of such theory is the spectacular increase of the computation power allowing the treatment of more and more complex objects and structures on computers. References on the subject are numerous. We refer to \cite{mardia1972}, \cite{bhattacharya2012}  and  \cite{ellingson2015} and the references therein for an overview. In this paper, our aim is to bridge this theory to sensitivity analysis. We build a general sensitivity index. On the one hand this index takes into account the whole characteristic of the involved distributions and not only their second order moments.  On the other hand the working frame is a general system with an output valued on a Riemannian manifold. One of the main ingredient of our work is a generalization of the Cram\'er von Mises criterion replacing, before integrating, the half lines by geodesic balls. Consequently, our new index has the nice property to be invariant with respect to the isometric maps of the Riemannian manifold. Moreover, we show that this generalized index can be easily estimated using a Pick Freeze like estimator. This estimator involves a $U$-statistics. Without loss of generality,  we will assume that the random input variables are real and independent. Sensitivity analysis for dependent inputs may be performed but is generally less readable (see for example  \cite{gamboa2015b} and references therein).

 The paper is organized as follows. In Section  \ref{conceptosbasicos}, in order to be self contained, we recall some key facts and properties on Riemannian manifolds. In Section \ref{indicedef}, we define our new sensitivity index for  Riemannian manifold. Further, we will build a Pick an Freeze estimator based on a $U$-statistic and study its asymptotic (consistency), and non asymptotic (concentration inequality) properties.    In Section \ref{datosreales}, we study the sensitivity of a real life model coming from mechanic. The output of the system is the rigidness matrix. This matrix is involved to model the linear elasticity of a solid object.  The set of these matrices is as a sub-manifold of the symmetric positive matrices. In the last section, we bold the main advantages of our method and discuss its possible extensions. The codes in  Julia and R languages (see \cite{julia2012} and \cite{team2017} respectively) are available upon request to the authors. All proofs are postponed to the Appendix.

\section{Basic concepts}

\label{conceptosbasicos}
Let us begin with  some useful tools and facts in Riemannian geometry.  A \textit{Riemannian metric} $g$ on a manifold $\mbox{$\mathcal{M}$}$ allows to define for every point $p \in \mbox{$\mathcal{M}$}$  a scalar product $g_p(\cdot,\cdot)$ acting on the tangent space (at $p$), $T_p \mbox{$\mathcal{M}$}$. This scalar product depends smoothly  on $p$.   The \textit{Riemannian manifold} $(\mbox{$\mathcal{M}$},g)$ is the manifold equipped with the Riemannian metric  $g$.
For any $v \in T_p \mbox{$\mathcal{M}$}$ the Riemannian norm of $v$ is given by $\Vert v \Vert \mathrel{\mathop:}= \sqrt{g_p(v,v)}$.
Let $x,y\in M$ and  $\gamma: I:=[0, 1] \rightarrow \mbox{$\mathcal{M}$}$ be a continuously differentiable curve contained in the Riemannian  manifold $\mbox{$\mathcal{M}$}$ (that is assumed to be connected) with
$\gamma(0)=x$ and $\gamma(1)=y$. Let denote by  $L(\gamma)$ its length. We may now define the \textit{induced distance}  $d_g(x, y)$ between the points $x$ and $y$ of $\mbox{$\mathcal{M}$}$  by setting
    $$d_g(x,y) = \inf \{ L(\gamma) : \gamma \,
 \}.$$
$d_g$ is the Riemannian distance on $\mbox{$\mathcal{M}$}$  with respect to the $g$ metric. A \textit{geodesic} (with speed $s \in \mathbb{R}^{+}_0$) is a smooth map $\alpha:I \rightarrow \mbox{$\mathcal{M}$}$, such that $\Vert \alpha^{'}(t) \Vert= s$  for all $t \in I$ and which is {\it locally length minimizing}. For $p \in \mbox{$\mathcal{M}$}$ and  $v\in T_p \mbox{$\mathcal{M}$}$, there exists a unique geodesic  $\alpha_{(p,v)}(t)$ starting from that point with initial tangent vector $v$. The \textit{exponential map} is the map $\exp_p$ given by $\exp_p (v) := \alpha_{(p,v)}(1)$.

Notice further that the geodesic that joins two points is not necessarily unique.
The \textit{cut locus}  of $p$ in the tangent space is defined to be the set of all vectors $v \in T_{p} \mbox{$\mathcal{M}$}$ such that  $\exp _{p}(tv)$ is a minimizing geodesic for  $0 \leq t \leq 1$  but fails to be one for $0 \leq t \leq 1+ \epsilon$, $\epsilon >0$. The \textit{cut locus of $p$ in $\mbox{$\mathcal{M}$}$}, denoted  by $C_{\mbox{$\mathcal{M}$}}(p)$, is defined as the image of the cut locus of  $p$ in the tangent space under the exponential map at $p$. The \textit{injectivity radius of $p$} is the maximal radius of centered balls on which the exponential map is a diffeomorphism. The \textit{injectivity radius of the manifold}  $r_{iny}$ is the infimum of the injectivity over the manifold. For example, in the sphere $S_{d}$ of $\mbox{$\mathbb{R}$}^{d+1}$ the cut locus of a point $p$ is its antipodal point  $-p$. There are infinite minimizing geodesics that connect a point with its antipodal (Figure  \ref{bola_diametro} shows two geodesics) but this configuration has zero probability of being drawn when dealing with the uniform probability measure.

\begin{figure}[!ht]
\centering
\subfigure{\includegraphics[width=50mm]{bola_diametro3.png}}
\subfigure{\includegraphics[width=40mm]{asas2.pdf}}
\caption{\textit{Left Panel}: Representation of a geodesic ball of diameter $\overline{PQ}$ and center $O$ on a hyperbolic paraboloid. \textit{Right Panel}: A graph of two minimizing geodesics that join two antipodal points of a sphere.}
\label{bola_diametro}
\end{figure}

Let $(\mbox{$\mathcal{M}$},g)$ be a connected and orientable Riemannian manifold,(see \cite{docarmo1992}, page $18$). We will assume that  $(\mbox{$\mathcal{M}$},d_g)$ is a complete  separable metric space. Since $(\mbox{$\mathcal{M}$},d_g)$ is complete, the Hopf--Rinow theorem (see \cite{docarmo1992}, p.146) implies that for any pair of points $p,q \in \mbox{$\mathcal{M}$}$  there exist at least one geodesic path in $\mbox{$\mathcal{M}$}$ connecting $p$ and $q$.
If the manifold $\mbox{$\mathcal{M}$}$ fulfills the assumptions of the Cartan--Hadamard Theorem (see \cite{petersen2006}, p.162)--that is, if $\mbox{$\mathcal{M}$}$ is a simply connected, complete Riemannian manifold with non--positive curvature (a Hadamard manifold), then the geodesic is unique.
 Let $X$ be a random element taking values in $\mbox{$\mathcal{M}$}$, with distribution $P$. In many useful examples (eg. the sphere),  it occurs that geodesic uniqueness fails.  But in what follows, we will only need that this failure not occurs too often. Roughly speaking, we will assume a condition on $P$ ensuring that there is a unique geodesic between each pair of points $p,q \in \mbox{$\mathcal{M}$}$ with  probability one. More precisely, we will assume that the random element $Y$ has a density $f_Y$ with respect to the volume measure $d \nu(y)$ on $\mbox{$\mathcal{M}$}$, (which exists since $\mbox{$\mathcal{M}$}$ is orientable - see  the last section of \cite{folland2013}) fulfilling the following condition:
Given $q \in \mbox{$\mathcal{M}$}$ let
\begin{equation}
A_q\mathrel{\mathop:}= \left\{y \in \mbox{$\mathcal{M}$}  \Big/
\begin{array}{ll}
\textrm{ there are more than one different } \\
\textrm{minimizing geodesic connecting $y$ and $q$}
\end{array} \right \}.
\end{equation}
Then, there exist a Borel set $B_q \subset \mbox{$\mathcal{M}$}$ with $A_q \subset B_q$ and such that for any  $q \in \mbox{$\mathcal{M}$}$
\begin{equation}
\label{unicidadd}
\int_{B_q} f_Y(y) d \nu(y)=0.
\end{equation}
\begin{remark}
Let $C_{\mbox{$\mathcal{M}$}}(p)$ stand for the  cut locus of $p$ in a complete manifold $\mbox{$\mathcal{M}$}$, (see \cite{docarmo1992}, p. $267$).  If $P \left( Y \in C_{\mbox{$\mathcal{M}$}}(p)\right) = 0$,  for all   $p \in \mbox{$\mathcal{M}$}$, then condition  (\ref{unicidadd}) is fulfilled (see \cite{pennec2006}). For instance, if the Riemannian manifold is Hadamard, then the cut locus of any point $p$ is empty and the condition is obviously fulfilled. In the unit sphere $S_{d} \mathrel{\mathop:}= \{u \in \mathbb{R}^d / \Vert u \Vert= 1 \}$  the cut locus of a point $p$ is its opposite point  $-p$ and any  probability measure having a density supported by the sphere fulfills condition  (\ref{unicidadd}).
\end{remark}
In this paper,  we will assume that the Riemannian manifold $(\mbox{$\mathcal{M}$},g)$ with the induced distance $d_g$ is connected and oriented, and that the metric space $(\mbox{$\mathcal{M}$},d_g)$ is separable and complete. We will also assume  that given two points $p,q \in \mbox{$\mathcal{M}$}$  there is a unique geodesic determined by $p,q$ with probability one with respect to the tensorial probability measure $d\xi(p,q) \mathrel{\mathop:}= f_Y(p)f_Y(q)d \nu(p)d\nu(q)$, (see Figure \ref{geodesica}).

\begin{figure}[!ht]
\centering
\subfigure{\includegraphics[width=38mm]{cara.png}}
\subfigure{\includegraphics[width=50mm]{geodesica_dona.png}}
\caption{\textit{Left Panel}: A graph of a geodesic that determines two points on the representation
of a face. \textit{Right Panel}: A graph of a geodesic that determines two points on the representation
of a Torus.}
\label{geodesica}
\end{figure}

\subsection{Geodesic balls with diameter $\overline{pq}$}
\label{bolasgeodesicas}
For any pair $p,q \in \mbox{$\mathcal{M}$}$ that determine a unique geodesic   $\overline{pq}$,  we define the ball of diameter  $\overline{pq}$  as the closed ball whose center is the middle point of the geodesic joining $p$ and $q$ with radius $d_g(p,q)/2$. It will be denoted by $B_{pq}$. Figure  \ref{bola_diametro}  depicts  a geodesic ball of diameter $p$ and $q$, where $o$ is the midpoint of the geodesic.
In \cite{christensen1970} it is shown that the family  $\mathcal B_p:= \{{B}_{p}, p \in \mbox{$\mathcal{M}$}\}$ of geodesic balls with center $p$ is a determining class; that is, if two probability measures $\eta, \nu$ on $\mbox{$\mathcal{M}$}$  coincide on any event of this class, then $\eta=\nu$.
We assume that $\mbox{$\mathcal{M}$}$ is a compact Riemannian manifold. We will start by proving that the family of balls $$\mathcal{B}_{pq}:= \{B_{pq}, p,q \in \mbox{$\mathcal{M}$}\ \textrm{and there is a unique geodesic between $p$ and $q$} \}$$ is also a determining class if the radius of injectivity  $r_{iny}$ of $\mbox{$\mathcal{M}$}$ is positive.  We will need the following weaker assumption:

\textbf{HB) B-continuity} A probability measure $\nu$ defined on a Riemannian manifold $\mbox{$\mathcal{M}$}$ fulfils \bf HB \rm if $\nu(\partial A )=0$  for all closed Borel sets  $A$ on $\mbox{$\mathcal{M}$}$ (Here $\partial A$ denotes the topological boundary of the set $A$).

\begin{property}
\label{propiedad1}
Let $(\mbox{$\mathcal{M}$},g)$ be a compact Riemannian manifold and $\nu$ be a probability distribution fulfilling \textbf{HB}, such that $\nu \left( \partial \mbox{$\mathcal{M}$} \right)=0$. If the injectivity radius is positive  then  $\mathcal{B}_{pq}$ is a determining class of $\nu$.
\end{property}
 In  \cite{fraiman2017} we showed what if $K \in \mbox{$\mathcal{M}$}$ is a compact set, then the family of balls
$\mathcal{B}_{pq}$ is also a Glivenko--Cantelli class in $K$; that is, $$\sup_{p,q \in K}  \vert  P \left(B_{pq} \right) - P_n\left(B_{pq} \right) \vert
 \rightarrow 0  \quad \textrm{a.s. \ \  as  $n \rightarrow + \infty$}.$$

In the next section, we define a sensitivity index using this family of sets.

\section{Sensitivity Index in Geodesic Balls}
\label{indicedef}
In order to build a sensitivity index, we use an idea previously developed in \cite{gamboa2015}. In our frame we  compare the distribution $P_Y$ and the conditional distribution $P_{Y \vert X_i}$ on the family of geodesic balls on the manifold (instead of half--lines in \cite{gamboa2015}. So that, we build an intrinsic index related to the manifold. This index depends only on the dimension and structure of the manifold and not on the dimension of the space where it is immersed, unlike in \cite{klein2014} and \cite{gamboa2015}.

\subsection{Constructing the index}

Let $\mathbb{X}=(X_1, \ldots, X_d) \in \mbox{$\mathbb{R}$}^d$ be an random vector and assume that $\mathbb{P} \mathrel{\mathop:}=  \mathbb{P}_1 \times \ldots \times \mathbb{P}_d$ is the probability law of $\mathbb{X}$. Further, let
\begin{displaymath}
Z=f(X_1, \ldots, X_d),
\end{displaymath}
\noindent where $f: \mathbb{R}^d \rightarrow  \mbox{$\mathcal{M}$}$ is a continuous function. Here, $\mbox{$\mathcal{M}$}$ a Riemannian manifold of dimension $k$. We wish to understand how sensitive is the  output $Z$ to perturbations in some of the input variables $(X_1,\ldots, X_d)$. As any manifold can be immersed in $\mathbb{R}^p$ for a sufficiently large $p$ (see \cite{nash1954}), this problem could be considered as a particular case of a multivariate output. Nevertheless,  in such a way the geometry and \textit{``minor dimensionality''} of the output are not taken into account. We begin by studying a very particular case and then extend our results to any Riemannian manifold, that satisfies   the previous conditions. Let $h$  a measurable function such that $h(X_1,X_2)$ is integrable, we set $$\mbox{$\mathbb{E}$}_{X_2}\left(h(X_1,X_2) \right) \mathrel{\mathop:}= \mbox{$\mathbb{E}$} \left( h(X_1,X_2)  / X_1 \right).$$

\subsection{A very particular case: The real line.}

Let us first focus on the very particular case  $\mbox{$\mathcal{M}$}=\mbox{$\mathbb{R}$}$. Let $F$ be the distribution function of $Z$, $F(t) \mathrel{\mathop:}= P(Z \leq t) \, (t\in \mathbb{R})$, and $F^{\nu}$ be  the distribution of $Z$ conditioned by a subset of the variable $(X_1, \ldots, X_d)$. That is, let $\nu= \{i_1,\ldots, i_k\}  \subset \{1, \dots, d\}$, and

\begin{align*}
& \hspace{-0.5cm} X_{\nu}\mathrel{\mathop:}= (X_{i_1}, \ldots, X_{i_k}),\\
&\hspace{-0.5cm} F^{\nu}(t)\mathrel{\mathop:}=   \mbox{$\mathbb{E}$} \left( \mathbbm{1}_{(\infty,t]} (Z) \big \vert X_{\nu} \right).
\end{align*}

\noindent In \cite{gamboa2015} the following Cram\'er von Mises sensitivity index is considered. Assume that $F$ is absolutely continuous. This normalized index (denoted  by  $C^{\nu}_{2}$) is defined as

\begin{equation}
C^{\nu}_{2} \mathrel{\mathop:}= \frac{N^{\nu}_{2,}}{\int_{\mathbb R} F(x)\left(1-F(x) \right) dF(x)} = 6 N^{\nu}_{2},
\end{equation}

where
\begin{align*}
N^{\nu}_{2} =\int  \mbox{$\mathbb{E}$} \left( \left[  F(t) -F^{\nu}(t) \right ]^2 \right) dF(t) \\
&\hspace{-4.7cm}= \int \left ( \mbox{$\mathbb{E}$}  \left[   \int \mathbbm{1}_{(\infty,t]}(z)  d(F -F^{\nu})(z) \right ]^2\right )  dF(t)\\
&\hspace{-4.7cm}= \mbox{$\mathbb{E}$} \left\{ \mbox{$\mathbb{E}$}  \left(  \left[ \mbox{$\mathbb{E}$} \left(\mathbbm{1}_{(\infty,t]}(Z) \right) - \mbox{$\mathbb{E}$} \left(\mathbbm{1}_{(\infty,t]}(Z) \big \vert X_{\nu} \right)  \right]^2 \right) \right\}.
\end{align*}

\noindent To build  of our new index, we replace the  function $\mathbbm{1}_{(\infty,t]}$  by the indicator function of the interval $[\min(s,t), \max(s,t)]$ (denoted by $h_{s,t}$),

\begin{displaymath}
h_{s,t}(x) \mathrel{\mathop:}= \mathbbm{1}_{ \{s \leq x \leq t\} }   +   \mathbbm{1}_{\{t \leq x \leq s\}}=   \mathbbm{1}_{ \{\min\{s,t\} \leq x \leq \max\{s,t\}\}}.
\end{displaymath}

\noindent The function $h_{s,t}$, can be thought of as the indicator function of the ball of diameter $\overline{st}$ in $\mathbb{R}$. Further, let

\begin{equation}
\label{bla}
H(s,t)\mathrel{\mathop:}= \mbox{$\mathbb{E}$}\left[ h_{s,t} (Z)\right] \quad \textrm{and}  \quad  H^{\nu}(s,t)= \mbox{$\mathbb{E}$}\left[ h_{s,t} (Z) \Big \vert X_{\nu} \right].
\end{equation}
\noindent  Obviously, $\mbox{$\mathbb{E}$}_{X_{\nu}}\left[ H^{\nu}(s,t) \right]= H(s,t)$.

\begin{definition}
\label{sensibilidad}
The normalized ball sensitivity index (denoted by $B^{\nu}_{2}$) is then defined as,
  \begin{equation}
  B^{\nu}_{2} \mathrel{\mathop:}= \frac{S^{\nu}_{2}}{\int_{\mathbb R^2} H(x,y)\left(1-H(x,y) \right) dF(y)dF(x)} = 6 S^{\nu}_{2},
  \end{equation}
  where
  \begin{equation}
  S^{\nu}_{2} \mathrel{\mathop:}= \mbox{$\mathbb{E}$}_{Z_1,Z_2} \left[ \mbox{$\mathbb{E}$}_{X_{\nu}} \left \{  \left[ H(Z_1,Z_2) - H^{\nu}(Z_1,Z_2)\right]^2 \right \} \right],
  \end{equation}
\noindent  and $Z_1$ and $Z_2$ are two independent copies of $Z$.
(The constant 6 follows from the well known fact that $\int_{t>s} (t-s)(1-(t-s))dtds = \int_{t\leq s} (t-s)(1-(t-s))dtds=1/12$.)

  \end{definition}

We may can rewrite $S^{\nu}_{2}$ as

\begin{displaymath}
S^{\nu}_{2}=\int_{ \mbox{$\mathcal{M}$} \times \mbox{$\mathcal{M}$}}  \mbox{$\mathbb{E}$}_{X_{\nu}} \left \{  \left( H(z_1, z_{2})-  H^{\nu}(z_1,z_{2})\right)^2 \right \} dF(z_1) dF(z_{2}).
\end{displaymath}

\subsection{Generalization for a  Riemannian manifold}

As discussed before, $\mbox{$\mathcal{M}$}$ is a Riemannian manifold of dimension $k$. So that, under some regularity conditions, given two points  $\{z_1, z_2\} \subset \mbox{$\mathcal{M}$}$, the following  function $h_{z_1, z_2 }: \mbox{$\mathcal{M}$} \rightarrow \{0,1\}$ is generically well defined,

\begin{displaymath}
h_{ z_1, z_2  }(t) \mathrel{\mathop:}= \mathbbm{1}_{ B_{z_1z_2}} (t).
\end{displaymath}

\noindent Let $Z_1,Z_2$ be independent copies of $Z$. Our normalized  sensitivity index is built by means of balls as in Definition \ref{sensibilidad},

 \begin{equation}
 \label{indic}
  B^{\nu}_{2} \mathrel{\mathop:}= \frac{S^{\nu}_{2}}{D^{\nu}_{2}},
  \end{equation}
where
\begin{displaymath}
S^{\nu}_{2} \mathrel{\mathop:}=  \mbox{$\mathbb{E}$}_{Z_1,Z_2} \left[ \mbox{Var}_{X_{\nu}} \left \{ \mbox{$\mathbb{E}$}_Z \left[ h_{Z_1,Z_2}(Z) \Big \vert X_{\nu} \right] \right\} \right],
\end{displaymath}

and

\begin{displaymath}
D^{\nu}_{2} \mathrel{\mathop:}= \mbox{$\mathbb{E}$} \left[ H(Z_1,Z_2) \left(1-H(Z_1,Z_2) \right) \right].
\end{displaymath}

\noindent  Notice that the index defined in the previous subsection is a particular case with  $\mbox{$\mathcal{M}$}=\mbox{$\mathbb{R}$}$.

\begin{remark}
 If   $B^{\nu}_{2}=0$ we have that  $\mbox{$\mathbb{E}$}_Z \left[ h_{z_1,z_2}(Z) \Big \vert X_{\nu} \right]= \mbox{$\mathbb{E}$}_Z \left[ h_{z_1,z_2}(Z) \right] \, a.s.$ for all $(z_1,z_2) \in \Omega \subset \mathcal{M}\times\mathcal{M}$ with $P(\Omega)=1$. Therefore, under the assumptions of Property \ref{propiedad1},  the probability measures  $P_{Z \vert X_{\nu}}$ and $P_{Z}$ are the same.
\end{remark}

\subsection{Estimation}

In the next subsection, an estimator of the index $B^{\nu}_{2}$ is proposed.   The variance of the conditional mean is estimated using the pick and freeze method (see \cite{sobol1993} and \cite{sobol2001}). Further, as the expected value  $\mbox{$\mathbb{E}$}_{Z_1,Z_2}(\cdot)$  is a symmetric function of $Z_1$ and $Z_2$ we will use an $U$--statistic. At the end  of the next subsection, the steps for the construction of the estimator will be detailed.

\subsubsection{Estimation by \textit{``Pick and Freeze''} method}
The method consists in writing the variance of the conditional mean as a covariance. For $\nu \in \{1, \ldots, d \}$, let $X^{\nu}$ be the random vector that coincides with $X$ on its $\nu$ components and is independently regenerated on the others components. That is, $X^{\nu}_{\nu}=  X_{\nu}$  and $X^{\nu}_{\bar{\nu}}=  X^{'}_{\bar{\nu}}$, where $X^{'}$ is an independent copy of $X$ and $\bar{\nu}$  is the complementary set of $\nu$ ($\bar{\nu}= \{1, \ldots, d \} \setminus \nu $). We set
\begin{displaymath}
Z^{\nu} \mathrel{\mathop:}= f(X^{\nu}).
\end{displaymath}

\noindent For sake of simplicity assume first that  $\mbox{$\mathbb{E}$}(Z)=0$. Then,

\begin{equation}
\mbox{Var}\left(  \mbox{$\mathbb{E}$}(Z \big \vert X_{\nu}) \right)= \mbox{$\mathbb{E}$} \left(  \mbox{$\mathbb{E}$}^2(Z \big \vert X_{\nu}) \right).
\label{pick}
\end{equation}

\noindent Since $Z$ and $Z^{\nu}$ are conditionally independent of $X_{\nu}$ (see \cite{klein2013}), so that

 \begin{align*}
 \mbox{cov} (Z, Z^{\nu})=  \mbox{$\mathbb{E}$}(Z Z^{\nu}) - \mbox{$\mathbb{E}$}(Z)\mbox{$\mathbb{E}$}( Z^{\nu})=  \mbox{$\mathbb{E}$}(Z Z^{\nu}) - \mbox{$\mathbb{E}$}^2(Z)\\
 &\hspace{-6.95cm} = \mbox{$\mathbb{E}$} \left[ \mbox{$\mathbb{E}$}(Z Z^{\nu} \big \vert X_{\nu}) \right] = \mbox{$\mathbb{E}$} \left[ \mbox{$\mathbb{E}$}( Z \big \vert X_{\nu})  \mbox{$\mathbb{E}$}( Z^{\nu} \big \vert X_{\nu}) \right] \\
  &\hspace{-6.95cm} =  \mbox{$\mathbb{E}$} \left[  \mbox{$\mathbb{E}$}^2  (Z \vert X_{\nu}) \right].
 \end{align*}

\noindent  Following  \cite{klein2013}, estimating the covariance  by a Monte Carlo Method we obtain the following  Pick Freeze estimator of  $\mbox{Var} \left[  \mbox{$\mathbb{E}$}  (Z \vert X_{\nu}) \right]$,

\begin{equation}
T^{\nu} \mathrel{\mathop:}= \frac{1}{N}\sum_{j=1}^N Z_j Z_j^{\nu} - \left( \frac{1}{2N} \sum_{j=1}^N (Z_j + Z_j^{\nu}) \right)^2.
\end{equation}

\noindent Here, $Z_j$ and $Z^{\nu}_j$,  are $N$ independent copies of $Z$  and  $Z^{\nu}$ respectively.

\vspace{0.2 cm}
\underline{\emph{First step}.} Using the previous idea, we build  a consistent estimator of the numerator $S^{\nu}_{2}$.
Let $\mathcal{P}_{N,p}$ be the set of all possible subsets of $\{1, \ldots,N\}$ with $p$ elements, that is,

\begin{displaymath}
\mathcal{P}_{N,p}= \{(i_1, \ldots, i_p) \in \{1, \ldots, N \}^p / i_1< \ldots < i_p  \}.
\end{displaymath}

\noindent Let $\tau=(i_1, \ldots, i_p) \in \mathcal{P}_{N,p}$. We denote by $\textbf{W}_{\tau}= \left( W_{i_1},\ldots ,W_{i_{p}} \right)$ a sample of $p$ independent copy of $Z$.
The estimator is obtained as follows,

\noindent  Further,  $D^{\nu}_{2}$  is  estimated analogously setting

\begin{align*}
\hat{D}^{\nu}_{2} \mathrel{\mathop:}= \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2}} \Big \{ \frac{1}{2N} \sum_{j=1}^N   \left( h_{\textbf{W}_{\tau}} \left( Z_j\right)+ h_{\textbf{W}_{\tau}} \left( Z_j^{\nu}\right) \right)  -  \\
& \hspace{-7.5cm}  - \left( \frac{1}{2N} \sum_{i=1}^N \left[ h_{\textbf{W}_{\tau}} \left( Z_i\right)+ h_{\textbf{W}_{\tau}} \left( Z_i^{\nu}\right) \right] \right)^2   \Big \}.
\end{align*}

\noindent The computation  of the  estimator is simple. Indeed, it only involves computation of indicator functions. Nevertheless, the computational time to perform  this estimator could be high. Indeed,  the number of sums is of  order $N^3$. Each term involves the determination of a geodesic. In the following section, we give some asymptotic properties of this estimator. We provide an exponential inequality that leads, from the Borel--Cantelli Lemma, to strong consistency.

\subsection{Asymptotic properties of  $\hat{B}^{\nu}_{2}$.}
We analyzed separately the strong convergence of $\hat{S}^{\nu}_{2}$ and $\hat{D}^{\nu}_{2}$. The strong convergence of  $\hat{S}^{\nu}_{2}$  follows from the next Theorem. Notice that another proof may be built using the fourth moment and the Rosenthal inequality for $U$--statistic of order $2$ developed in \cite{ibragimov1999}.

\textbf{Theorem (Exponential inequality )}
\label{exponencial}
Let $s>0$, there exist $N_0$ such that if $N>N_0$,
\begin{equation}
P \left( \left \vert \hat{S}^{\nu}_{2} -  S^{\nu}_{2} \right \vert > s   \right) \leq 16 \exp \left\{  -\frac{N \left(\frac{s}{9}\right)^2}{8} \right\}.
 \end{equation}

From the previous Theorem and the Borel--Cantelli Lemma the strong consistency of $\hat{S}^{\nu}_{2}$ follows.

\textbf{Corollary (Consistency of the estimator)}
\label{consistencia}
$\hat{S}^{\nu}_{2}$ is a consistent estimator of $S^{\nu}_{2}$.

Notice that we may analogously show the convergence of the denominator involved in (\ref{indic}).
\section{Simulations}
\label{simulaciones}
This simulation   section is made up on three examples. The first one shows the accuracy of the estimator  when the output is real valued. In the second example the output lies on a circle immersed in $\mathbb{R}^2$. The estimates are compared with those obtained in \cite{gamboa2015}. The last example shows that the index proposed in \cite{gamboa2015} may fail as sensitivity indicator.

\subsection{Example 1: Output on the real line}
\label{ejemplo1}
We study here an example where the Sobol index does not provide information on sensitivity, while $D^{\nu}_{2}$  and $B^{\nu}_{2}$ do. Let $Z=\alpha X_1 + X_2$  where $\alpha>0$, and $X_1,X_2$ are independent. Assume further that $X_1 \sim \textrm{Bernoulli}(p)$, $X_2 \sim F$, with $m=\mbox{$\mathbb{E}$}(X_2)$, $\sigma^2=\textrm{Var}(X_2)=\alpha^2 p(1-p)$.

We calculate $B^{1}_{2}$ in the Appendix while  $C^{1}_{2}$ has been calculated in \cite{gamboa2015}. In the case when $X_2 \sim U(0,b)$ with $b= \sqrt{12\alpha^2 p(1-p)}$. We find
\begin{displaymath}
B^{1}_{2}= 12p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^3  \left(\frac{1}{3}- \frac{1}{4}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/12 & \text{if } \alpha > b,
\end{cases}
\end{displaymath}
\begin{displaymath}
C^{1}_{2}= 6p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^2  \left(1- \frac{2}{3}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/3 & \text{if } \alpha > b,
\end{cases}
\end{displaymath}

In the Figure \ref{index}, the estimates of both indices are compared for different values of $p \in [0,1]$. Samples sizes $100$, $500$ and $1000$ are considered.  For each $p$, through resampling Bootstrap for $U$--statistics (see \cite{gine1992b}), $95\%$ confidence intervals are also given. In all cases, the mean squared deviation (MSD) of the estimators are computed.
The MSD of $\hat{B}^{1}_{2}$ is slightly smaller than the one for $\hat{C}^{1}_{2}$ (see Table \ref{MSDrecta}). Similar results are obtained when comparing  $\hat{B}^{2}_{2}$  and  $\hat{C}^{2}_{2}$.

\begin{table}[ht]
\centering
\begin{tabular}[t]{ccc}
\toprule
 Size &   $MSD_{\hat{B}^{1}_{2}}$ &  $MSD_{\hat{C}^{1}_{2}}$ \\
\midrule
$N=100$    & 0.051  & 0.067 \\
$N=500$    & 0.022 & 0.028 \\
$N=1000$   & 0.013  & 0.018  \\
\bottomrule
\end{tabular}
\caption{The MSD of the  $\hat{B}^{1}_{2}$ and $\hat{C}^{1}_{2}$ estimators for sample sizes $N=100,500$ and $1000$.}
\label{MSDrecta}
\end{table}%

We may conclude that the new proposed index has a similar behavior as the Cram\'er von Mises one for  the real line.

\begin{figure}[ht]\centering
\includegraphics[width=.3\textwidth]{toyN100.pdf}
\includegraphics[width=.3\textwidth]{toyN500.pdf}
\includegraphics[width=.3\textwidth]{toyN1000.pdf}
\caption{Cram\'er von Mises indexes and our new indexes for Example \ref{ejemplo1}. In all cases,   ($\textendash \, \, \textendash$) are the true index  and ($\edge$) their estimate. Violet color depicts CVM sensitivity index while orange color depicts Ball sensitivity index. The $95\%$ bootstrap confidence intervals are represented with shading.  Left--hand Panel: $N=100$. Center--hand Panel: $N=500$. Right--hand Panel: $N=1000$.}
\label{index}
\end{figure}

\subsection{Example 2: Output on a simple manifold immersed in $\mathbb{R}^2$}
\label{ejemplo2}
The case where $Z$ ranges on the unit circle $S_1$ of  $\mathbb{R}^2$ is considered. We assume that the input vector has distribution  $$X=\binom{X_1}{X_2}\sim N \left[ \binom{\mu_1}{\mu_2},  \left(
\begin{array}{cc}
\sigma^2_1 & 0 \\
0 & \sigma^2_2
\end{array} \right)  \right].$$

Here, we study the case when the output $Z$ is the normalized version of $X$

\begin{displaymath}
\label{circu}
Z \mathrel{\mathop:}= \frac{X}{\Vert X \Vert}.
\end{displaymath}

The distribution of $Z$ has been widely studied, we refer for example to \cite{mardia1972}, \cite{kendall1974} and \cite{watson1983}. In Figure \ref{circulo} the  index proposed  in \cite{gamboa2015} and the one studied here are depicted for the last system. These sensitivity indices are computed for $\mu_1 \in [-5,0]$. The other values are setted to  $\mu_2=0$ and $\sigma^2_1=\sigma^2_2=1$. We observe that the variability of $\hat{B}^{\nu}_{2}$ is  smaller that  the one of $\hat{D}^{\nu}_{2}$  for a sample size $N=300$.

\begin{figure}[ht]\centering
\includegraphics[width=.49\textwidth]{circuloB1.pdf}
\includegraphics[width=.49\textwidth]{circuloB2.pdf}
\caption{Calculation of the indices $\hat{B}^{\nu}_{2}$ and $D^{\nu}_{2}$ in the equation \ref{circu} with $\mu_1 \in [-5,0]$ and $N=300$. For all cases, it is represented using  ($\edge$) the index estimation functions for $N=300$. In violet color depicts CVM sensitivity index and orange color depicts Ball sensitivity index. The $95\%$ bootstrap confidence intervals are represented with shading. Left--hand Panel: $\nu=1$. Right--hand Panel: $\nu=2$.}
\label{circulo}
\end{figure}

\subsection{Example 3: Output on a manifold immersed in $\mathbb{R}^3$ }

Let us now consider a system with output in the following manifold $\mbox{$\mathcal{M}$}$ immersed in $\mathbb{R}^3$

\begin{displaymath}
\mbox{$\mathcal{M}$}= \left \{(x,y,z) \in \mathbb{R}^3 /  xyz=1, \, \, x,y,z >0  \right \}.
\end{displaymath}

More precisely, let $Z=f(X,Y)\mathrel{\mathop:}= (X+Y, \frac{1}{X}, \frac{X}{X+Y})$  where $X$ and $Y$ are iid random variables with Gamma distribution with parameters $(\mu_1,1)$, $(\mu_1>0)$. Since the function    $\mathbbm{1}_{\{z \leq w\}} = 0$ for all $z, w \in M, z\neq w$, the index estimation $\hat{C}^{\nu}_{2}$ based on the left lower quadrants introduced in \cite{gamboa2015}, does not provide any information about the sensitivity. Also the second moment of $1/X$ does not exist, therefore it is not possible to compute the Sobol' index defined in \cite{klein2014}. By varying the parameter $\mu_1$, we calculate the estimation  of the ball sensitivity index $B^{1}_{2}$ and $B^{2}_{2}$ . For this purpose,  $N=1000$ pick-freeze samples were generated for each value of $Z$, $(Z_j, Z^{\nu}_j)$, $j=1,\ldots,N$; and other $1000$ samples of $Z$,    independent of $(Z_j, Z^{\nu}_j)$ for the corresponding $W_k$.
In Figure \ref{ej3} the values of the indices are observed by varying the parameter $\mu_1$ between  $0$ and $5$. It is clear from Figure \ref{ej3} that the new index is able to detect the effect of each variable on the output for different values of the parameter $\mu$.

\begin{figure}[ht]\centering
\includegraphics[width=.80\textwidth]{ejemplo3.pdf}
\caption{Calculation of the index $\hat{B}^{\nu}_{2}$ }
\label{ej3}
\end{figure}

\section{Sensitivity for isotropic matrix}
\label{datosreales}
We consider now the stiffness matrix $Z$ of an  isotropic materials as a function of the constants of Lam\'e  $\lambda$ and $\mu$ (see \cite{landau1965}, p. 13). Hence, if the temperature is constant we have

\begin{equation*}
Z = \left(
\begin{array}{cccccc}
K+ 4\mu/3   & K-2\mu/3  & K-2\mu/3 & 0 & 0 & 0
 \\
K-2\mu/3  & K+ 4\mu/3  & K-2\mu/3 & 0 & 0 & 0 \\
K-2\mu/3  & K-2\mu/3  & K+ 4\mu/3 & 0 & 0 & 0 \\
0&0&0&\mu&0&0 \\
0&0&0&0&\mu&0 \\
0&0&0&0&0& \mu\\
\end{array} \right)
\end{equation*}

\noindent where $K= \lambda +2\mu/3$ is the volumetric modulus. The parameter $\mu$ is called the stiffness modulus, $K$ and $\mu$ take nonnegative values and are are modeled with a distribution supported by $\mathbb{R}^{+}$. The set of stiffness matrices is considered as a sub--manifold of the Riemannian manifold of the symmetric positive-definite matrices with the metric $g$. This manifold is usually denoted by $(\mathbb{P}_d,g)$, we refer to \cite{moakher2005} for more on the subject. Given two matrices $A$ and $B$  there exists an unique geodesic joining $A$ and $B$ given by,

\begin{equation}
\gamma(t) \mathrel{\mathop:}=  A ^{1/2} \left(A ^{-1/2} B A ^{-1/2}  \right) ^ t A ^{1/2}.
\end{equation}

\noindent So that, we can calculate the midpoint between $A$ and $B$. We will denote by $A\#B$ this midpoint. We have

\begin{align}
A\#B = A ^{1/2} \left(A ^{-1/2} B A ^{-1/2}  \right) ^{1/2} A ^{1/2} \\
& \hspace{-6.5 cm} d(A,B)= \Vert \textrm{ln}  \left(A ^{-1/2} B A ^{-1/2}  \right) \Vert,
\end{align}
where $\Vert \cdot \Vert$  is the Hilbert--Schmidt norm. Note that the midpoint is simply the geometric mean of the matrices.

\vspace{2 mm}

We will focus on the case when  $K$ and $\mu$ are independent random variables and on two scenarios:

\begin{tabular}{lll}
\textbf{Case 1} & $\quad K \sim \gamma( 1/\lambda_K, \lambda_K )$ & $ \mu \sim \gamma(1/\lambda_{\mu}, \lambda_{\mu} )$ \\
\textbf{Case 2} & $\quad K \sim \textrm{U}(1-\lambda_K , 1+\lambda_K )$ & $\mu \sim \textrm{U}(1-\lambda_{\mu},1+\lambda_{\mu})$ \\
\end{tabular}

Tables \ref{indice_ug1} and \ref{indice_ug2}  show the values of  $\hat{B}^{1}_{2}$ and $\hat{B}^{2}_{2}$  for  various values of $(\lambda_K,\lambda_{\mu})$. The sample size was $N=500$. In both tables we can see that the new index provides significant information regarding the sensibility of $Z$ with respect to the input variables $K$ and $\mu$.  In particular, when the variance of the input variables increases (larger $\lambda_{\mu}$ or $\lambda_{K}$) the value of the index increases, i.e. a a greater effect is observed on the output.

\begin{table}[ht]
\begin{center}
\caption{Index on isotropic matrix with Gamma distributions parameter}
\scalebox{0.8}{
\begin{tabular}{c | cccc|| cccc }
  \toprule[0.6 mm]
  Distribution  &  \multicolumn{4}{c}{Case 1: $B^1$ } & \multicolumn{4}{c}{Case 1: $B^2$ }             \\
\toprule[0.3 mm]
 $ \lambda_{\mu} \Big \backslash \lambda_K$ & \textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} &\textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} \\
  \toprule[0.3 mm]
  \textbf{0.001} & 0.625   & 0.212 & 0.016 & 0.001           & 0.083 & 0.435 & 0.855 & 0.980  \\
  \textbf{0.01} &0.925 & 0.593 & 0.215 & 0.033           & 0.004 & 0.072 & 0.458 & 0.865  \\
 \textbf{0.1} & 0.987 & 0.912 & 0.587 & 0.184          & 0.001 &0.006& 0.137 & 0.518         \\
\textbf{1} & 0.999 & 0.990 & 0.930 & 0.600       & 0.000 & 0.007 & 0.210 & 0.311   \\

  \toprule[0.6 mm]
\end{tabular}}
\label{indice_ug1}
\end{center}
\end{table}

\begin{table}[ht]
\begin{center}
\caption{ Index on isotropic matrix with Uniform distributions parameter}
\scalebox{0.8}{
\begin{tabular}{c | cccc || cccc }
  \toprule[0.6 mm]
  Distribution  &  \multicolumn{4}{c}{Case 1: $B^1$ } & \multicolumn{4}{c}{Case 1: $B^2$ }             \\
\toprule[0.3 mm]
  $ \lambda_{\mu} \Big \backslash \lambda_K$ & \textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} &\textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} \\
  \toprule[0.3 mm]
  \textbf{0.001} & 0.620   & 0.008  & 0.001 & 0.000          & 0.109 & 0.848 & 0.997 & 1.000  \\
  \textbf{0.01}  &0.989     & 0.623 & 0.003 & 0.001           & 0.018 & 0.092 & 0.849 & 0.997  \\
 \textbf{0.1} & 1.000      & 0.989 & 0.623 & 0.624          & 0.016 &0.016& 0.102 & 0.846        \\
\textbf{1} & 1.000         & 1.000 & 0.990 & 0.987            & 0.015 & 0.016 & 0.100 & 0.211     \\

  \toprule[0.6 mm]
\end{tabular}}
\label{indice_ug2}
\end{center}
\end{table}

\section{Conclusions and Future Work}
\label{conclusiones}

We introduce a new framework to measure the  global sensitivity  for system outputs valued on a Riemannian manifolds. This sensitivity index has nice properties:
\begin{enumerate}
\item It is a generalization of the one proposed by \cite{gamboa2015}.
\item  It is built on the whole distributions and not only on moments.

\item The construction of the index lies on the geometry of the support of the output through the  geodesic distance.

\item A pick-freeze like estimator of the index based on an $U$--statistic is  proposed. It is easy to compute. The consistency of the estimator is shown.

\item Considering  various simulation scenarios, the desirable properties of the estimator are illustrated. Using a Bootstrap resampling method for $U$--statistics we construct confidence intervals.

\item  The sensitivity index for an isotropic matrix model is studied.

\item In future research we will challenge to study new applications, for instance, the impact of inputs generated by a dynamic system, when the output is supported by a Riemannian Lie  Group.

\end{enumerate}

\section*{Appendix}
\begin{enumerate}[label=(\Alph*)]
\item \textbf{Proof of Property \ref{propiedad1}.}

In \cite{christensen1970} it is shown that both the sets

\begin{equation}\label{determ}
 \mathcal B \mathrel{\mathop:}= \{{B}_{p,r} \big/ p \in \mbox{$\mathcal{M}$}  \, \textrm{and} \, r>0 \}, \ \mbox{and} \ \mathcal B_{\delta} \mathrel{\mathop:}= \{{B}_{p,r} \big/ p \in \mbox{$\mathcal{M}$}  \, \textrm{and} \, 0<r<\delta \},
 \end{equation}
are determining classes. We aim to  show that the family of balls $$\mathcal{B}_{pq}:= \{B_{pq}, p,q \in \mbox{$\mathcal{M}$}\ \textrm{and there is a unique geodesic between $p$ and $q$} \},$$ is also a determining class if the radius of injectivity  $r_{iny}$ of $\mbox{$\mathcal{M}$}$ is positive.
Recall that in a metric space the family of closed sets $\mathcal{C}$ is a determining class (\cite{billingsley1999}, page $7$). Therefore, it suffices to show that if two probabilities $\nu$ and $\eta$ coincides on $\mathcal{B}_{pq}$, then $\nu(A)=\eta(A)$ for all $A \in \mathcal{C}$. Let  $\partial A$ the topological boundary of $A$ and $\partial A^{\epsilon} \mathrel{\mathop:}= \bigcup_{x \in \partial A} B(x, \epsilon)$ for $\epsilon >0$. So,

$$\nu (A)=  \nu \left( A \setminus (\partial A^{\epsilon} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon} ) \right) + \nu \left( A \cap (\partial A^{\epsilon} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon} ) \right).$$

Let  $\epsilon^{*} \mathrel{\mathop:}= \min (\epsilon/2, r_{iny})$.  Since by (\ref{determ}) $\mathcal B_{\epsilon^{*}}$ is a determining class, using the exponential map  we derive that the set of balls in $\mathcal B_{\epsilon^{*}}$
that determine the $\nu$-probability of $A \setminus (\partial A^{\epsilon^{*}} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon^{*}})$  are in
  $\mathcal{B}_{pq}$. Therefore, $$\nu \left( A \setminus (\partial A^{\epsilon^{*}} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon^{*}} ) \right)= \eta \left( A \setminus (\partial A^{\epsilon^{*}} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon^{*}}) \right).$$
Finally, from the dominated convergence theorem  and the assumptions, we conclude that
$$\nu (A)=  \eta(A)- \eta(\partial A  \cup \partial \mbox{$\mathcal{M}$}) + \nu(\partial A  \cup \partial \mbox{$\mathcal{M}$}) =\eta(A).$$

\QEDB

\item \textbf{Proof of Property \ref{exponencial}.}

We will show that given $s>0$, there is an $N_0$ such that for any $N>N_0$,
\begin{equation}
\label{exponencial2}
P \left( \left \vert \hat{S}^{\nu}_{2} -  S^{\nu}_{2} \right \vert > 9s   \right) \leq 16 \exp \left\{  -\frac{N^2}{8} \right\}.
 \end{equation}

We will prove that $$P \left(  \hat{S}^{\nu}_{2} -  S^{\nu}_{2}  > 9s   \right)  \leq 8 \exp \left\{  -\frac{Ns^2}{8} \right\},$$  and it is analogous for the other tail.

For $1 \leq j,k \leq N$ and $\tau \in \mathcal{P}_{N,2}$ let,
\begin{itemize}
\item $\textbf{W}_{\tau}= \left( W_{k_1},W_{k_2} \right)$
\item $\textbf{Z}_{j}= \left( Z_j, Z_j^{\nu}  \right)$
\item $G(\textbf{Z}_{j},\textbf{W}_{\tau})= h_{\textbf{W}_{\tau}} (Z_j) h_{\textbf{W}_{\tau}} (Z_j^{\nu})$
\item $J(\textbf{Z}_{j},\textbf{W}_{\tau})= \frac{1}{2} \left[  h_{\textbf{W}_{\tau}} (Z_j)+ h_{\textbf{W}_{\tau}} (Z_j^{\nu}) \right]$
\item  $H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau})=J(\textbf{Z}_{i},\textbf{W}_{\tau})J(\textbf{Z}_{j},\textbf{W}_{\tau})$
\end{itemize}

The proof is built on three steps:

\begin{itemize}
\item \textbf{Step 1} [Rewrite the difference $\hat{S}^{\nu}_{2}-S^{\nu}_{2}$]
Similarly to  \cite{gamboa2015}  but for an $U$--statistics, we have that

\begin{align*}
\hat{S}^{\nu}_{2} = \frac{1}{N\binom{N}{2}} \sum_{\substack{j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau})  \\
&\hspace{-11.7cm} =\frac{1}{N\binom{N}{2}} \sum_{\substack{j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right\}-\\
 &\hspace{-11.7cm}  -\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right \}  +  \\
&\hspace{-11.7cm} + \mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{1},\textbf{W}_{\tau_1}) \right] - \left( 1-\frac{1}{N}\right) \mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{2},\textbf{W}_{\tau_1}) \right]  -   \frac{1}{N}\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{1},\textbf{W}_{\tau_1}) \right],
\end{align*}
and
\begin{align*}
S^{\nu}_{2}= \mbox{$\mathbb{E}$}_{\textbf{W}_1} \left\{ \textrm{Var}_{X_{\nu}}\left( H^{\nu}(\textbf{W}_1)\right) \right \}&\\
&\hspace{-4.5cm}  =  \mbox{$\mathbb{E}$}_{\textbf{W}_1} \left\{  \textrm{Var}_{X_{\nu}}\left( \mbox{$\mathbb{E}$}_Z (h_{\textbf{W}_1} (Z)/ X_{\nu}) \right) \right \}=  \mbox{$\mathbb{E}$}_{\textbf{W}_{1}} \left\{  \textrm{cov} \left( h_{\textbf{W}_{1}} (Z_1), h_{\textbf{W}_{1}} (Z_1^{\nu}) \right) \right \} \\
&\hspace{-4.5cm}= \mbox{$\mathbb{E}$}_{\textbf{W}_{1}} \left\{  \mbox{$\mathbb{E}$}_{\textbf{Z}_1} \left( h_{\textbf{W}_{1}} (Z_1) h_{\textbf{W}_{1}} (Z_1^{\nu}) \right) \right \}-  \mbox{$\mathbb{E}$}_{\textbf{W}_{1}}\left\{ \left[  \mbox{$\mathbb{E}$}_{\textbf{Z}_1}  \left( h_{\textbf{W}_{1}} (Z_1) \right)\right]^2 \right \}.
\end{align*}

Now, we may  decompose the error   into three terms
\begin{align*}
\hat{S}^{\nu}_{2} -  S^{\nu}_{2}= \underbrace{\frac{1}{N\binom{N}{2}} \sum_{\substack{j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right\}}_{(A)}-\\
 &\hspace{-9cm} \underbrace{ -\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right \}  }_{(B)} + \\
&\hspace{-9cm} \underbrace{+\frac{1}{N}  \left \{ \mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{2},\textbf{W}_{\tau_1}) \right]  - \mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{1},\textbf{W}_{\tau_1}) \right] \right \}.}_{(C)}
\end{align*}

In the following steps, we bounded the terms $(A)$ and $(B)$ for a sufficiently large $N$.

In  steps  $2$ and $3$ we use Hoeffding inequality (see for example \cite{hoeffding1963}) for independent and real variables. That is, if  $W_1, \ldots, W_N$ are  independent centred random variables, supported by $[a_i,b_i] \quad i=1, \ldots,N$,
\begin{equation}
P\left( \sum_{i=1}^{N} W_i > s \right) \leq  \exp \left\{ - \frac{2s^2}{\sum_{i=1}^{N}(b_i-a_i)^2} \right\} \quad (s>0).
\end{equation}

We also use an extension of the previous inequality for $U$--statistics of order $2$, see \cite{serfling1980}, Theorem A, (5.6).  If  $s=s(X_1, X_2)$ is the $U$--statistics kernel of $U_n$ such that  $\mbox{$\mathbb{E}$}(s(X_1,X_2))=\theta$ and  $a \leq s(x_1, x_2) \leq b$. So, for  $s>0$ and  $N>2$,
\begin{equation}
P\left( U_n -\theta > s \right) \leq  \exp \left\{ - \frac{Ns^2}{(b-a)^2} \right\}
\end{equation}
For us, the kernels are centred bounded indicator functions.

\item \textbf{Step 2} [Bounds for $(A)$.]

Let $G_c(\textbf{Z}_{j},\textbf{W}_{\tau})= G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{j},\textbf{W}_{\tau}) \right]$ and $\tilde{G}(Z_i, \textbf{W}_{\tau})= G_c(Z_i, \textbf{W}_{\tau})  - \mbox{$\mathbb{E}$}_{Z_i} (G_c (Z_i, \textbf{W}_{\tau}))$,
\begin{displaymath}
A= \underbrace{\frac{1}{N\binom{N}{2}} \sum_{\substack{ i \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }}  \tilde{G}(Z_i,\textbf{W}_{\tau})}_{A_1} + \underbrace{\frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ G_c(Z,\textbf{W}_{\tau}) \right] }_{A_2}
\end{displaymath}

\begin{align*}
P(A_1>s)= P \left(\frac{1}{N\binom{N}{2}} \sum_{\substack{ i \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }}  \tilde{G}(Z_i,\textbf{W}_{\tau}) >s \right) \\
&\hspace{-7 cm} = P \left( \sum_{i=1}^N  \frac{1}{\binom{N}{2}}\sum_{\tau \in \mathcal{P}_{N,2} } \tilde{G}(Z_i,\textbf{W}_{\tau}) >Ns\right)\\
&\hspace{-7 cm} \leq \exp \left\{  -\frac{s^2N}{2} \right\} \leq  \exp \left\{ -\frac{Ns^2}{8} \right\},
\end{align*}
while
\begin{align*}
P(A_2>s)= P \left( \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ G_c(Z,\textbf{W}_{\tau}) \right] >s \right) \\
&\hspace{-7 cm} \leq   \exp \left\{\frac{-2Ns^2}{8} \right\}\leq \exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}

But, $\left\{ A_1 + A_2> 2s \right\} \subset \{A_1>s\} \cup \{A_2>s \} $, so
\begin{equation}
P(A>2s) \leq 2\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{equation}

\item \textbf{Step 3} [Bounds for $(B)$.]

Let $H_c(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) = H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) \right]$ and  $\tilde{H}(Z_i, Z_j , \textbf{W}_{\tau})= H_c (Z_i,Z_j, \textbf{W}_{\tau})  - \mbox{$\mathbb{E}$}_{Z_j} (H_c (Z_i, Z_j,\textbf{W}_{\tau}) ) $,

\begin{align*}
B= \underbrace{\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }}  \tilde{H}(Z_i,Z_j, \textbf{W}_{\tau})}_{B_1}+ \\
&\hspace{-6.5cm}  + \underbrace{\frac{1}{N\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z,\textbf{W}_{\tau}) \right] }_{B_2} + \\
&\hspace{-6.5 cm}+ \underbrace{\frac{N-1}{N^2\binom{N}{2}} \sum_{\substack{ j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} } } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] }_{B_3}.
\end{align*}

We have $\left\{ \sum_{i=1}^3B_i> 6s \right\} \subset \{B_1>3s\} \cup \{B_2>s \}\cup \{B_3>2s \} $.

So that,
\begin{align*}
P(B_1>3s) \leq P \left(\sum_{j} \frac{1}{\binom{N}{2}}  \sum_{ \tau \in \mathcal{P}_{N,2} }  \tilde{H}(Z_j,Z_j,\textbf{W}_{\tau}) >N^2s \right)+ \\
&\hspace{-8.5 cm}+ P \left( \sum_{j>i} \frac{1}{\binom{N}{2}}  \sum_{ \tau \in \mathcal{P}_{N,2} }   \tilde{H}(Z_i,Z_j,\textbf{W}_{\tau}) >N^2s  \right)+ \\
&\hspace{-8.5 cm}+ P \left( \sum_{j<i} \frac{1}{\binom{N}{2}} \sum_{ \tau \in \mathcal{P}_{N,2} }   \tilde{H}(Z_i,Z_j,\textbf{W}_{\tau}) >N^2s \right)  \leq \\
&\hspace{-8.5cm} \leq \exp \left\{  -\frac{s^2N^3}{2} \right\} + 2 \exp \left\{  \frac{-Ns^2N^2}{4} \right\}\leq 3\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}
Further,
\begin{align*}
P(B_2>s)  \leq P \left( \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z,\textbf{W}_{\tau}) \right]  > Ns \right) \\
&\hspace{-8 cm} \leq  \exp \left\{  -\frac{Ns^2N^2}{4} \right\} \leq \exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}

Finally,
\begin{align*}
P(B_3>2s) \leq P \left( \frac{N-1}{N^2\binom{N}{2}} \sum_{\substack{ j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} } } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] >2s \right)  \\
&\hspace{-10.1 cm} \leq P \left( \sum_{ j \in \{1,\ldots, N\}  } \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} }\mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] >2s \frac{N^2}{N-1}\right) \\
&\hspace{-10.1 cm}  \leq P \Bigg( \sum_{ j \in \{1,\ldots, N\}  } \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \Big \{ \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right]- \\
&\hspace{-10.1 cm} -\mbox{$\mathbb{E}$}_{Z,Z_j} \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] \Big \} >\frac{s N^2}{N-1}\Bigg) + \\
&\hspace{-10.1 cm} + P \left(  \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} }\mbox{$\mathbb{E}$}_{Z_1,Z_2} \left[ H_c(Z_1,Z_2,\textbf{W}_{\tau}) \right]  >s \frac{N}{N-1}\right) \\
&\hspace{-10.1 cm} \leq  \exp \left\{  -\frac{s^2N^4}{2N(N-1)^2} \right\} + \exp \left\{  -\frac{2Ns^2 N^2}{8(N-1)^2} \right\} \leq 2\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}

Collecting the previous inequality we obtain,
\begin{equation}
P(B>6s) \leq 6\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{equation}
\end{itemize}

\noindent  Now, using the bounds obtained in steps $2$ and $3$, we may conclude that,  for large enough  large $N$,  inequality (\ref{exponencial2}) holds.

\QEDB

\item \textbf{Proof to Example \ref{ejemplo1}.}

Let  $P(z_1, z_2)= F(z_2)-F(z_1)$, then we have
 \begin{align*}
H(s,t)=\mathbbm{1}_{t <s} +F_Z(t) -F_Z(s))\\
 &\hspace{-3.4cm} =  \mathbbm{1}_{t <s} +(1-p)P(s,t)+pP(s-\alpha,t-\alpha),
 \end{align*}
 and,
\begin{displaymath}
H^1(s,t)=\begin{cases}
\mathbbm{1}_{t <s} + P(s,t) & \text{if } X_1=0,\\
\mathbbm{1}_{t <s} +P(s-\alpha,t-\alpha) & \text{if } X_1=1.
\end{cases}
\end{displaymath}
Therefore
\begin{small}
\begin{displaymath}
H(s,t)-H^1(s,t)=\begin{cases}
p \left[ P(s,t)-P(s-\alpha,t-\alpha)\right] & \text{if } X_1=0,\\
(1-p)\left[P(s-\alpha,t-\alpha)- P(s,t)\right]& \text{if } X_1=1,
\end{cases}
\end{displaymath}
\begin{displaymath}
\mbox{$\mathbb{E}$}_1 \left[\left(H(s,t)-H^1(s,t)\right)^2 \right] =p(1-p)\left[ P(s,t)-P(s-\alpha,t-\alpha)\right]^2,
\end{displaymath}
\end{small}
and,
\begin{align*}
S^{1}_{2}=p(1-p) \mbox{$\mathbb{E}$}\left[ \left(P(Z_1-\alpha,Z_1)-P(Z_2-\alpha,Z_2) \right)^2 \right]\\
 &\hspace{-6.9cm} = 2p(1-p) \textrm{Var} \left[ F(Z)-F(Z-\alpha) \right].
\end{align*}

Let  $X_2 \sim U(0,b)$, with $b= \sqrt{12\alpha^2 p(1-p)}$. We  obtain

\begin{displaymath}
\mbox{$\mathbb{E}$}\left(  F(Z)- F(Z-\alpha)  \right)=\begin{cases}\lim_{N\rightarrow \infty} \frac{1}{N^2}
  \frac{\alpha}{b}   \left(1- \frac{1}{2}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/2 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}

\begin{displaymath}
\mbox{$\mathbb{E}$}\left[ \left( F(Z)- F(Z-\alpha) \right)^2 \right]=\begin{cases}
 \left( \frac{\alpha}{b} \right)^2  \left(1- \frac{2}{3}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/3 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}

\begin{displaymath}
\textrm{Var} \left( F(Z)- F(Z-\alpha) \right)=\begin{cases}
 \left( \frac{\alpha}{b} \right)^3  \left(\frac{1}{3}- \frac{1}{4}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/12 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}

So,

\begin{displaymath}
D^{1}_{2,CVM}= p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^2  \left(1- \frac{2}{3}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/3 & \text{if } \alpha > b,
\end{cases}
\end{displaymath}

and

\begin{displaymath}
S^{1}_{2}= 2p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^3  \left(\frac{1}{3}- \frac{1}{4}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/12 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}

\end{enumerate}

\QEDB


\section{Introduction}

After formulating the \textit{cosmic censorship conjecture},
Penrose proposed \cite{Penrose1973} that, when considering
collapsing matter,  if the conjecture is valid, the mass $m$ and the area $A$ of the resulting black hole must satisfy the relation:

\begin{equation}\label{penroseineq}
m\geq \sqrt{\frac{A}{16\pi}}
\end{equation}

Given that the validity of this relation gives one of the most important tests to evaluate the cosmic censorship conjecture, there are plenty of works that study this relation, and one can find very exhaustive review articles \cite{Mars2009, Bray2003} as well as general approaches to the problem  \cite{Frauendiener:01, BrayKhuri:11, BrayKhuri:09}. Moreover, since the Penrose heuristic argument one can strengthen the to include charge and angular momentum and charge (see \cite{Dain:2014Geo} \cite{Dain:2012} \cite{Mars2009} for more details), this topic has become an active area of research.
Good progress has been made in considering the case of a charged black hole without angular momentum \cite{Weinstein:2005,Disconzi:2012, KhuriYamada:2013, Khuri:2013, Khuri:2014, Khuri:2015}, but regarding the case with angular momentum, there are only a few results exploring the relation between the angular momentum, the size and the mass of a compact object \cite{Anglada:2016dbu, Anglada:2017ryp, Jaracz-Khuri18}. For further details and references regarding geometrical inequalities bounding angular momentum see the review article \cite{Dain:2017jkj}. In this work we present an extension of our previous work \cite{Anglada:2017ryp} for compact and connected general horizons.

Take an axially symmetric initial data $M$ connecting the black hole region with spatial infinity , such that the collapse has already occurred, and calculate the mass $m$, the area $A$, and the angular momentum $J$ of the black hole. Then, from the Penrose heuristic argument for rotating black holes, see for example \cite{Anglada:2017ryp}, we expect that:
\begin{equation} \label{PenJ}
m^2 \geq \frac{A }{16 \pi} + \frac{4 \pi J^2}{A}
\end{equation}
Note that this version of the Penrose inequality admits a rigidity case which states that the equality can only occur for the Kerr black hole.

In \cite{Anglada:2017ryp} we studied this problem in the particular case that the apparent horizon is a compact and connected minimal surface.
We use the monotonicity properties of the Geroch energy \cite{Geroch1973} along the IMCF and proved the following version of \eqref{PenJ}:
\begin{equation}\label{PenJms}
m_{ADM}^2\geq \frac{A }{16 \pi} +  \frac{ J^2}{ 2\mathcal{R}^2}
\end{equation}
where $m_{ADM}$ is the ADM mass, \cite{ADM}, and $\mathcal{R}$ is a specific measure of size defined in terms of the norm of the axial Killing vector. This measure has reasonably nice properties, see \cite{Anglada:2016dbu, Anglada:2017ryp}, and under certain conditions can be related to usual measures of size.
In this previous work, in order to have a non-negativity scalar curvature we need to assume some special conditions for the extrinsic curvature, for example that the initial data is maximal.
In this work we use Hawking energy, instead of Geroch energy, similar techniques and the same measure of size to extend and improve \eqref{PenJms} for a compact and connected general horizon.

Regarding the mentioned papers, \cite{Anglada:2016dbu, Anglada:2017ryp, Jaracz-Khuri18}, that explore the relation between the angular momentum, the size and the mass of a compact object it is important to remark the following issue. Although the technique used to relate the angular momentum with the surface integral of the extrinsic curvature are very similar, and the particular measure of the axial radii is the same, in \cite{Anglada:2016dbu}, \cite{ Anglada:2017ryp}, and in the present work, one need to impose some particular conditions over the extrinsic curvature in order to assure that the energy is non-decreasing along the flow, and also one need to assume there exist a smooth solution of the IMCF for the initial data. In \cite{Jaracz-Khuri18} the author present a different approach to this problem, they use an embellished version of the Jang equation and they study a Jang/IMCF system of equations. This allows the authors to obtain a result very similar to the one presented on this work but assuming only a few general conditions, but they need assume that there exist a smooth solution of a more complex system of equations, where the IMCF equation and a particular version of Jang equation are coupled.

\section{Background}
\label{Background}

We consider an \textit{asymptotically flat} and \textit{axially symmetric} initial data set $(M, \partial M,\bar g, K; \mu, j^i)$ with boundary $\partial M$, where $M$ is a 3-manifold with positive definite metric $\bar g$ and extrinsic curvature $K$, $\partial M$ is a connected and compact 2-surface, $\mu$ is the energy density and $j^i$ is the matter current density. This set must satisfy the constraint equations
\begin{align}
 \label{const1}
   \bar D_j   K^{i j} -  \bar D^i   k= -8\pi j^i,\\
 \label{const2}
   \bar R -  K_{i j}   K^{i j}+  k^2=16\pi \mu,
\end{align}
where ${\bar D}$ and $\bar R$ are the Levi-Civita connection and the curvature scalar associated with $  \bar g$, and $k=\mbox{tr}_{\bar g}K$.
We assume the matter fields satisfy the \textit{Dominant Energy Condition} (DEC), $\mu \geq |j|$,  and that $\partial M$ is an outermost future apparent horizon and there are no other trapped surfaces on $M$. With these assumptions $M$ is an \textit{exterior region} and has the topology $\mathbb{R}^3$ minus a ball \cite{Huisken2001}.

Assume there exists a smooth inverse mean curvature flow (IMCF) of surfaces $S_t$ starting from $S_{0}=\partial M$ and having spherical topology, and take $\nu$ to be unit normal vector of $S_t$ (see equation \eqref{eqIMCF} on Appendix A). Then one can write the metric $\bar g$ in the form:
\begin{equation}\label{barg}
ds^2_{\bar g}=\frac{dt^2}{H^2} + g_{ij}dx^i dx^j
\end{equation}
where $H$ and $g_{ij}$ are the mean curvature and the induced metric of $S_t$ respectively, and $(x^1, x^2)$ are and the induced coordinates.
See \cite{Huiskenevol} \cite{Szabados04} for a review of the basic properties of the IMCF.
In this context the extrinsic curvature can be decomposed \cite{Malec:2002ki}:
\begin{equation}\label{K}
K_{ij}=z\nu_i \nu_j + \nu_i s_j + s_i \nu_j  + g^k_i g^l_j \chi_{lk} + \frac{q}{2}g_{ij}
\end{equation}
where q is the trace with respect to $g_{ij}$ of $K$, $q=K_{ij}g^{ij}$ and
\begin{equation}\label{Kparts}
z=K_{ij}\nu^i \nu^j \qquad s_i=g^j_i K_{jl} \nu^l \qquad \chi_{ij}=g^l_i g^n_j K_{ln} - \frac{q}{2} g_{ij}
\end{equation}
then the trace of the extrinsic curvature takes the form $k=\mbox{tr}(K)=z+q$ and its norm is
\begin{equation}\label{KK}
K_{ij}K^{ij}=z^2+2s_i s^i + \chi_{ij} \chi^{ij} + \frac{q^2}{2}.
\end{equation}
and then from equation \eqref{const2} we have that the scalar curvature of the initial data can be written in the following way:
\begin{equation}
\label{barR}
\bar  R= 16\pi \mu + 2s_i s^i + \chi_{ij} \chi^{ij} -  \frac{q}{2}(q + 4z).
\end{equation}

Let $\vartheta^+$ and $ \vartheta^-$ be the expansions of the outgoing null geodesics orthogonal to $S_t$, future directed and past directed respectively, then because $\partial M$ is a future apparent horizon we have $\vartheta^+|_{\partial M}=0$ and $\vartheta^-|_{\partial M}\geq0$.
From the previous decomposition $ \vartheta^+|_{S_t}=H+q$ and  $\vartheta^-|_{S_t}=H-q$ and then
if $M$ has no other trapped surface than $\partial M$, the expansions satisfy satisfy $\vartheta^+|_{S_t} \geq 0$, $\vartheta^-|_{S_t}\geq 0$ $\forall t >0$, and then
\begin{equation}\label{NTS}
(\vartheta^+\vartheta^-)|_{S_t}=H^2-q^2>0 \quad \forall t > 0,
\end{equation}

Following \cite{Malec:2002ki} we are going to use a functional proposed by Hawking \cite{Hawking:1968qt},the Hawking energy of a surface $E_H(S)$:
\begin{equation} \label{Hmass-gral}
E_H(S):=\frac{A^{1/2}}{(16\pi)^{3/2}}\left(16\pi-\int_{S}\vartheta^+\vartheta^-dS\right)
\end{equation}
where $A$ is the area of $S$.
This energy, under certain conditions (see \cite{Malec:2002ki}) is monotonic under a smooth inverse mean curvature flow, and has the interesting properties that it tends to the ADM mass of $M$ at infinity and for a future apparent horizon is equal to $\sqrt{\frac{A}{16 \pi}}$.
From \cite{Malec:2002ki} we have that the derivate of Hawking energy along the IMCF can be written in the following way (see appendix \ref{sec:IMCF} for more details and definitions of used quantities):
\begin{multline}
\label{dtEgralf}
\frac{d}{dt}E_H=\frac{A_t^{1/2}}{(16\pi)^{3/2}}\left[ \int_{S_t} \left(16 \pi ( \mu +\frac{q}{H}\nu^i j_i)  \right) + \int_{S_t} 2 \frac{q}{H} g^{ij} \bar \nabla_i s_j dS
\right.\\
\left.              +\int_{S_t} \left(\chi_{ij} \chi^{ij} -  2\frac{q}{H}\chi^{ij} t_{ij}  + t_{ij}   t^{ij}\right)dS \right.\\
\left.              +\int_{S_t} 2\left(s_i s^i - 2 \frac{q}{H} s^j \frac{ \bar \nabla_j H}{H} + \frac{ g^{ij} \bar \nabla_i H \bar \nabla_j H}{H^2} \right)dS \right]
\end{multline}

\section{Main result}

Now from \cite{Anglada:2016dbu} we know that when considering the IMCF in axially symmetric initial data,
the IMCF equation preserves axial symmetry. Then from now on, when we discuss the IMCF flow, we always consider it consisting of axially symmetric surfaces $S_t$. Then for each surface of the flow we can define orthogonal coordinates $\theta, \varphi$ such that $\eta^i=\frac{\partial}{\partial \varphi}^i$. One can always choose this for axially symmetric 2-surfaces that are diffeomorphic to $S^2$, see for example
\cite{Dain:2011pi}. Hence we have:
\begin{equation} \label{gSt}
  ds^2_g=\Psi^4 d\theta^2+ \eta d\varphi^2
\end{equation}
where $\eta=g_{ij}\eta^i \eta^j$ is the square norm of the axial Killing vector.

The physical and geometrical quantities we are interested in are the ADM mass $m_{ADM}$ and the Komar angular momentum $J(S_t)$:
\begin{equation}\label{defangmom}
 J(S_t)=\frac{1}{8\pi}\int_{S_t} K_{ij} \eta^i \nu^jdS,
\end{equation}
where we use that $\bar g_{ij} \nu^i \eta^j=0$.

To measure the size of the surface $S_t$ we will use the areal and circumferential radii of a surface $S_t$ in $M$:
\begin{equation}\label{size}
\mathcal{R}_A(S_t):=\sqrt{\frac{A_t}{4\pi}},\qquad \mathcal{R}_C(S_t):=\frac{\mathcal C (S_t)}{2\pi}
\end{equation}
where $A_t$ is the area of $S_t$ and $\mathcal C(S_t)$ is the length of the greatest axially symmetric circle in $S_t$.
It is also useful to consider the following size measure studied in \cite{Anglada:2016dbu, Anglada:2017ryp}:
\begin{equation}\label{R}
\frac{1}{\mathcal{R}(S_t)^2} := A_t^{1/2} \int_{t}^\infty \frac{A_{t'}^{1/2}}{\int_{S_{t'}}  \eta dS }dt'
\end{equation}

This measure of size of a surface $S_t$, based on the behavior of the norm of the Killing vector along the IMCF from $S_t$ to infinity, is positive and well defined provided the flow remains smooth. Moreover, as shown in \cite{Anglada:2017ryp}, in some cases, $\mathcal{R}$ can be related to $R_A$ and  $R_C$. In particular assuming that the IMCF is convex we have:
\begin{equation}
\mathcal{R}^2(S_t) \leq \frac{5}{2}\frac{\int_{S_t} \eta dS}{A_t} \leq \frac{5}{2} \mathcal{R}_C^2(S_t),
\end{equation}

Using the previous tools and this definition of size, and assuming the same conditions of the main theorem in \cite{Malec:2002ki}, we prove the following theorem.

\begin{theorem}
\label{theo1}
Let $(M,\partial M, \bar g, K)$ be a vacuum, asymptotically flat, and axially symmetric initial data, such that $\partial M$ is a compact and connected outermost apparent horizon and there are no other trapped surfaces on $M$. Assume there exists a smooth IMCF of surfaces $S_t$ starting from $\partial M$ and having spherical topology, then if the initial data satisfies either:
\begin{enumerate}[a)]
        \item \label{a} $\frac{q}{H}$ is constant for each surface $S_t$ , or
        \item \label{b} $g^{ij} \bar \nabla_i s_j=0$,
\end{enumerate}
then:
\begin{equation}\label{teo1}
m_{ADM}^2 \geq \frac{A}{16 \pi} + \frac{ J^2}{\mathcal{R}^2}
\end{equation}
where $J$ and $A$ are the angular momentum and the area of $\partial M$ respectively, and $\mathcal{R}=\mathcal{R}(\partial M)$ is defined by \eqref{R}.
\end{theorem}

\begin{proof}

From the fact that the surfaces of the flow are axially symmetric we have that mean curvature $H$ does not depend on the coordinate $\varphi$:
\begin{equation}
s^j \bar \nabla_j H = s^\theta \bar \nabla_\theta H \qquad  g^{ij} \bar \nabla_i H \bar \nabla_j H= g^{\theta \theta} \bar \nabla_\theta H \bar \nabla_\theta H
\end{equation}
hence in this case:
\begin{multline}
\frac{d}{dt}E_H=\frac{A_t^{1/2}}{(16\pi)^{3/2}}\left[ \int_{S_t} 2 \frac{q}{H} g^{ij} \bar \nabla_i s_j dS
+\int_{S_t} \left(\chi_{ij} \chi^{ij} -  2\frac{q}{H}\chi^{ij} t_{ij}  + t_{ij}   t^{ij}\right)dS
\right.\\
\left.    +\int_{S_t} 2\left(s_\theta s^\theta - 2 \frac{q}{H^2} s^\theta  \bar \nabla_\theta H + \frac{1}{H^2} g^{\theta\theta} \bar \nabla_\theta H \bar \nabla_\theta H  \right)dS \right] + 2\frac{A_t^{1/2}}{(16\pi)^{3/2}}\int_{S_t}  s_\varphi s^\varphi dS
\end{multline}
where we have used that the initial data is vacuum.

In order to include the angular momentum into the inequality we know from \cite{Anglada:2017ryp} that we can relate the angular momentum of any surface $S_t$ to the surface integral of the norm of $s_i$ and the norm of $\eta^i$.  In this work we need to improve the previous calculation in order to relate the angular momentum, not with the norm of $s_i$, but only with the component of $s_i$ along the axial Killing vector. First note that $K_{ij} \eta^i \nu^j=s_i \eta^i=s_\varphi$, then using the Cauchy-Schwarz inequality in the definition of $J_t:=J(S_t)$:
\begin{equation}\label{angmom}
\begin{split}
 J_t^2&=\frac{1}{(8\pi)^2} \left( \int_{S_t} s_i \eta^i dS\right)^2 = \frac{1}{(8\pi)^2} \left( \int_{S_t} s_\varphi dS\right)^2 \\
 &\leq \frac{1}{(8\pi)^2} \left( \int_{S_t} |s_\varphi| \frac{\sqrt{ \eta}}{\sqrt{\eta}} dS\right)^2 \leq \frac{1}{(8\pi)^2} \int_{S_t} \frac{s_\varphi^2}{\eta}dS  \int_{S_t}\eta dS \\
 &= \frac{1}{(8\pi)^2} \int_{S_t} s_\varphi s^\varphi dS  \int_{S_t}\eta dS
\end{split}
\end{equation}
where in the fourth step we have used the H\"older inequality with $p_1=p_2=2$. Then we have:
\begin{equation}
\int_{S_t} s_\varphi s^\varphi dS \geq \frac{(8\pi)^2 J_t^2}{\int_{S_t}\eta dS}
\end{equation}
and hence we could include explicitly the angular momentum on the derivate of Hawking energy:

\begin{multline}\label{dtEvsJgral}
\frac{d}{dt}E_H \geq\frac{A_t^{1/2}}{(16\pi)^{3/2}}\left[ \int_{S_t} 2 \frac{q}{H} g^{ij} \bar \nabla_i s_j dS
+\int_{S_t} \left(\chi_{ij} \chi^{ij} -  2\frac{q}{H}\chi^{ij} t_{ij}  + t_{ij}   t^{ij}\right)dS
\right.\\
\left.    +\int_{S_t} 2\left(s_\theta s^\theta - 2 \frac{q}{H^2} s^\theta  \bar \nabla_\theta H + \frac{1}{H^2} g^{\theta\theta} \bar \nabla_\theta H \bar \nabla_\theta H \right)dS \right] +2\sqrt{\pi} J_t^2 \frac{A_t^{1/2}}{\int_{S_t}\eta dS}
\end{multline}

Now first note that from the hypothesis that there are no other trapped surfaces in $M$ than $\partial M$ we have that $H^2>q^2$ and thus $|\frac{H}{q}| \geq 1$, hence the integrands on the second  and third terms $$\left(\chi_{ij} \chi^{ij} -  2\frac{q}{H}\chi^{ij} t_{ij}  + t_{ij}   t^{ij}\right), \quad \left( s_\theta s^\theta - 2 \frac{q}{H^2} s^\theta  \bar \nabla_\theta H + \frac{1}{H^2} g^{\theta\theta} \bar \nabla_\theta H \bar \nabla_\theta H\right) $$ are positive quadratic forms, and thus the second and third integrals in \eqref{dtEvsJgral} are positive. Then note that if we assume condition \ref{b} the integrand on the first term in \eqref{dtEvsJgral} is equal to zero, and if we assume condition \ref{a} then by partial integration on $S_t$ the first term in \eqref{dtEvsJgral} also vanishes.

Thus assuming the hypothesis of the theorem and either of the conditions \ref{a} or \ref{b}, the first term in \eqref{dtEvsJgral} vanishes and the second and third terms are positive, hence:
\begin{equation}\label{dtEvsJ}
\frac{d}{dt}E_H\geq \sqrt{4\pi} J^2 \frac{A_t^{1/2}}{\int_{S_t}\eta dS}
\end{equation}
where we have used that $M$ is a vacuum exterior region, thus $J_t=J$.

From these arguments we have that $E_H$ is monotonically increasing along the flow, hence $E_H(S_t)\geq E_H(\partial M) \quad \forall t\geq 0 $, then because $\partial M$ is an apparent horizon $H^2=q^2$ we have $E_H(\partial M)=\sqrt{\frac{A}{16\pi}}$ and thus:
\begin{equation}
\label{Penrose}
E_H(S_t)\geq \sqrt{\frac{A}{16\pi}} \quad \forall t\geq 0
\end{equation}

Now we calculate the derivate along the flow of the functional $E_H^2$ and use equation \eqref{dtEvsJ} to obtain a lower bound for it in terms of $J$:
\begin{equation}
\frac{d}{dt}E_H^2 =2E_H(S_t) \frac{d}{dt}E_H  \geq   2E_H(S_t) \sqrt{4\pi} J^2 \frac{A_t^{1/2}}{\int_{S_t}\eta dS}
\end{equation}
then using equation \eqref{Penrose} we have:
\begin{equation}\label{dtE2vsJ}
\frac{d}{dt}E_H^2  \geq    J^2 \sqrt{A} \frac{A_t^{1/2}}{\int_{S_t}\eta dS}
\end{equation}
Now integrating this expression along the flow from $\partial M$ to infinity and using the relation between Hawking energy and the ADM mass we obtain:
\begin{equation}\label{evol4}
m_{ADM}^2\geq\lim_{t \to \infty} E_H^2(S_t) \geq E_H^2(S_0) + J^2 \sqrt{A} \int_0^\infty \frac{A_t^{1/2}}{\int_{S_t}  \eta dS}dt.
\end{equation}

Finally we use the fact that $E_H^2(S_0)=\frac{A}{16\pi}$, and we write this expression in terms of $\mathcal{R}$ and obtain \eqref{teo1}.

\end{proof}

Inequality \eqref{teo1} is also valid for non-vacuum initial data, provided that the matter fields satisfy the DEC and that $j_i\eta^i=0$ everywhere in $M$.
Assuming the DEC we assure that Hawking energy remains monotonic for non-vacuum initial data. Condition $j_i\eta^i=0$ assures that the angular momentum is preserved along the flow $J(S_t)=J$, and that there is no contribution to $J$ coming from the matter fields $J=J(\partial M)$.

In case we have a non-zero contribution of the matter fields to the angular momentum, $j_i\eta^i \neq0$, we obtain an extension, for objects that contain a general horizon, of the results for ordinary objects presented in \cite{Anglada:2016dbu}. We assume that the matters fields satisfy the DEC and that both the matter density and the matter current have compact support.
In this case the angular momentum of a surface $S_t$ is
\begin{equation}
 J(S_t)=\frac{1}{8\pi}\int_{S_t} K_{ij} \eta^i \nu^jdS=  J(\partial M) - \int_{V(S_t)} j_{i} \eta^i dv,
\end{equation}
where $V(S_t)$ is the region enclosed between $\partial M$ and $S_t$. Thus the conservation of the angular momentum along the flow  is only satisfied when the surfaces $S_t$ are outside the compact support of the matter fields. Then the measures of size involved in the rotational contribution to the energy are not measures of size of the apparent horizon, but measures of size of the first surface of the flow $S_T$ that enclosed the object.

Then for a non-vacuum initial data with $j_i\eta^i \neq0$, assuming the same conditions of theorem \eqref{theo1}, we obtain the following result.

\begin{theorem}

Let $(M,\partial M, \bar g, K; \mu, j^i)$ be an initial data satisfying the same conditions of theorem  \ref{theo1}. Assume the matter fields satisfy the dominant energy condition and have compact support, and let $T$ such that for all $t\geq T$ the matter density and the matter current have compact support inside $S_t$, then:

\begin{equation}\label{teo2}
m_{ADM} \geq  m_T + \frac{\mathcal{R}_A}{2} + \frac{ J^2}{\mathcal{R}_A(T)\mathcal{R}^2(T)}
\end{equation}
where $J$ is the total angular momentum of the data, $\mathcal{R}_A$ and $\mathcal{R}_A$(T) are the areal radii of $\partial M$ and $S_T$ respectively, $\mathcal{R}(T)=\mathcal{R}(S_T)$ is defined  by \eqref{R}, and
\begin{equation}\label{mTbh}
m_{T}:=\int_{\mathcal{R}_A}^{\mathcal{R}_A(T)}\int_{S_\xi} \left( \mu +\frac{q}{H}\nu^i j_i \right) dS d\xi
\end{equation}
where $\xi$ stands for the  areal radius coordinate.
\end{theorem}

\begin{proof}
From \eqref{dtEgralf} and the previous calculations we have:
\begin{equation}\label{dtEvsJ2}
\frac{d}{dt}E_H\geq  \sqrt{\frac{A_t}{16\pi}} \int_{S_t} \left( \mu +\frac{q}{H}\nu^i j_i  \right) dS +\sqrt{4\pi} J_t^2 \frac{A_t^{1/2}}{\int_{S_t}\eta dS}
\end{equation}
Note that because the matter fields satisfy the DEC, and $|\frac{q}{H}|\leq 1$ the fist term in \eqref{dtEvsJ2} is also positive.
Then, integrating this expression along the flow from $\partial M$ to infinity and using the relation between Hawking energy and the ADM mass we have:
\begin{equation}
m_{ADM}\geq E_H(S_0) + \int_0^\infty \sqrt{\frac{A_t}{16\pi}} \int_{S_t} \left(  \mu +\frac{q}{H}\nu^i j_i\right) dS +\int_0^\infty \sqrt{4\pi} J_t^2 \frac{A_t^{1/2}}{\int_{S_t}\eta ds}dt
\end{equation}

Now because the matter fields have compact support inside $S_T$ the fist integral runs only from $0$ to $T$, then dividing the integral involving the angular momentum from $0$ to $T$ and from $T$ to infinity and using that $J_t=J \quad \forall t\geq T$ we obtain
\begin{multline}
m_{ADM} \geq E_H(S_0) + \int_0^T \sqrt{\frac{A_t}{16\pi}} \int_{S_t} \left( \mu +\frac{q}{H}\nu^i j_i  \right) dS dt \\
+ \sqrt{4\pi} J^2 \int_T^\infty \frac{A_t^{1/2}}{\int_{S_t}  \eta dS}dt + \sqrt{4\pi} \int_0^T J_t^2 \frac{A_t^{1/2}}{\int_{S_t}\eta dS} dt .
\end{multline}
Finally disregarding the last term, using that $E_H(S_0)=\frac{\mathcal{R}_A}{2}$, and writing this expression in terms of $\mathcal{R}_A(S_T)$, $\mathcal{R}(S_T)$ and $m_T$ we obtain \eqref{teo2}.

\end{proof}

\textbf{Remarks}

\vspace{0.5cm}
The notion of size we use, $\mathcal{R}$, albeit apparently artificial at first sight
have proved to be very useful to relate the angular momentum to the total mass in axially symmetric and asymptotically flat initial data \cite{Anglada:2017ryp}. It comes from the particular method we use to relate the angular momentum with the ADM
mass, and gives a good measure of how different the IMCF is from a spherical one.
These kind of measures based on the norm of
the Killing vector have been found to give an appropriate description of size of a
region when describing both regular objects and black holes with angular
momentum \cite{Anglada:2016dbu}, \cite{Reiris:2014tva}, \cite{Reiris:2013jaa},
\cite{Dain:2014}.

Assuming particular properties for the IMCF we can write \eqref{teo1} in terms of the usual measures of size.
The best situation is to have a spherical IMCF, thus $\mathcal{R}^2(\partial M)=\frac{A}{4\pi}$, in which case this proof implies the validity of \eqref{PenJ} in the vaccum case:
\begin{equation}
\label{inqA}
m_{ADM}^2\geq \frac{A}{16 \pi} + \frac{4 \pi J^2}{A}.
\end{equation}
In general we do not expect to have a spherical IMCF in the context we are considering.
For weaker conditions for the IMCF, for example assuming that the flow is convex, we obtain a weaker version of the Penrose inequality with angular momentum for vaccum initial data in terms of the axial radius of the apparent horizon $\mathcal{R}_C=\mathcal{R}_C(\partial M) $:
\begin{equation}
 m_{ADM}^2\geq \frac{A}{16 \pi} + \frac{2}{5}\frac{ J^2}{\mathcal{R}_C^2}.
\end{equation}

\vspace{0.5cm}

If we consider a spherically symmetric initial data then the angular momentum is zero, and the presented result implies the validity of the usual Penrose inequality for apparent horizons in spherical symmetry that was already proved by Malec and O'Murchadha in \cite{Malec:1994}. It is important to mention that in \cite{Ben-Dov:2004} Ben-Dov presents an important counterexample to some formulation of the Penrose inequality in spherical symmetry, but these counterexamples do not contradict either the result of Malec and O'Murchadha or the present result. In \cite{Ben-Dov:2004} the author explicitly constructs an asymptotically flat initial data that contains an apparent horizon and violates the Penrose inequality, but this particular initial data contains past-trapped surfaces outside the apparent horizon, and thus it does not satisfy the hypothesis of theorem \ref{teo1}. From this results one can infer that, if we consider a future apparent horizon $\partial M$, the hypothesis that states that there are no other trapped surfaces in $M$ other than $\partial M$ is necessary to assure the validity of the Penrose inequality.

\vspace{0.5cm}

In the previous work \cite{Anglada:2017ryp} we consider $\partial M$ to be a minimal surface and we use the monotonicity property of Geroch energy. To assure that this monotonicity property is valid we need to have a non-negativity scalar curvature. In order to do this we need to assume some special conditions for the extrinsic curvature, one possible choice is to take $K$ such that for every surface of the flow $q|_{S_t}=0$. This condition, known as the polar gauge condition, together with the assumption that there are no other trapped surfaces in $M$ than $\partial M$, also assures that $\partial M$ is an outermost future apparent horizon. Note that if we assume the condition $q|_{S_t}=0$ then Geroch and Hawking energies are equal along the IMCF, and in this case the previous and present results are the same.
The other possible condition, the most usual, is to consider that the initial data is maximal, $k=0$. With this condition we have a non-negative scalar curvature, and thus Geroch energy is monotonic, but in this case the a minimal surface is a future apparent horizon only if the extrinsic curvature also satisfies $q|_{\partial M}=z|_{\partial M}=0$.
In present work we study the case in which $\partial M$ is a general outermost future apparent horizon, in this case the Geroch energy of $\partial M$ is given by the area of the horizon plus a term that involves the surface integral of $q^2$. Thus Geroch energy is no longer useful to study the problem, we need to use Hawking energy and the maximal condition is not sufficient to assure its monotonicity. The main problem is that one need to control the sing of the second term in \eqref{dtEgralf} to assure that the derivate along the flow of Hawking energy is positive. One possible option, condition \ref{a}, is to control the behavior of $q$ along the IMCF, the other option, condition \ref{b}, is to control the $\nu, \eta$ components of the extrinsic curvature.
Although \ref{a} and \ref{b} are less general than the maximal condition, because both of them depend on a particular foliation of the initial data given by the IMCF, these conditions allows us to extend the previous result to the case in which $q|_{\partial M}\neq 0$. At the present time we do not have a clear geometrical or physical interpretation of these conditions and it will be interesting to study its meaning in detail.

\vspace{0.5cm}

Condition \ref{b} can be fulfilled by choosing a particular form for $s_i$. Note that this condition is a necessary condition to get the monotonicity of Hawking energy if we do not want to assume the very restricted condition \ref{a}. First if we take $s_i$ such that does not have any component on the $\theta$ direction, that is to say $s^i=s\eta^i$ condition \ref{b} can be written in the following way:
\begin{equation}
\begin{split}
g^{ij} \bar \nabla_i s_j =& g^{ij} \eta_j \bar \nabla_i s + s g^{ij}  \bar \nabla_i \eta_j \\
                         =& \eta^i\bar \nabla_i s +  \frac{ s}{2} \left( g^{ij} \bar \nabla_i \eta_j -  g^{ij} \bar \nabla_j \eta_i \right) \\
                         =& \partial_\varphi s \\
\end{split}
\end{equation}
where in the second step we use the Killing equation for $\eta^i$, $\bar \nabla_i \eta_j =- \bar \nabla_j \eta_i$. Then one of the possible choices to get condition \ref{b} is to assume that $s^i=s\eta^i$ and that $s$ does not depend on the coordinate $\varphi$.

Is important to note that if we assume that $s_i=0$, conditions \ref{a} and \ref{b} are not necessary to get the monotonicity of Hawking energy, but in this case we do not have angular momentum. Take an asymptotically flat and axially symmetric initial data that do not have any other trapped surface than $\partial M$, satisfy the DEC and have $s_i=0$. Then for this initial data the existence of a smooth solution of the IMCF is the only necessary condition one needs to prove the positivity of the derivative of Hawking energy along the flow, and thus the only necessary condition one needs to prove the Penrose inequality \eqref{penroseineq}. Then one can infer that the angular momentum generates difficulties in obtaining a foliation of $M$ for which it can be assured that Hawking energy is monotonically increasing.

\vspace{0.5cm}
In respect to the assumption of existence of a smooth solution of the IMCF, the conditions we assume to assure that Hawking energy is monotonic will probably not be fulfilled for a weak flow. Moreover the method we use to relate the angular momentum with the energy strongly depends on having a smooth IMCF.
In order to adapt the present formulation to a weak formulation of the flow one need to control the behavior of all the terms involved in equation on each discontinuity of the flow, and we do not know if it is possible to control the terms involving $q$. In this sense we think the method presented in the previous work \cite{Anglada:2017ryp} has better chances to be adapted to a weak formulation of the IMCF

\section*{Acknowledgments}
I thank Edward Malec and Marc Mars for enlightening discussions.
I would also like to thank Maria E. Gabach-Clement and Omar E. Ortiz for their encouraging support, and I thank the anonymous referees for their helpful and detailed reviews.
This work was partially supported by grants from CONICET and
SECyT, UNC.