# Generalized activity equations for spiking neural network dynamics

Michael A. Buice1 Carson C. Chow2 1Allen Institute for Brain Science
2Laboratory of Biological Modeling, NIDDK, NIH, Bethesda, MD

###### Abstract

Much progress has been made in uncovering the computational capabilities of
spiking neural networks. However, spiking neurons will always be more
expensive to simulate compared to rate neurons because of the inherent
disparity in time scales - the spike duration time is much shorter than the
inter-spike time, which is much shorter than any learning time scale. In
numerical analysis, this is a classic stiff problem. Spiking neurons are also
much more difficult to study analytically. One possible approach to making
spiking networks more tractable is to augment mean field activity models with
some information about spiking correlations. For example, such a generalized
activity model could carry information about spiking rates and correlations
between spikes self-consistently. Here, we will show how this can be
accomplished by constructing a complete formal probabilistic description of
the network and then expanding around a small parameter such as the inverse of
the number of neurons in the network. The mean field theory of the system
gives a rate-like description. The first order terms in the perturbation
expansion keep track of covariances.

## Introduction

Even with the rapid increase in computing power due to Moore’s law and
proposals to simulate the entire human brain notwithstanding Markram and
collaborators (2012), a realistic simulation of a functioning human brain
performing nontrivial tasks is remote. While it is plausible that a network
the size of the human brain could be simulated in real time Izhikevich and
Edelman (2008); Eliasmith et al. (2012) there are no systematic ways to
explore the parameter space. Technology to experimentally determine all the
parameters in a single brain simultaneously does not exist and any attempt to
infer parameters by fitting to data would require exponentially more computing
power than a single simulation. We also have no idea how much detail is
required. Is it sufficient to simulate a large number of single compartment
neurons or do we need multiple-compartments? How much molecular detail is
required? Do we even know all the important biochemical and biophysical
mechanisms? There are an exponential number of ways a simulation would not
work and figuring out which remains computationally intractable. Hence, an
alternative means to provide appropriate prior distributions for parameter
values and model detail is desirable. Current theoretical explorations of the
brain utilize either abstract mean field models or small numbers of more
biophysical spiking models. The regime of large but finite numbers of spiking
neurons remains largely unexplored. It is not fully known what role spike time
correlations play in the brain. It would thus be very useful if mean field
models could be augmented with some spike correlation information.

This paper outlines a scheme to derive generalized activity equations for the
mean and correlation dynamics of a fully deterministic system of coupled
spiking neurons. It synthesizes methods we have developed to solve two
different types of problems. The first problem was how to compute finite
system size effects in a network of coupled oscillators. We adapted the
methods of the kinetic theory of gases and plasmas Ichimaru (1973); Nicholson
(1993) to solve this problem. The method exploits the exchange symmetry of the
oscillators and characterizes the phases of all the oscillators in terms of a
phase density function $\eta(\theta,t)$, where each oscillator is represented
as a point mass in this density. We then write down a formal flux conservation
equation of this density, called the Klimontovich equation, which completely
characterizes the system. However, because the density is not differentiable,
the Klimontovich equation only exists in the weak or distributional sense.
Previously, e.g Desai and Zwanzig (1978); Strogatz and Mirollo (1991); Abbott
and van Vreeswijk (1993); Treves (1993) the equations were made usable by
taking the “mean field limit” of $N\rightarrow\infty$ and assuming that the
density is differentiable in that limit, resulting in what is called the
Vlasov equation. Instead of immediately taking the mean field limit, we
regularize the density by averaging over initial conditions and parameters and
then expand in the inverse system size $N^{-1}$ around the mean field limit.
This results in a system of coupled moment equations known as the BBGKY moment
hierarchy. In Hildebrand et al. (2007), we solved the moment equations for the
Kuramoto model perturbatively to compute the pair correlation function between
oscillators. However, the procedure was somewhat ad hoc and complicated. We
then subsequently showed in Buice and Chow (2007), that the BBGKY moment
hierarchy could be recast in terms of a density functional of the phase
density. This density functional could be written down explicitly as an
integral over all possible phase histories, i.e. a Feynman-Kac path integral.
The advantage of using this density functional formalism is that the moments
to arbitrary order in $1/N$ could be computed as a steepest-descent expansion
of the path integral, which can be expressed in terms of Feynman diagrams.
This made the calculation more systematic and mechanical. We later applied the
same formalism to synaptically coupled spiking models Buice and Chow (2013a).

Concurrently with this line of research, we also explored the question of how
to generalize population activity equations, such as the Wilson-Cowan
equations, to include the effects of correlations. The motivation for this
question is that the Wilson-Cowan equations are mean field equations and do
not capture the effects of spike-time correlations. For example, the gain in
the Wilson-Cowan equations is fixed, (which is a valid approximation when the
neurons fire asynchronously), but correlations in the firing times can change
the gain Salinas and Sejnowski (2000). Thus, it would be useful to develop a
systematic procedure to augment population activity equations to include spike
correlation effects. The approach we took was to posit plausible microscopic
stochastic dynamics, dubbed the spike model, that reduced to the Wilson-Cowan
equations in the mean field limit and compute the self-consistent moment
equations from that microscopic theory. Buice and Cowan Buice and Cowan (2009)
showed that the solution of the master equation of the spike model could be
expressed formally in terms of a path integral over all possible spiking
histories. The random variable in the path integral is a spike count whereas
in the path integral for the deterministic phase model we described above, the
random variable is a phase density. To generate a system of moment equations
for the microscopic stochastic system, we transformed the random spike count
variable in the path integral into moment variables Buice et al. (2010). This
is accomplished using the effective action approach of field theory, where the
exponent of the cumulant generating functional, called the action, which is a
function of the random variable is Legendre transformed into an effective
action of the cumulants. The desired generalized Wilson-Cowan activity
equations are then the equations of motion of the effective action. This is
analogous to the transformation from Lagrangian variables of position and
velocity to Hamiltonian variables of position and momentum. Here, we show how
to apply the effective action approach to a deterministic system of
synaptically coupled spiking neurons to derive a set of moment equations.

## Approach

Consider a network of single compartment conductance-based neurons

$\displaystyle C\frac{dV_{i}}{dt}$ $\displaystyle=$
$\displaystyle-\sum_{r=1}^{n}g_{r}(x_{i}^{r})(V_{i}-v_{r})+\sum_{j=1}^{N}g_{ij}s_{j}(t)$
$\displaystyle\tau_{i}^{r}\frac{dx_{i}^{r}}{dt}$ $\displaystyle=$
$\displaystyle f(V_{i},x_{i})$ $\displaystyle\tau_{j}\frac{ds_{j}}{dt}$
$\displaystyle=$ $\displaystyle h(V_{j},s_{j})$
$\displaystyle\tau_{g}\frac{dg_{ij}}{dt}$ $\displaystyle=$
$\displaystyle\phi(g_{ij},V)$

The equations are remarkably stiff with time scales spanning orders of
magnitude from milliseconds for ion channels, to seconds for adaptation, and
from hours to years for changes in synaptic weights and connections. Parameter
values must be assigned for $10^{11}$ neurons with $10^{4}$ connections each.
Here, we present a formalism to derive a set of reduced activity equations
directly from a network of deterministic spiking neurons that capture the
spike rate and spike correlation dynamics. The formalism first constructs a
density functional for the firing dynamics of all the neurons in a network. It
then systematically marginalizes the unwanted degrees of freedom to isolate a
set of self-consistent equations for the desired quantities. For heuristic
reasons, we derive an example set of generalized activity equations for the
first and second cumulants of the firing dynamics of a simple spiking model
but the method can be applied to any spiking model.

A convenient form to express spiking dynamics is with a phase oscillator.
Consider the quadratic integrate-and-fire neuron

$\frac{dV_{i}}{dt}=I_{i}+V_{i}^{2}+\alpha_{i}u(t)$ (1)

where $I$ is an external current and $u(t)$ are the synaptic currents with
some weight $\alpha_{i}$. The spike is said to occur when $V$ goes to infinity
whereupon it is reset to minus infinity. The quadratic nonlinearity ensures
that this transit will occur in a finite amount of time. The substitution
$V=\tan(\theta/2)$ yields the theta model Ermentrout and Kopell (1986):

$\frac{d\theta_{i}}{dt}=1-\cos\theta_{i}+(1+\cos\theta_{i})(I_{i}+\alpha_{i}u)$
(2)

which is the normal form of a Type I neuron near the bifurcation to firing
Ermentrout (1996). The phase neuron is an adequate approximation to spiking
dynamics provided the inputs are not overly strong as to disturb the limit
cycle. The phase neuron also includes realistic dynamics such as not firing
when the input is below threshold. Coupled phase models arise naturally in
weakly coupled neural networks Ermentrout and Kopell (1991); Golomb and Hansel
(2000); Hoppensteadt and Izhikevich (1997). They include the Kuramoto model
Kuramoto (1984), which we have previously analyzed Hildebrand et al. (2007);
Buice and Chow (2007).

Here, we consider the phase dynamics of a set of $N$ coupled phase neurons
obeying

$\displaystyle\dot{\theta}_{i}$ $\displaystyle=$ $\displaystyle
F(\theta,\gamma_{i},u(t))$ (3) $\displaystyle\dot{u}(t)$ $\displaystyle=$
$\displaystyle-\beta u(t)+\beta\nu(t)$ (4) $\displaystyle\nu(t)$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{j=1}^{N}\sum_{l}\delta(t-t^{l}_{j})$ (5)

where each neuron has a phase $\theta_{i}$ that is indexed by $i$, $u$ is a
global synaptic drive, $F(\theta,\gamma,u)$ is the phase and synaptic drive
dependent frequency, $\gamma_{i}$ represents all the parameters for neuron $i$
drawn from a distribution with density $g(\gamma)$, $\nu$ is the population
firing rate of the network,$t^{l}_{j}$ is the $l$th firing time of neuron $j$
and a neuron fires when its phase crosses $\pi$. In the present paper, we
consider all-to-all or global coupling through a synaptic drive variable
$u(t)$. However, our basic approach is not restricted to global coupling.

We can encapsulate the phase information of all the neurons into a neuron
density function Hildebrand et al. (2007); Buice and Chow (2007, 2011, 2013b,
2013a).

$\eta(\theta,\gamma,t)=\frac{1}{N}\sum_{i=1}^{N}\delta(\theta-\theta_{i}(t))\delta(\gamma-\gamma_{i})$
(6)

where $\delta(\cdot)$ is the Dirac delta functional, and $\theta_{i}(t)$ is a
solution to system (3)-(5). The neuron density gives a count of the number of
neurons with phase $\theta$ and synaptic strength $\gamma$ at time $t$. Using
the fact that the Dirac delta functional in (5) can be expressed as
$\sum_{l}\delta(t-t_{j}^{l})=\dot{\theta_{j}}\delta(\pi-\theta_{j})$, the
population firing rate can be rewritten as

$\nu(t)=\int d\gamma\,F(\pi,\gamma,u(t))\eta(\pi,\gamma,t)$ (7)

The neuron density formally obeys the conservation equation

$\frac{\partial}{\partial
t}\eta(\theta,\gamma,t)+\frac{\partial}{\partial\theta}\left[F\eta(\theta,\gamma,t)\right]=0$
(8)

with initial condition $\eta(\theta,\gamma,t_{0})=\eta_{0}(\theta,\gamma)$ and
$u(t_{0})=u_{0}$. Equation (8) is known as the Klimontovich equation Ichimaru
(1973); Liboff (2003). The Klimontovich equation, the equation for the
synaptic drive (4), and the firing rate expressed in terms of the neuron
density (7), fully define the system. The system is still fully deterministic
but is now in a form where various sets of reduced descriptions can be
derived. Here, we will produce an example of a set of reduced equations or
generalized activity equations that capture some aspects of the spiking
dynamics. The path we take towards the end will require the introduction of
some formal machinery that may obscure the intuition around the
approximations. However, we feel that it is useful because it provides a
systematic and controlled way of generating averaged quantities that can be
easily generalized.

For finite $N$, (8) is only valid in the weak or distributional sense since
$\eta$ is not differentiable. In the $N\rightarrow\infty$ limit, it has been
argued that $\eta$ will approach a smooth density $\rho$ that evolves
according to the Vlasov equation that has the same form as (8) but with $\eta$
replaced by $\rho$ Ichimaru (1973); Nicholson (1993); Desai and Zwanzig
(1978); Strogatz and Mirollo (1991); Hildebrand et al. (2007). This has been
proved rigorously in the case where noise is added using the theory of coupled
diffusions McKean Jr (1966); Faugeras et al. (2009); Touboul (2012); Baladron
et al. (2012). This $N\rightarrow\infty$ limit is called mean field theory. In
mean field theory, the original microscopic many body neuronal network is
represented by a smooth macroscopic density function. In other words, the
ensemble of networks prepared with different microscopic initial conditions is
sharply peaked at the mean field solution. For large but finite $N$, there
will be deviations away from mean field Hildebrand et al. (2007); Buice and
Chow (2007, 2013a, 2013b). These deviations can be characterized in terms of a
distribution over an ensemble of coupled networks that are all prepared with
different initial conditions and parameter values. Here, we show how a
perturbation theory in $N^{-1}$ can be developed to expand around the mean
field solution. This requires the construction of the probability density
functional over the ensemble of spiking neural networks. We adapt the tools of
statistical field theory to perform such a construction.

### Formalism

The complete description of the system given by equations (4), (7), and (8)
can be written as

$\displaystyle\dot{u}(t)+\beta u(t)-\beta\int
d\gamma\,F(\pi,\gamma,u(t))\eta(\pi,\gamma,t)=0$ (9)
$\displaystyle\frac{\partial}{\partial
t}\eta(\theta,\gamma,t)+\frac{\partial}{\partial\theta}\left[F(\theta,\gamma,u(t))\eta(\theta,\gamma,t)\right]\equiv{\cal
L}\eta=0$ (10)

The probability density functional governing the system specified by the
synaptic drive and Klimontovich equations (9) and (10) given initial
conditions $(\eta_{0},u_{0})$ can be written as

$\displaystyle P[\eta,u]=$ $\displaystyle\int{\cal D}u_{0}(t)\,{\cal
D}\eta_{0}(\theta,\gamma)\,P[\eta,u|\eta_{0},u_{0}]\,P_{0}[\eta_{0},u_{0},\gamma]$
(11)

where $P[\eta,u|\eta_{0},u_{0}]$ is the conditional probability density
functional of the functions $(\eta,u)$, and $P_{0}[\eta_{0},u_{0}]$ is the
density functional over initial conditions of the system. The integral is a
Feynman-Kac path integral over all allowed initial condition functions.
Formally we can write $P[\eta,u|\eta_{0},u_{0}]$ as a point mass (Dirac delta)
located at the solutions of (9) and (10) given the initial conditions:

$\displaystyle\delta\left[{\cal
L}\eta-\eta_{0}\delta(t-t_{0})\right]\delta\left[\dot{u}+\beta u-\beta\int
d\gamma\,F(\pi,\gamma,u(t))\eta(\pi,\gamma,t)-u_{0}\delta(t-t_{0})\right]$

The probability density functional (11) is then

$\displaystyle P[\eta,u]=$ $\displaystyle\int{\cal D}u_{0}(t)\,{\cal
D}\eta_{0}(\theta,\gamma)\,\delta\left[{\cal
L}\eta-\eta_{0}\delta(t-t_{0})\right]$
$\displaystyle\times\delta\left[\dot{u}+\beta u-\beta\int
d\gamma\,F(\pi,\gamma,u(t))\eta(\pi,\gamma,t)-u_{0}\delta(t-t_{0})\right]P_{0}[\eta_{0},u_{0},\gamma]$
(12)

Equation (12) can be made useful by noting that the Fourier representation of
a Dirac delta is given by $\delta(x)\propto\int dk\,e^{ikx}$. Using the
infinite dimensional Fourier functional transform then gives

$P[\eta,u]=\int{\cal D}\tilde{\eta}{\cal
D}\tilde{u}\,e^{-NS[\eta,\tilde{\eta},u,\tilde{u}]}.$

The exponent $S[\eta,u]$ in the probability density functional is called the
action and has the form

$S=S_{u}+S_{\varphi}+S_{0}$ (13)

where

$\displaystyle S_{\varphi}=\int d\theta d\gamma
dt\,\tilde{\varphi}(x)\left[\partial_{t}\varphi(x)+\partial_{\theta}F(\theta,\gamma,u(t))\varphi(x)\right]$
(14)

represents the contribution of the transformed neuron density to the action,

$\displaystyle S_{u}=\frac{1}{N}\int dt\,\tilde{u}(t)\left(\dot{u}(t)+\beta
u(t)-\beta\int d\gamma
F(\pi,\gamma,u(t))[\tilde{\varphi}(\pi,\gamma,t)+1]\varphi(\pi,\gamma,t)\right)$
(15)

represents the global synaptic drive,
$S_{0}[\tilde{\varphi}_{0}(x_{0}),u_{0}(t_{0})]$ represents the initial
conditions, and $x=(\theta,\gamma,t)$. For the case where the neurons are
considered to be independent in the initial state, we have

$\displaystyle S_{0}[\tilde{\varphi}_{0}(x_{0}),u_{0}(t_{0})]$
$\displaystyle=-\frac{1}{N}\tilde{u}(t_{0})u_{0}-\ln\left(1+\int d\theta
d\gamma\tilde{\varphi}_{0}(\theta,\gamma,t_{0})\rho_{0}(\theta,\gamma,t_{0})\right)$
(16)

where $u_{0}$ is the initial value of the coupling variable and
$\rho_{0}(\theta,\gamma,t)$ is the distribution from which the initial
configuration is drawn for each neuron. The action includes two imaginary
auxiliary response fields (indicated with a tilde), which are the infinite
dimensional Fourier transform variables. The factor of $1/N$ appears to ensure
correct scaling between the $u$ and $\varphi$ variables since $u$ applies to a
single neuron while $\varphi$ applies to the entire population. The full
derivation is given in Buice and Chow (2013a) and a review of path integral
methods applied to differential equations is given in Buice and Chow (2010).
In the course of the derivation we have made a Doi-Peliti-Jannsen
transformation Janssen and Täuber (2005); Buice and Chow (2013a), given by

$\displaystyle\varphi(x)$ $\displaystyle=$
$\displaystyle\eta(x)e^{-\tilde{\eta}(x)}$ $\displaystyle\tilde{\varphi}(x)$
$\displaystyle=$ $\displaystyle e^{\tilde{\eta}(x)}-1$

In deriving the action, we have explicitly chosen the Ito convention so that
the auxiliary variables only depend on variables in the past. The action (13)
contains all the information about the statistics of the network.

The moments for this distribution can be obtained by taking functional
derivatives of a moment generating functional. Generally, the moment
generating function for a random variable is given by the expectation value of
the exponential of that variable with a single parameter. Because our goal is
to transform to new variables for the first and second cumulants, we form a
“two-field” moment generating functional, which includes a second parameter
for pairs of random variables,

$\displaystyle\exp(N\,$ $\displaystyle W[J,K])=$ $\displaystyle\int{\cal
D}\xi\,\exp\left[-NS[\xi]+N\int dx\,J^{i}(x)\xi_{i}(x)+\frac{N}{2}\int
dxdx^{\prime}\xi_{i}(x)K^{ij}(x,x^{\prime})\xi_{j}(x^{\prime})\right]$ (17)

where $J$ and $K$ are moment generating fields, $\xi_{1}(x)=u(t)$,
$\xi_{2}(x)=\tilde{u}(t)$, $\xi_{3}(x)=\varphi(x)$,
$\xi_{4}(x)=\tilde{\varphi}(x)$, and $x=(\theta,\gamma,t)$. Einstein summation
convention is observed beween upper and lower indices. Unindexed variables
represent vectors. The integration measure $dx$ is assumed to be $dt$ when
involving indices 1 and 2. Covariances between an odd and even index
corresponds to a covariance between a field and an auxiliary field. Based on
the structure of the action $S$ and (17) we see that this represents a linear
propagator and by causality and the choice of the Ito convention is only
nonzero if the time of the auxiliary field is evaluated at an earlier time
than the field. Covariances between two even indices correspond to that
between two auxiliary fields and are always zero because of the Ito
convention.

The mean and covariances of $\xi$ can be obtained by taking derivatives of the
action $W[J,K]$ in (17), with respect to $J$ and $K$ and setting $J$ and $K$
to zero:

$\displaystyle\frac{\delta W}{\delta J^{i}}$
$\displaystyle=\left.\langle\xi_{i}\rangle\right|_{J,K=0}$
$\displaystyle\frac{\delta W}{\delta K^{ij}}$
$\displaystyle=\left.\frac{1}{2}\langle\xi_{i}\xi_{j}\rangle\right|_{J,K=0}$

Expressions for these moments can be computed by expanding the path integral
in (17) perturbatively around some mean field solution. However, this can be
unwieldy if closed form expressions for the mean field equations do not exist.
Alternatively, the moments at any order can be expressed as self-consistent
dynamical equations that can be analyzed or simulated numerically. Such
equations form a set of generalized activity equations for the means
$a_{i}=\langle\xi_{i}\rangle$, and covariances
$C_{ij}=N[\langle\xi_{i}\xi_{j}\rangle-a_{i}a_{j}]$.

We derive the generalized activity equations by Legendre transforming the
action $W$, which is a function of $J$ and $K$, to an effective action
$\Gamma$ that is a function of $a$ and $C$. Just as a Fourier transform
expresses a function in terms of its frequencies, a Legendre transform
expresses a convex function in terms of its derivatives. This is appropriate
for our case because the moments are derivatives of the action. The Legendre
transform of $W[J,K]$ is

$\Gamma[a,C]=-W[J,K]+\int dxJ^{i}a_{i}+\frac{1}{2}\int
dxdx^{\prime}\left[a_{i}a_{j}+\frac{1}{N}C_{ij}\right]K^{ij}$ (18)

which must obey the constraints

$\displaystyle\frac{\delta W}{\delta J^{i}}$ $\displaystyle=a_{i}$
$\displaystyle\frac{\delta W}{\delta K^{ij}}$
$\displaystyle=\frac{1}{2}\left[a_{i}a_{j}+\frac{1}{N}C_{ij}\right]$

and

$\displaystyle\frac{\delta\Gamma}{\delta a_{i}}$
$\displaystyle\equiv\Gamma^{i,00}=J^{i}+\frac{1}{2}a_{j}\left[K^{ij}+K^{ji}\right]$
$\displaystyle\frac{\delta\Gamma}{\delta C_{ij}}$
$\displaystyle\equiv\Gamma^{0,ij}=\frac{1}{2N}K^{ij}$ (19)

The generalized activity equations are given by the equations of motion of the
effective action, in direct analogy to the Euler-Lagrange equations of
classical mechanics, and are obtained by setting $J^{i}=0$ and $K^{ij}=0$ in
(19).

In essence, what the effective action does is to take a probabilistic
(statistical mechanical) system in the variables $\xi$ with action $S$ and
transform them to a deterministic (classical mechanical) system with an action
$\Gamma$. Our approach here follows that used in Buice et al. (2010) to
construct generalized activity equations for the Wilson Cowan model. However,
there are major differences between that system and this one. In Buice et al.
(2010), the microscopic equations were for the spike counts of an inherently
probabilistic model so the effective action and ensuing generalized activity
equations could be constructed directly from the Markovian spike count
dynamics. Here, we start from deterministically firing individual neurons and
get to a probabilistic description through the Klimontovich equation. It would
be straightforward to include stochastic effects into the spiking dynamics.

Using (18) in (17) gives

$\displaystyle\exp(-N\,\Gamma[a,C])=\int{\cal D}\psi\,\exp\left[-NS[\xi]+N\int
dx\,J^{i}(\xi_{i}-a_{i})\right.$ $\displaystyle+\left.\frac{N}{2}\int
dxdx^{\prime}\left[\xi_{i}\xi_{j}-a_{i}a_{j}-\frac{1}{N}C_{ij}\right]K^{ij}\right]$
(20)

where $J$ and $K$ are constrainted by (19). We cannot compute the effective
action explicitly but we can compute it perturbatively in $N^{-1}$. We first
perform a shift $\xi_{i}=a_{i}+\psi_{i}$, expand the action as
$S[a+\psi]=S[a]+\int dx(L^{i}[a]\psi_{i}+(1/2)\int
dx^{\prime}L^{ij}[a]\psi_{i}\psi_{j})+\cdots$ and substitute for $J$ and $K$
with the constraints (19) to obtain

$\displaystyle\exp(-N\,\Gamma[a,C])=$
$\displaystyle\exp(-NS[a]-N\operatorname{Tr}\Gamma^{0,ij}C_{ij})\int{\cal
D}\psi\,\exp\left[-N\int dx\bigg{(}L^{i}[a]\psi_{i}\right.$
$\displaystyle\left.+\frac{1}{2}\int
dx^{\prime}L^{ij}[a]\psi_{i}\psi_{j}\bigg{)}+N\int
dx\,\Gamma^{i,00}\psi_{i}+N^{2}\int
dxdx^{\prime}\psi_{i}\psi_{j}\Gamma^{0,ij}\right]$ (21)

where

$\operatorname{Tr}A^{ij}B_{ij}=\int
dxdx^{\prime}A^{ij}(x,x^{\prime})B_{ij}(x,x^{\prime})$ (22)

Our goal is to construct an expansion for $\Gamma$ by collecting terms in
successive orders of $N^{-1}$ in the path integral of (21). Expanding $\Gamma$
as $\Gamma[a,C]=\Gamma_{0}+N^{-1}\Gamma_{1}+N^{-2}\Gamma_{2}$ and equating
coefficients of $N$ in (21) immediately leads to the conclusion that
$\Gamma_{0}=S[a]$, which gives

$\displaystyle\exp(-N\,\Gamma[a,C])=\exp\left(-NS[a]-\operatorname{Tr}\Gamma_{1}^{0,ij}C_{ij}\right)\int{\cal
D}\psi\,\exp\left[-\frac{N}{2}\int dxL^{ij}[a]\psi_{i}\psi_{j}\right.$
$\displaystyle+\left.N\int dx\,\Gamma_{1}^{0,ij}\psi_{i}\psi_{j}\right]$

where higher order terms in $N^{-1}$ are not included. To lowest nonzero order
$\Gamma^{0,ij}=N^{-1}\Gamma_{1}^{0,ij}$ since $\Gamma_{0}$ is only a function
of $a$ and not $C$. If we set

$\Gamma_{1}^{0,ij}=(1/2)L^{ij}-(1/2)Q^{ij},$ (23)

we obtain

$\displaystyle\exp(-N\,\Gamma[a,C])=\exp\left(-NS[a]-\frac{1}{2}\operatorname{Tr}L^{ij}C_{ij}+\frac{1}{2}\operatorname{Tr}Q^{ij}C_{ij}\right)$
$\displaystyle\times\int{\cal D}\psi\,\exp\left[-\frac{N}{2}\int
dx\,Q^{ij}[a]\psi_{i}\psi_{j}\right]$ (24)

to order $N^{-1}$. $Q^{ij}$ is an unknown function of $a$ and $C$, which we
will deduce using self-consistency. The path integral in (24), which is an
infinite dimensional Gaussian that can be explicitly integrated, is
proportional to $1/\sqrt{\det Q^{ij}}=\exp(-(1/2)\ln\det
Q^{ij})=\exp(-(1/2)\operatorname{Tr}\ln Q^{ij})$, using properties of
matrices. Hence, (24) becomes

$\displaystyle\exp(-N\,\Gamma[a,C])=\exp\left(-NS[a]-\frac{1}{2}\operatorname{Tr}L^{ij}C_{ij}-\frac{1}{2}\operatorname{Tr}Q^{ij}C_{ij}+\frac{1}{2}\operatorname{Tr}\ln
Q^{ij}\right)$

and

$\displaystyle\Gamma[a,C]=S[a]+\frac{1}{2N}\operatorname{Tr}L^{ij}C_{ij}+\frac{1}{2N}\operatorname{Tr}\ln
Q^{ij}-\frac{1}{2N}\operatorname{Tr}Q^{ij}C_{ij}$

Taking the derivative of $\Gamma$ with respect to $C_{ij}$ yields

$\displaystyle\Gamma^{0,ij}=\frac{1}{2N}\left(L^{ij}+(Q^{-1})^{kl}\frac{\partial}{\partial
C_{ij}}Q^{lk}-\frac{\partial}{\partial C_{ij}}(Q^{kl}C_{lk})\right)$

Self consistency with (23) then requires that $Q^{ij}=(C^{-1})^{ij}$ which
leads to the effective action

$\Gamma[a,C]=S[a]+\frac{1}{2N}\operatorname{Tr}\ln(C^{-1})^{ij}+\frac{1}{2N}\operatorname{Tr}L^{ij}C_{ij}$
(25)

where

$\int
dx^{\prime}\,(C^{-1})^{ik}(x,x^{\prime})C_{kj}(x^{\prime},x_{0})=\delta_{ij}\delta(x-x_{0})$

and we have dropped the irrelevant constant terms.

The equations of motion to order $N^{-1}$ are obtained from (19) with $J^{i}$
and $K^{ij}$ set to zero:

$\displaystyle\frac{\delta S[a]}{\delta
a_{i}}+\frac{1}{2N}\frac{\delta}{\delta a_{i}}\operatorname{Tr}L^{ij}C_{ij}=0$
(26) $\displaystyle\frac{1}{2N}[-(C^{-1})^{ij}+L^{ij}]=0$ (27)

and (27) can be rewritten as

$\displaystyle\int
dx^{\prime}L^{ik}(x,x^{\prime})C_{kj}(x^{\prime},x_{0})=\delta_{ij}\delta(x-x_{0})$
(28)

Hence, given any network of spiking neurons, we can write down the a set of
generalized activity equations for the mean and covariance functions by 1)
constructing a neuron density function, 2) writing down the conservation law
(Klimontovich equation), 3) constructing the action and 4) using formulas (26)
and (28). We could have constructed these equations directly by multiplying
the Klimontovich and synaptic drive equations by various factors of $u$ and
$\eta$ and recombining. However, as we saw in Buice et al. (2010) this is not
a straightforward calculation. The effective action approach makes this much
more systematic and mechanical.

### Phase Model Example

We now present a simple example to demonstrate the concepts and approximations
involved in our expansion. Our goal is not to analyze the system per se but
only to demonstrate the application of our method in a heuristic setting. We
begin with a simple nonleaky integrate-and-fire neuron model, which responds
to a global coupling variable. This is a special case of the dynamics given
above, with $F$ given by

$\displaystyle F[\theta,\gamma,u]=I(t)+\gamma u$ (29)

The action from (14) and (15) is

$\displaystyle S[a]=$ $\displaystyle\int d\theta d\gamma
dt\,a_{4}(x)\left[\partial_{t}a_{3}(x)+\partial_{\theta}(I+\gamma
a_{1}(t))a_{3}(x)\right]$ $\displaystyle+\frac{1}{N}\int
dt\,a_{2}(t)\left(\dot{a}_{1}(t)+\beta a_{1}(t)-\beta\int d\gamma\,(I+\gamma
a_{1}(t))[a_{4}(\pi,\gamma,t)+1]a_{3}(\pi,\gamma,t)\right)$ (30)

and we ignore initial conditions for now.

In order to construct the generalized activity equations we need to compute
the first and second derivatives of the action $L^{i}$ and $L^{ij}$. Taking
the first derivative of (30) gives

$\displaystyle L^{1}[a](x,x^{\prime})=\frac{\delta S[a(x)]}{\delta
a_{1}(t^{\prime})}$ $\displaystyle=\int d\theta d\gamma\,dt\gamma
a_{4}(x)\partial_{\theta}a_{3}(x)\delta(t-t^{\prime})$
$\displaystyle+\frac{1}{N}\left[\int
dt\,a_{2}(t)\frac{d}{dt}\delta(t-t^{\prime})+\beta
a_{2}(t^{\prime})-a_{2}(t^{\prime})\beta\int
d\gamma\,\gamma[a_{4}(\pi,\gamma,t^{\prime})+1]a_{3}(\pi,\gamma,t^{\prime})\right]$
$\displaystyle L^{2}[a](x,x^{\prime})=\frac{\delta S[a(x)]}{\delta
a_{2}(t^{\prime})}$
$\displaystyle=\frac{1}{N}\left[\frac{da_{1}}{dt^{\prime}}+\beta
a_{1}(t^{\prime})-\beta\int d\gamma(I+\gamma
a_{1}(t^{\prime}))[a_{4}(\pi,\gamma,t^{\prime})+1]a_{3}(\pi,\gamma,t^{\prime})\right]$
$\displaystyle L^{3}[a](x,x^{\prime})=\frac{\delta S[a(x)]}{\delta
a_{3}(x^{\prime})}$ $\displaystyle=\int
dt\,a_{4}(\theta^{\prime},\gamma^{\prime},t)\partial_{t}\delta(t-t^{\prime})+\int
d\theta
a_{4}(\theta,\gamma^{\prime},t^{\prime})\partial_{\theta}(I+\gamma^{\prime}a_{1}(t^{\prime}))\delta(\theta-\theta^{\prime})$
$\displaystyle-\frac{\beta}{N}a_{2}(t^{\prime})(I+\gamma^{\prime}a_{1}(t^{\prime}))(a_{4}(\pi,\gamma^{\prime},t^{\prime})+1)\delta(\pi-\theta^{\prime})$
$\displaystyle L^{4}[a](x,x^{\prime})=\frac{\delta S[a(x)]}{\delta
a_{4}(x^{\prime})}$
$\displaystyle=\partial_{t^{\prime}}a_{3}(x^{\prime})+\partial_{\theta^{\prime}}(I+\gamma^{\prime}a_{1}(t^{\prime}))a_{3}(x^{\prime})-\frac{\beta}{N}a_{2}(t^{\prime})(I+\gamma^{\prime}a_{1}(t^{\prime}))a_{3}(\pi,\gamma^{\prime},t^{\prime})\delta(\pi-\theta^{\prime})$

The mean field equations are obtained by solving $L^{i}=0$ using
(LABEL:firstderivs). We immediately see that $a_{2}=a_{4}=0$ are solutions,
which leaves us with

$\displaystyle\dot{a_{1}}+\beta a_{1}-\beta\int d\gamma(I+\gamma
a_{1})a_{3}(\pi,\gamma,t)=0$ (32) $\displaystyle\partial_{t}a_{3}+(I+\gamma
a_{1})\partial_{\theta}a_{3}=0$ (33)

The mean field equations should be compared to those of the spike response
model Gerstner (1995, 2000). We can also solve (33) directly to obtain

$a_{3}(x,t)=\rho_{0}\left(\theta-\int_{t_{0}}^{t}dt^{\prime}\left[I_{\Omega}(t^{\prime})+\gamma
a_{1}(t^{\prime})\right],\gamma,\Omega\right)$

where $\rho_{0}$ is the initial distribution. If the neurons are initially
distributed uniformly in phase, then $\rho_{0}=g(\gamma)/2\pi$ and the mean
field equations reduce to

$\dot{a_{1}}(t)+\beta
a_{1}(t)-\frac{\beta}{2\pi}\left(I+\bar{\gamma}a_{1}(t)\right)=0$ (34)

which has the form of the Wilson-Cowan equation, with
$(\beta/2\pi)\left(I+\bar{\gamma}a_{1}\right)$ acting as a gain function.
Hence, the Wilson-Cowan equation is a full description of the infinitely large
system limit of a network of globally coupled simple phase oscillators in the
asynchronous state. For all other initial conditions, the one-neuron
conservation equation (called the Vlasov equation in kinetic theory) must be
included in mean field theory.

To go beyond mean field theory we need to compute
$L^{ij}(x,x^{\prime},x^{\prime\prime})=\delta L^{i}(x,x^{\prime})/\delta
a_{j}(x^{\prime\prime})$:

$\displaystyle L^{11}[a]$ $\displaystyle=0$ $\displaystyle L^{12}[a]$
$\displaystyle=\frac{1}{N}\left[-\frac{d}{dt^{\prime\prime}}+\beta-\beta\int
d\gamma\,\gamma[a_{4}(\pi,\gamma,t^{\prime\prime})+1]a_{3}(\pi,\gamma,t^{\prime\prime})]\right]\delta(t^{\prime\prime}-t^{\prime})$
$\displaystyle L^{13}[a]$ $\displaystyle=\left[\gamma^{\prime\prime}\int
d\theta\,a_{4}(x)\delta(\gamma-\gamma^{\prime\prime})\partial_{\theta}\delta(\theta-\theta^{\prime\prime})-\frac{\beta}{N}\gamma^{\prime\prime}a_{2}(t^{\prime})[a_{4}(\pi,\gamma^{\prime\prime},t^{\prime\prime})+1]\delta(\pi-\theta^{\prime\prime})\right]\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{14}[a]$
$\displaystyle=\left[\gamma^{\prime\prime}\partial_{\theta^{\prime\prime}}a_{3}(x^{\prime\prime})-\frac{\beta}{N}\gamma^{\prime\prime}a_{2}(t^{\prime\prime})a_{3}(\pi,\gamma^{\prime\prime},t^{\prime\prime})\delta(\pi-\theta^{\prime\prime})\right]\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{21}[a]$
$\displaystyle=\frac{1}{N}\left[\frac{d}{dt^{\prime}}+\beta-\beta\int
d\gamma\,\gamma[a_{4}(\pi,\gamma,t^{\prime})+1]a_{3}(\pi,\gamma,t^{\prime})\right]\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{22}[a]$ $\displaystyle=0$ $\displaystyle L^{23}[a]$
$\displaystyle=-\frac{\beta}{N}(I+\gamma^{\prime\prime}a_{1}(t^{\prime}))[a_{4}(\pi,\gamma^{\prime\prime},t^{\prime}))+1]\delta(\pi-\theta^{\prime})\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{24}[a]$
$\displaystyle=-\frac{\beta}{N}(I+\gamma^{\prime\prime}a_{1}(t^{\prime}))a_{3}(\pi,\gamma^{\prime\prime},t^{\prime})]\delta(\pi-\theta^{\prime\prime})\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{31}[a]$ $\displaystyle=\left[\int
d\theta\,a_{4}(\theta,\gamma^{\prime},t^{\prime})\gamma^{\prime}\partial_{\theta}\delta(\theta-\theta^{\prime})-\frac{\beta}{N}a_{2}(t^{\prime})\gamma^{\prime}[a_{4}(\pi,\gamma^{\prime},t^{\prime})+1]\delta(\pi-\theta^{\prime})\right]\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{32}[a]$
$\displaystyle=-\frac{\beta}{N}(I+\gamma^{\prime}a_{1}(t^{\prime}))(a_{4}(\pi,\gamma^{\prime},t^{\prime})+1)\delta(\pi-\theta^{\prime})\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{33}[a]$ $\displaystyle=0$ $\displaystyle L^{34}[a]$
$\displaystyle=\left[\delta(\theta^{\prime}-\theta^{\prime\prime})\partial_{t^{\prime\prime}}-\partial_{\theta^{\prime\prime}}(I+\gamma^{\prime}a_{1}(t^{\prime}))\delta(\theta^{\prime\prime}-\theta^{\prime})\right.$
$\displaystyle\left.-\frac{\beta}{N}a_{2}(t^{\prime})(I+\gamma^{\prime}a_{1}(t^{\prime}))\delta(\pi-\theta^{\prime})\delta(\pi-\theta^{\prime\prime})\right]\delta(\gamma^{\prime}-\gamma^{\prime\prime})\delta(t^{\prime\prime}-t^{\prime})$
$\displaystyle L^{41}[a]$
$\displaystyle=\left[\partial_{\theta^{\prime}}\gamma^{\prime}a_{3}(x^{\prime})-\frac{\beta}{N}a_{2}(t^{\prime})\gamma^{\prime}a_{3}(\pi,\gamma^{\prime},t^{\prime})\delta(\pi-\theta^{\prime})\right]\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{42}[a]$
$\displaystyle=-\frac{\beta}{N}(I+\gamma^{\prime}a_{1}(t^{\prime}))a_{3}(\pi,\gamma^{\prime},t^{\prime})\delta(\pi-\theta^{\prime})\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{43}[a]$
$\displaystyle=\partial_{t^{\prime}}\delta(x^{\prime}-x^{\prime\prime})+\partial_{\theta^{\prime}}(I+\gamma^{\prime}a_{1}(t^{\prime}))\delta(x^{\prime}-x^{\prime\prime})$
$\displaystyle-\frac{\beta}{N}a_{2}(t^{\prime})(I+\gamma
a_{1}(t^{\prime}))\delta(\pi-\theta^{\prime})\delta(\pi-\theta^{\prime\prime})\delta(\gamma^{\prime}-\gamma^{\prime\prime})\delta(t^{\prime}-t^{\prime\prime})$
$\displaystyle L^{44}[a]$ $\displaystyle=0$

The activity equations for the means to order $N^{-1}$ are given by (26). The
only nonzero contributions are given by $L^{13}$ and $L^{31}$ resulting in

$\displaystyle L^{2}+\frac{1}{2N}\frac{\delta}{\delta a_{2}}\int
dxdx^{\prime}(L^{13}C_{13}+L^{31}C_{31})=0$ $\displaystyle
L^{4}+\frac{1}{2N}\frac{\delta}{\delta a_{4}}\int
dxdx^{\prime}(L^{13}C_{13}+L^{31}C_{31})=0$

since $a_{2}=a_{4}=0$ and correlations involving response variables (even
indices) will be zero for equal times. The full activity equations for the
means are thus

$\displaystyle\dot{a_{1}}+\beta a_{1}-\beta\int d\gamma(I+\gamma
a_{1})a_{3}(\pi,\gamma,t)-\frac{\beta}{N}\int d\gamma\,\gamma
C(\pi,\gamma,t)=0$ (35) $\displaystyle\partial_{t}a_{3}+(I+\gamma
a_{1})\partial_{\theta}a_{3}+\frac{1}{N}\gamma\partial_{\theta}C(\theta,\gamma,t)=0$
(36)

where
$C(\theta,\gamma,t)=C_{13}(t;\theta,\gamma,t)=C_{31}(\theta,\gamma,t;t)$.

We can now use the $L^{ij}$ in (28) to obtain activity equations for $C_{ij}$.
There will be sixteen coupled equations in total but the applicable nonzero
ones are

$\displaystyle\left[\frac{d}{dt}+\beta-\beta\int d\gamma\,\gamma
a_{3}(\pi,\gamma,t)]\right]$ $\displaystyle C_{11}(t;t_{0})-\beta\int
d\gamma\,(I+\gamma a_{1})C_{31}(\pi,\gamma,t;t_{0})$ $\displaystyle-\beta\int
d\gamma\,(I+\gamma a_{1}(t))a_{3}(\pi,\gamma,t)C_{41}(\pi,\gamma,t;t_{0})=0$
(37) $\displaystyle\left[\frac{d}{dt}+\beta-\beta\int d\gamma\,\gamma
a_{3}(\pi,\gamma,t)\right]$ $\displaystyle C_{13}(t;x_{0})-\beta\int
d\gamma\,(I+\gamma a_{1})C_{33}(\pi,\gamma,t;x_{0})$ $\displaystyle-\beta\int
d\gamma\,(I+\gamma a_{1}(t))a_{3}(\pi,\gamma,t)C_{43}(\pi,\gamma,t;x_{0})=0$
(38)
$\displaystyle\gamma\partial_{\theta}a_{3}(x)C_{11}(t;t_{0})+[\partial_{t}+$
$\displaystyle(I+\gamma a_{1})\partial_{\theta}]C_{31}(x;t_{0})$
$\displaystyle-\frac{\beta}{N}(I+\gamma
a_{1}(t))a_{3}(\pi,\gamma,t)\delta(\pi-\theta)C_{21}(t,t_{0})=0$ (39)
$\displaystyle\gamma\partial_{\theta}a_{3}(x)C_{13}(t;x_{0})+[\partial_{t}+$
$\displaystyle(I+\gamma a_{1}(t))\partial_{\theta}]C_{33}(x,x_{0})$
$\displaystyle-\frac{\beta}{N}(I+\gamma
a_{1}(t))a_{3}(\pi,\gamma,t)\delta(\pi-\theta)C_{23}(t,x_{0})=0$ (40)

Adding (38) and (39) and taking the limit $t_{0}\rightarrow t$ and setting
$\theta_{0}=\theta$, $\gamma_{0}=\gamma$ gives

$\displaystyle\partial_{t}$ $\displaystyle
C(\theta,\gamma,t)+\left[\beta-\beta\int
d\gamma^{\prime}\,\gamma^{\prime}a_{3}(\pi,\gamma^{\prime},t)+(I+\gamma
a_{1})\partial_{\theta}\right]C(\theta,\gamma,t)-\beta\int
d\gamma^{\prime}\,(I+\gamma^{\prime}a_{1})C_{33}(\pi,\gamma^{\prime},t;x)$
$\displaystyle-2\beta(I+\gamma
a_{1}(t))a_{3}(\pi,\gamma,t)\delta(\pi-\theta)+\gamma\partial_{\theta}a_{3}(x)C_{11}(t;t)=0$

where we use the fact that $C_{21}(t,t^{\prime})=N$ and
$C_{43}(x;x^{\prime})=\delta(\theta-\theta^{\prime})\delta(\gamma-\gamma^{\prime})$
in the limit of $t^{\prime}$ approaching $t$ from below and equal to zero when
approaching from above. Adding (37) and (40) to themselves with $t$ and
$t_{0}$ interchanged and taking the limit of $t_{0}$ approaching $t$ gives

$\displaystyle\left[\frac{d}{dt}+2\beta-2\beta\int d\gamma\,\gamma
a_{3}(\pi,\gamma,t)]\right]C_{11}(t;t)-2\beta\int d\gamma\,(I+\gamma
a_{1})C(\pi,\gamma,t)=0$ $\displaystyle\left[\partial_{t}+(I+\gamma
a_{1}(t))\partial_{\theta}\right]C_{33}(x;x)+2\gamma[\partial_{\theta}a_{3}(x)]C(x)=0$

because $C_{41}(x;t)=0$ and $C_{23}(t;x)=0$. Putting this all together, we get
the generalized activity equations

$\displaystyle\frac{{da_{1}}}{dt}+\beta a_{1}(t)-\beta\int d\gamma(I+\gamma
a_{1}(t))a_{3}(\pi,\gamma,t)-\frac{\beta}{N}\int d\gamma\,\gamma
C(\pi,\gamma,t)=0$ (41)
$\displaystyle\partial_{t}a_{3}(\theta,\gamma,t)+(I+\gamma
a_{1})\partial_{\theta}a_{3}(\theta,\gamma,t)+\frac{1}{N}\gamma\partial_{\theta}C(\theta,\gamma,t)=0$
(42) $\displaystyle\partial_{t}C(\theta,\gamma,t)+\left[\beta-\beta\int
d\gamma^{\prime}\,\gamma^{\prime}a_{3}(\pi,\gamma^{\prime},t)+(I+\gamma
a_{1})\partial_{\theta}\right]C(\theta,\gamma,t)$ $\displaystyle-\beta\int
d\gamma^{\prime}\,(I+\gamma^{\prime}a_{1})C_{33}(\pi,\gamma^{\prime},t;\theta,\gamma,t)-2\beta(I+\gamma
a_{1}(t))a_{3}(\theta,\gamma,t)\delta(\pi-\theta)$
$\displaystyle+\gamma\partial_{\theta}a_{3}(\theta,\gamma,t)C_{11}(t;t)=0$
(43) $\displaystyle\left[\frac{d}{dt}+2\beta-2\beta\int d\gamma\,\gamma
a_{3}(\pi,\gamma,t)]\right]C_{11}(t;t)-2\beta\int d\gamma\,(I+\gamma
a_{1})C(\pi,\gamma,t)=0$ (44) $\displaystyle\left[\partial_{t}+(I+\gamma
a_{1}(t))\partial_{\theta}\right]C_{33}(\theta,\gamma,t;\theta,\gamma,t)+2\gamma\partial_{\theta}a_{3}(\theta,\gamma,t)C(\theta,\gamma,t)=0$
(45)

Initial conditions, which are specified in the action, are required for each
of these equations. The derivation of these equations using classical means
require careful consideration for each particular model. Our method provides a
blanket mechanistic algorithm. We propose that these equations represent a new
scheme for studying neural networks.

Equations (41)-(45) are the complete self-consistent generalized activity
equations for the mean and correlations to order $N^{-1}$. It is a system of
partial differential equations in $t$ and $\theta$. These equations can be
directly analyzed or numerically simulated. Although the equations seem
complicated, one must bear in mind that they represent the dynamics of the
system averaged over initial conditions and unknown parameters. Hence, the
solution of this PDE system replaces multiple simulations of the original
system. In previous work, we required over a million simulations of the
original system to obtained adequate statistics Buice and Chow (2013a). There
is also a possibility that simplifying approximations can be applied to such
systems. The system has complete phase memory because the original system was
fully deterministic. However, the inclusion of stochastic effects will shorten
the memory and possibly simplify the dynamics. It will pose no problem to
include such stochastic effects. In fact, the formalism is actually more
suited for stochastic systems Buice et al. (2010).

## Discussion

The main goal of this paper was to show how to systematically derive
generalized activity equations for the ensemble averaged moments of a
deterministically coupled network of spiking neurons. Our method utilizes a
path integral formalism that makes the process algorithmic. The resulting
equations could be derived using more conventional perturbative methods
although possibly with more calculational difficulty as we found before Buice
et al. (2010). For example, for the case of the stochastic spike model, Buice
et al. Buice et al. (2010) presumed that the Wilson-Cowan activity variable
was the rate of a Poisson process and derived a system of generalized activity
equations that corresponded to deviations around Poisson firing. Bressloff
Bressloff (2010), on the other hand, assumed that the Wilson-Cowan activity
variable was a mean density and used a system-size expansion to derive an
alternative set of generalized activity equations for the spike model. The
classical derivations of these two interpretations look quite different and
the differences and similarities between them are not readily apparent.
However, the connections between the two types of expansions are very
transparent using the path integral formalism.

Here, we derived equations for the rate and covariances (first and second
cumulants) of a deterministic synaptically coupled spiking network as a system
size expansion to first order. However, our method is not restricted to these
choices. What is particularly advantageous about the path integral formalism
is that it is straightforward to generalize to include higher order cumulants,
extend to higher orders in the inverse system size, or to expand in other
small parameters such as the inverse of a slow time scale. The action fully
specifies the system and all questions regarding the system can be addressed
with it.

To give a concrete illustration of the method, we derived the self-consistent
generalized activity equations for the rates and covariances to order $N^{-1}$
for a simple phase model. The resulting equations consist of ordinary and
partial differential equations. This is to be expected since the original
system was fully deterministic and memory cannot be lost. Even mean field
theory requires the solution of an advective partial differential equation.
The properties of these and similar equations remain to be explored
computationally and analytically. The system is possibly simpler near the
asynchronous state, which is marginally stable in mean field theory like the
Kuramoto model Strogatz and Mirollo (1991) and like the Kuramoto model, we
conjecture that the finite size effects will stabilize the asynchronous state
Hildebrand et al. (2007); Buice and Chow (2007). The addition of noise will
also stabilize the asynchronous state. Near asynchrony could be exploited to
generate simplified versions of the asynchronous state.

We considered a globally connected network, which allowed us to assume that
networks for different parameter values and initial conditions converge
towards a “typical” system in the large $N$ limit. However, this property may
not hold for more realistic networks. While the formalism describing the
ensemble average will hold regardless of this assumption, the utility of the
equations as descriptions of a particular network behavior may suffer. For
example, heterogeneity in the connectivity (as opposed to the global
connectivity we consider here) may threaten this assumption. This is the case
with so called “chaotic random networks” Sompolinsky et al. (1988) in which
there is a spin-glass transition owing to the variance of the connectivity
crossing a critical threshold. This results in the loss of a “typical” system
in the large $N$ limit requiring an effective stochastic equation which
incorporates the noise induced by the network heterogeneity. Whether the
expansion we present here is useful without further consideration depends upon
whether the network heterogeneity induces this sort of effect. This is an area
for future work. A simpler issue arises when there are a small discrete number
of “typical” systems (such as with bistable solutions to the continuity
equation). In this case, there are noise induced transitions between states.
While the formalism has a means of computing this transition Elgart and
Kamenev (2004), we do not consider this case here.

An alternative means to incorporate heterogeneous connections is to consider a
network of coupled systems. In such a network, a set of generalized activity
equations, such as those derived here or simplified versions, would be derived
for each local system, together with equations governing the covariances
between the local systems. Correlation based learning dynamics could then be
imposed on the connections between the local systems. Such a network could
serve as a generalization of current rate based neural networks to include the
effects of spike correlations with applications to both neuroscience and
machine learning.

## Acknowledgments

This work was supported by the Intramural Research Program of the NIH, NIDDK.

## References

  * Markram and collaborators (2012) H. Markram and collaborators, _The Human Brain Project: a report to the European commission_ (The HBP-PS Consortium, Lausanne, 2012).
  * Izhikevich and Edelman (2008) E. M. Izhikevich and G. M. Edelman, Proceedings of the national Academy, USA 105, 3593 (2008).
  * Eliasmith et al. (2012) C. Eliasmith, T. C. Stewart, X. Choo, T. Bekolay, T. DeWolf, Y. Tang, and D. Rasmussen, Science 338, 1202 (2012).
  * Ichimaru (1973) S. Ichimaru, _Basic principles of plasma physics, a statistical approach._ (Benjamin, New York, 1973).
  * Nicholson (1993) D. R. Nicholson, _Introduction to plasma theory_ (Krieger Publishing Co, Malabar, FL, 1993).
  * Desai and Zwanzig (1978) R. Desai and R. Zwanzig, Journal of Statistical Physics 19, 1 (1978).
  * Strogatz and Mirollo (1991) S. Strogatz and R. Mirollo, Journal of Statistical Physics 63, 613 (1991).
  * Abbott and van Vreeswijk (1993) L. Abbott and C. van Vreeswijk, Physical Review E 48, 1483 (1993).
  * Treves (1993) A. Treves, Network: Computation in Neural Systems 4, 259 (1993).
  * Hildebrand et al. (2007) E. J. Hildebrand, M. A. Buice, and C. C. Chow, Physical Review Letters 98, 054101 (2007).
  * Buice and Chow (2007) M. A. Buice and C. C. Chow, Physical Review E 76, 031118 (2007).
  * Buice and Chow (2013a) M. A. Buice and C. C. Chow, PLoS Computational Biology 9, e1002872 (2013a).
  * Salinas and Sejnowski (2000) E. Salinas and T. J. Sejnowski, The Journal of neuroscience : the official journal of the Society for Neuroscience 20, 6193 (2000).
  * Buice and Cowan (2009) M. Buice and J. Cowan, Progress in Biophysics and Molecular Biology 99, 53 (2009).
  * Buice et al. (2010) M. A. Buice, J. D. Cowan, and C. C. Chow, Neural Computation 22, 377 (2010).
  * Ermentrout and Kopell (1986) G. Ermentrout and N. Kopell, SIAM Journal on Applied Mathematics pp. 233–253 (1986).
  * Ermentrout (1996) B. Ermentrout, Neural computation 8, 979 (1996).
  * Ermentrout and Kopell (1991) G. Ermentrout and N. Kopell, Journal of Mathematical Biology 29, 195 (1991).
  * Golomb and Hansel (2000) D. Golomb and D. Hansel, Neural computation 12, 1095 (2000).
  * Hoppensteadt and Izhikevich (1997) F. Hoppensteadt and E. Izhikevich, _Weakly Connected Neural Networks_ (Springer-Verlag, New York, 1997).
  * Kuramoto (1984) Y. Kuramoto, _Chemical Oscillations, Waves, and Turbulence_ (Springer, Berlin, 1984).
  * Buice and Chow (2011) M. A. Buice and C. C. Chow, Physical Review E 84, 051120 (2011).
  * Buice and Chow (2013b) M. A. Buice and C. C. Chow, Journal of Statistical Mechanics: Theory and Experiment 2013, P03003 (2013b).
  * Liboff (2003) R. L. Liboff, _Kinetic Theory_ (Springer, New York, 2003).
  * McKean Jr (1966) H. McKean Jr, Proceedings of the National Academy of Sciences of the United States of America 56, 1907 (1966).
  * Faugeras et al. (2009) O. Faugeras, J. Touboul, and B. Cessac, Frontiers in computational neuroscience 3, 1 (2009).
  * Touboul (2012) J. Touboul, Physica D: Nonlinear Phenomena 241, 1223 (2012).
  * Baladron et al. (2012) J. Baladron, D. Fasoli, O. , and J. Touboul, J. Math. Neurosci. 2, 10 (2012).
  * Buice and Chow (2010) M. A. Buice and C. C. Chow, arXiv.org p. 1009.5966 (2010).
  * Janssen and Täuber (2005) H. Janssen and U. Täuber, Annals of Physics 315, 147 (2005).
  * Gerstner (1995) W. Gerstner, Physical Review E 51, 738 (1995).
  * Gerstner (2000) W. Gerstner, Neural computation 12, 43 (2000).
  * Bressloff (2010) P. C. Bressloff, SIAM Journal on Applied Mathematics 70, 1488 (2010).
  * Sompolinsky et al. (1988) H. Sompolinsky, A. Crisanti, and H. Sommers, Physical Review Letters 61, 259 (1988).
  * Elgart and Kamenev (2004) V. Elgart and A. Kamenev, arXiv cond-mat.stat-mech (2004).


# New spinorial mass-quasilocal angular momentum inequality for initial data
with marginally outer trapped surface

Jarosław Kopiński  Center for Quantum Mathematics and Physics (QMAP),
Department of Mathematics, University of California, Davis, CA95616, USA
<EMAIL_ADDRESS>,  Alberto Soria  Escuela Técnica Superior de Ingenieros
de Telecomunicación (ETSIT), Departamento de Matemática Aplicada a las TIC,
Universidad Politécnica de Madrid, Avenida Complutense 30, 28040, Madrid,
Spain<EMAIL_ADDRESS>and  Juan A. Valiente Kroon  School of
Mathematical Sciences, Queen Mary University of London, E1 4NS, London, UK
<EMAIL_ADDRESS>
###### Abstract.

We prove a new geometric inequality that relates the Arnowitt-Deser-Misner
(ADM) mass of initial data to a quasilocal angular momentum of a marginally
outer trapped surface (MOTS) inner boundary. The inequality is expressed in
terms of a 1-spinor, which satisfies an intrinsic first-order Dirac-type
equation. Furthermore, we show that if the initial data is axisymmetric, then
the divergence-free vector used to define the quasilocal angular momentum
cannot be a Killing field of the generic boundary.

##  1\. Introduction

Geometric inequalities arise naturally in General Relativity (GR) as relations
involving quantities characterizing black holes, like mass, angular momentum,
and horizon area. One of the most significant example of a geometric
inequality is the positive (ADM) mass theorem. It was first proven by Schoen
and Yau in [21] for manifolds of dimension $d\leq 7$, and later by Witten with
the use of spinor formalism for arbitraty $d$ [27]. The latter proof was
extended to the case of initial data with trapped surfaces by Gibbons,
Hawking, Horowitz, and Perry in [10]. The spinorial approach was further
adapted by Ludvigsen and Vickers in the context of the Bondi mass in [17]. A
refined version of the positivity of the ADM mass has been formulated by
Penrose in the form of a lower bound on this quantity in terms of the horizon
area of a black hole. If true, it would provide further evidence in favor of
the weak cosmic censorship conjecture [19]. There exists a stronger version of
the Penrose inequality involving the angular momentum of the initial data,
namely

(1) $m\geq\left(\frac{|S|}{16\pi}+\frac{4\pi
J^{2}}{|S|}\right)^{\frac{1}{2}},$

where $m$, $J$ and $|S|$ are the ADM mass, angular momentum and the outermost
apparent horizon area respectively—see e.g. [18, 9, 8] for more details. It
admits a rigidity case, where equality exclusively occurs for the Kerr black
hole. A quasilocal version of this relation states that

$m\geq\left(\frac{|S|}{16\pi}+\frac{4\pi
J^{2}_{BH}}{|S|}\right)^{\frac{1}{2}},$

where $J_{BH}$ is the quasilocal angular momentum of the horizon.

Geometric inequalities for black holes remain a very active area of research
with new interesting results being obtained. Among them is a bound on the ADM
energy in terms of horizon area, angular momentum, and charge obtained by
Jaracz and Khuri in [11]. In [2, 3] a different approach was considered by
Anglada, who used the monotonic properties of the Geroch and Hawking energy
along the inverse mean curvature flow in order to prove a Penrose-like
inequality with angular momentum. The first author and Tafel considered
perturbations of Schwarzschild data and showed that (1) holds in this setting
[14, 15]. Another refinement to the Penrose inequality with angular momentum
has been proven by Alaee and Kunduri for 4-dimensional biaxially symmetric
maximal initial data [1]. Additionally, recent numerical results such as the
ones obtained in [16] give support to the validity of (1) in its full
generality. The examples presented above are far from exhaustive, providing a
glimpse into contemporary research in geometric inequalities in GR. We refer
the reader to [18, 9] for further references.

In the present work a spinorial approach is used to obtain a geometric
inequality involving the ADM mass of the initial data and a quasilocal angular
momentum (à la Szabados [24, 25]) of the MOTS inner boundary. It generalises
the result presented in [13] to the case of non-vanishing connection 1-form on
the normal bundle of the boundary. The solvability of the boundary value
problem for the so-called _approximate twistor equation_ is still an essential
ingredient for deriving the main result. The existence of solution is used to
obtain a basic mass inequality

(2) $4\pi
m\geq\sqrt{2}\oint_{\partial\mathcal{S}}\widehat{\phi}^{A}\gamma_{A}{}^{B}\not{\mathcal{D}}_{BC}\phi^{C}dS,$

where $\not{\mathcal{D}}_{AB}$ and $\gamma_{AB}$ are the 2-dimensional Sen
connection and the complex metric on the boundary respectively (see below for
details), while $\phi_{A}$ is a valence 1 spinor on $\partial\mathcal{S}$. The
right-hand side of (2) can be rewritten in terms of the inner null expansion
$\theta^{-}$ of the boundary and the aforementioned angular momentum, provided
that $\phi_{A}$ satisfies a certain first-order Dirac-type equation.
Ultimately,

(3) $4\pi
m\geq\sqrt{2}\oint\displaylimits_{\mathbb{S}^{2}}\rho^{\prime}|\widetilde{\phi}_{0}|^{2}\Omega
d\mathbb{S}^{2}+\frac{\kappa}{\sqrt{2}}O\left[\widetilde{\phi},U\right],$

where $\rho^{\prime}=-\tfrac{\theta^{-}}{2}$, $\widetilde{\phi}_{A}$ is a
Dirac eigenspinor on $\mathbb{S}^{2}$, $\Omega$ is a conformal factor relating
the metrics of $\partial\mathcal{S}$ and $\mathbb{S}^{2}$ and
$O\left[\widetilde{\phi},U\right]$ a quasilocal angular momentum depending on
$\widetilde{\phi}_{A}$ and a rotation potential $U$ defined below. It should
be noted that the integrals in (3) are now taken with respect to the 2-sphere
volume element.

A natural symmetry associated with the angular momentum is the existence of
axial Killing vector. Therefore, with such assumption we analyze a scenario
where the quasilocal angular momentum is generated by such vector (on top of
arising from a spinor $\phi_{A}$) and show that it is in fact impossible for a
generic MOTS inner boundary $\partial\mathcal{S}$.

The article is structured as follows: Section 2 provides a discussion of our
main mathematical tools, in particular a new formalism for the $1+1+2$
decomposition of spinors. Section 3 is an adaptation of the result of [13] to
the case of non-vanishing connection 1-form on the normal bundle of
$\partial\mathcal{S}$. In Section 4 we present the main result of this work, a
new mass-quasilocal angular momentum inequality for the initial data with a
MOTS. In the last section we particularise our analysis to the axisymmetric
setting and show that the divergence-free vector generating the quasilocal
angular momentum cannot arise simultaneously from a first-order Dirac-type
equation and be a Killing vector of the boundary.

In the following, 4-dimensional metrics are considered to have the signature
$(+---)$. As a result, Riemannian 3- and 2-dimensional metrics will be
negative definite. Whenever appropriate, we will expand spinorial expressions
using either the Geroch-Held-Penrose (GHP) or Newman-Penrose (NP) formalism,
following the conventions outlined in [20]. Throughout this paper, we employ
abstract index notation, with lowercase letters representing tensorial indices
and uppercase letters representing spinorial indices. Bold font will be used
to denote components in a basis.

##  2\. Preliminaries

###  2.1. Basic setting

An initial data set $(\mathcal{S},h_{ab},K_{ab})$ for the vacuum Einstein
field equations is said to be _asymptotically Schwarzschildean_ if the metric
$h_{ab}$ and the second fundamental form $K_{ab}$ satisfy the decay conditions

(4a) $\displaystyle
h_{ab}=-\left(1+\frac{2m}{r}\right)\delta_{ab}+o_{\infty}(r^{-3/2}),$ (4b)
$\displaystyle K_{ab}=o_{\infty}(r^{-5/2}),$

with $r^{2}\equiv(x^{1})^{2}+(x^{2})^{2}+(x^{3})^{2}$,
$(x^{\alpha})=(x^{1},x^{2},x^{3})$ being asymptotically Cartesian coordinates
and $m$ the ADM mass. In this work we assume that $\mathcal{S}$ has an inner
boundary $\partial\mathcal{S}$ which is a topological 2-sphere and is equipped
with a metric $\sigma_{ab}$. We consider a $1+1+2$ spinor formalism, first
proposed in [23] by Szabados, and based on the use of $SL(2,\mathbb{C})$
spinors. Maintaining the same philosophy as in [13], the so-called
$SU(2,\mathbb{C})$ spinors (or space spinors) introduced in [22] will be
essential to our purposes since they allow to work efficiently on spacelike
hypersurfaces. For more information on the spinor formalism, we refer the
reader to [26, 4].

Let $\tau^{AA^{\prime}}$ be a spinorial counterpart of the orthogonal future
vector $\tau^{a}$ to $\mathcal{S}$ such that
$\tau_{AA^{\prime}}\tau^{AA^{\prime}}=2$. Likewise, we will denote a spinorial
counterpart of the normal vector $\rho^{a}$ to $\partial\mathcal{S}$ on
$\mathcal{S}$ as $\rho^{AA^{\prime}}$ and assume that
$\rho_{AA^{\prime}}\rho^{AA^{\prime}}=-2$. The spinors $\tau_{AA^{\prime}}$
and $\rho_{AA^{\prime}}$ are orthogonal —i.e.
$\tau_{AA^{\prime}}\rho^{AA^{\prime}}=0$. We consider dyads
$\\{o^{A},\iota^{A}\\}$ such that

$\displaystyle\tau_{AA^{\prime}}=o_{A}\overline{o}_{A^{\prime}}+\iota_{A}\overline{\iota}_{A^{\prime}},$
$\displaystyle\rho_{AA^{\prime}}=o_{A}\overline{o}_{A^{\prime}}-\iota_{A}\overline{\iota}_{A^{\prime}}.$

The spinor $\tau_{AA^{\prime}}$ is used to construct a space-spinor version of
a given spinor. In particular,

$\gamma_{AB}\equiv\tau_{(B}{}^{A^{\prime}}\rho_{A)A^{\prime}}$

is the space-spinor version of $\rho_{AA^{\prime}}$, also called the _complex
metric_. It satisfies $\gamma_{A}{}^{B}\gamma_{B}{}^{C}=\delta_{A}{}^{C}$ and
can be expressed as

$\gamma_{AB}=o_{A}\iota_{B}+o_{B}\iota_{A}$

with the use of spin dyad.

The spinorial counterpart of the projection operator $\Pi_{a}{}^{b}$ onto the
2-dimensional surface $\partial\mathcal{S}$ can now be defined as

(5)
$\begin{split}\Pi_{AA^{\prime}}{}^{BB^{\prime}}&\equiv\delta_{A}{}^{B}\delta_{A^{\prime}}{}^{B^{\prime}}-\tfrac{1}{2}\tau_{AA^{\prime}}\tau^{BB^{\prime}}+\tfrac{1}{2}\rho_{AA^{\prime}}\rho^{BB^{\prime}}=\tfrac{1}{2}\left(\delta_{A}{}^{B}\delta_{A^{\prime}}{}^{B^{\prime}}-\gamma_{A}{}^{B}\overline{\gamma}_{A^{\prime}}{}^{B^{\prime}}\right).\end{split}$

Similarly, the spinorial counterpart of the projector
$T_{AA^{\prime}BB^{\prime}}$ onto $\mathcal{S}$ reads

$T_{AA^{\prime}}{}^{BB^{\prime}}\equiv\tfrac{1}{2}\left(\delta_{A}{}^{B}\delta_{A^{\prime}}{}^{B^{\prime}}-\tfrac{1}{2}\tau_{AA^{\prime}}\tau^{BB^{\prime}}\right).$

Let $\nabla_{AA^{\prime}}$ be the spinorial counterpart of the spacetime
covariant derivative $\nabla_{a}$. The $T_{AA^{\prime}}{}^{BB^{\prime}}$
projector allows to define the $3$-dimensional Sen connection
$\mathcal{D}_{AA^{\prime}}$ associated to $\nabla_{AA^{\prime}}$ as

$\mathcal{D}_{AA^{\prime}}\pi_{C}\equiv
T_{AA^{\prime}}{}^{BB^{\prime}}\nabla_{BB^{\prime}}\pi_{C}.$

As mentioned above, the $SU(2,\mathbb{C})$ (i.e. space-spinor version) of
$\mathcal{D}_{AA^{\prime}}$ can be constructed by means of
$\tau^{AA^{\prime}}$ as

$\mathcal{D}_{AB}=\tau_{(B}{}^{A^{\prime}}\mathcal{D}_{A)A^{\prime}}.$

The space-spinor version of the 3-dimensional Levi-Civita connection on
$\mathcal{S}$ can be recovered form $\mathcal{D}_{AB}$ via

$\nabla_{AB}\pi_{C}=\mathcal{D}_{AB}\pi_{C}-\tfrac{1}{2}K_{ABC}{}^{Q}\pi_{Q},$

where
$K_{ABCD}\equiv\tau_{D}{}^{C^{\prime}}\mathcal{D}_{AB}\tau_{CC^{\prime}}$ is
the Weingarten spinor. It decomposes as

$K_{ABCD}=\Omega_{ABCD}-\tfrac{1}{3}K\epsilon_{A(C}\epsilon_{D)B},$

where $\Omega_{ABCD}\equiv K_{(ABCD)}$ is its fully symmetrized part, and
$K\equiv K_{AB}{}^{AB}$ is the mean curvature of $\mathcal{S}$. The
3-dimensional Levi-Civita operator satisfies $\nabla_{AB}\epsilon_{CD}=0$.

Given a spinor $\pi_{A_{1}{}\dots A_{K}}$, its Hermitian conjugate is defined
as follows,

$\widehat{\pi}_{A_{1}{}\dots
A_{K}}\equiv\tau_{A_{1}}{}^{A_{1}^{\prime}}\dots\tau_{A_{k}}{}^{A_{k}^{\prime}}\,\overline{\pi}_{A_{1}^{\prime}\dots
A_{k}^{\prime}}.$

A spinor $\pi_{A_{1}{}\dots A_{K}}$ is said to be real if

$\widehat{\pi}_{A_{1}B_{1}{}\dots A_{k}B_{k}}{}^{C_{1}D_{1}{}\dots
C_{m}D_{m}}=(-1)^{(k+m)}\pi_{A_{1}B_{1}{}\dots
A_{k}B_{k}}{}^{C_{1}D_{1}{}\dots C_{m}D_{m}}.$

The space counterpart of the Levi-Civita connection $\nabla_{AB}$ is real in
the sense that $\widehat{\nabla_{AB}\pi_{C}}=-\nabla_{AB}\widehat{\pi}_{C}$,
while

$\widehat{\mathcal{D}_{AB}\pi_{C}}=-\mathcal{D}_{AB}\widehat{\pi}_{C}+K_{ABC}{}^{D}\widehat{\pi}_{D}.$

###  2.2. On the inner boundary.

A 2-dimensional Sen connection $\not{\mathcal{D}}_{AA^{\prime}}$ on
$\partial\mathcal{S}$ arises as a $\Pi$-projection of $\nabla_{AA^{\prime}}$,
i.e.

(6)
$\not{\mathcal{D}}_{AA^{\prime}}\equiv\Pi_{AA^{\prime}}{}^{BB^{\prime}}\nabla_{BB^{\prime}},$

and its associated $SU(2,\mathbb{C})$ version is given by
$\not{\mathcal{D}}_{AB}\equiv\tau_{(B}{}^{A^{\prime}}\not{\mathcal{D}}_{A)A^{\prime}}$.
It can be promoted to the 2-dimensional Levi-Civita connection
$\not{\nabla}_{AA^{\prime}}$ with the use of the transition spinor
$Q_{AA^{\prime}BC}$,

(7)
$\not{\nabla}_{AA^{\prime}}v_{BB^{\prime}}=\not{\mathcal{D}}_{AA^{\prime}}v_{BB^{\prime}}-Q_{AA^{\prime}B}{}^{C}v_{CB^{\prime}}-\overline{Q}_{AA^{\prime}B^{\prime}}{}^{C^{\prime}}v_{BC^{\prime}},$

where

(8)
$Q_{AA^{\prime}BC}\equiv-\tfrac{1}{2}\gamma_{C}{}^{D}\not{\mathcal{D}}_{AA^{\prime}}\gamma_{BD}.$

The $\not{\nabla}_{AA^{\prime}}$ connection is torsion-free by definition,
i.e.
$\left(\not{\nabla}_{AA^{\prime}}\not{\nabla}_{BB^{\prime}}-\not{\nabla}_{BB^{\prime}}\not{\nabla}_{AA^{\prime}}\right)\phi=0$,
and its curvature spinor
$\not{r}_{CC^{\prime}DD^{\prime}AA^{\prime}BB^{\prime}}$ can be defined
defined with the use of the following relation

(9)
$\left(\not{\nabla}_{AA^{\prime}}\not{\nabla}_{BB^{\prime}}-\not{\nabla}_{BB^{\prime}}\not{\nabla}_{AA^{\prime}}\right)\pi^{C}=\not{r}^{C}{}_{QAA^{\prime}BB^{\prime}}\pi^{Q}=\left(m_{a}\overline{m}_{b}-\overline{m}_{a}m_{b}\right)\left(\rho\rho^{\prime}-\sigma\sigma^{\prime}+\Psi_{2}\right)\gamma^{CD}\pi_{D},$

where $m^{AA^{\prime}}\equiv o^{A}\overline{\iota}^{A^{\prime}}$,
$\overline{m}^{AA^{\prime}}\equiv\iota^{A}\overline{o}^{A^{\prime}}$, and
$\Psi_{2}\equiv\Psi_{ABCD}o^{A}o^{B}\iota^{C}\iota^{D}$ is a component of the
Weyl spinor $\Psi_{ABCD}$.

Another 2-dimensional connection ($\not{D}_{AB}$) can be obtained by
considering a space-spinor counterpart of $\not{\nabla}_{AA^{\prime}}$, i.e.

(10) $\not{D}_{AB}\equiv\tau_{(B}{}^{A^{\prime}}\not{\nabla}_{A)A^{\prime}}.$

It is particularly useful in some calculations and can be related to
$\not{\mathcal{D}}_{AB}$ via

(11)
$\not{D}_{AB}\pi_{C}=\not{\mathcal{D}}_{AB}\pi_{C}-Q_{AB}{}^{Q}{}_{C}\pi_{Q},$

where the transition spinor is now given by

(12)
$Q_{AB}{}^{C}{}_{D}\equiv-\tfrac{1}{2}\gamma_{D}{}^{Q}\not{\mathcal{D}}_{AB}\gamma_{Q}{}^{C}=\sigma^{\prime}o_{A}o_{B}o_{C}o_{D}+\sigma\iota_{A}\iota_{B}\iota_{C}\iota_{D}-\rho
o_{A}o_{B}\iota_{C}\iota_{D}-\rho^{\prime}\iota_{A}\iota_{B}o_{C}o_{D}.$

A natural choice for the ingoing and outgoing null vectors $k^{a}$ and $l^{a}$
spanning the normal bundle to $\partial\mathcal{S}$ is given by

$k^{a}=\tfrac{1}{2}(\tau^{a}-\rho^{a}),\qquad
l^{a}=\tfrac{1}{2}(\tau^{a}+\rho^{a}).$

The nature of a trapped surfaces is determined by the causal character and
orientation of its mean curvature vector, or equivalently, by the signs of the
associated inner and outer null expansions,

$\theta^{-}=\sigma^{ab}\nabla_{a}k_{b},\qquad\theta^{+}=\sigma^{ab}\nabla_{a}l_{b}.$

Making use of Proposition $4.14.2$ in [20] we can express $\theta^{-}$ and
$\theta^{+}$ in terms of a GHP spin coefficients,

$\theta^{-}=-2\rho^{\prime},\qquad\theta^{+}=-2\rho.$

We are now ready to define a MOTS.

######  Definition 1.

The boundary $\partial\mathcal{S}$ is said to be a marginally outer trapped
surface (MOTS) if $\theta^{+}=0$ and $\theta^{-}\leq 0$, i.e. if $\rho=0$ and
$\rho^{\prime}\geq 0$ on $\partial\mathcal{S}$.

The 2-dimensional connection $\not{D}_{AB}$ annihilates $\epsilon_{AB}$ and
$\gamma_{AB}$:

$\not{D}_{AB}\epsilon_{CD}=0,\quad\not{D}_{AB}\gamma_{CD}=0.$

However, $\not{D}_{AB}$ is not a Levi-Civita connection on
$\partial\mathcal{S}$ as it has a non-vanishing torsion,

(13)
$\begin{split}\not{D}_{AB}\not{D}_{CD}\phi-\not{D}_{CD}\not{D}_{AB}\phi&=\tfrac{1}{2}\left(A_{AB}\gamma_{C}{}^{X}\delta_{D}{}^{Y}-A_{CD}\gamma_{A}{}^{X}\delta_{B}{}^{Y}\right)\not{D}_{XY}\phi,\end{split}$

where $A_{AB}$ is defined as

(14)
$A_{AB}\equiv\tau^{CC^{\prime}}\not{\mathcal{D}}_{AB}\rho_{CC^{\prime}}=2\left(\alpha+\overline{\beta}\right)o_{A}o_{B}-2\left(\overline{\alpha}+\beta\right)\iota_{A}\iota_{B}.$

This spinor is real, $\widehat{A}_{AB}=-A_{AB}$, and satisfies
$\gamma^{AB}A_{AB}=0$. We can use it to recover a space-spinor version of the
Levi-Civita connection $\not{\nabla}_{AB}$,

(15)
$\not{\nabla}_{AB}\pi_{C}\equiv\not{D}_{AB}\pi_{C}-\tfrac{1}{4}A_{AB}\gamma_{C}{}^{D}\pi_{D}.$

Indeed,

(16)
$\not{\nabla}_{AB}\epsilon_{CD}=\not{\nabla}_{AB}\gamma_{CD}=0,\quad\gamma^{AB}\not{\nabla}_{AB}\pi_{C}=0,$

and $\not{\nabla}_{AB}$ has vanishing torsion. Moreover,

$\left(\not{\nabla}_{AC}\not{\nabla}_{B}{}^{C}-\not{\nabla}_{B}{}^{C}\not{\nabla}_{AC}\right)\pi^{B}=\tfrac{1}{2}\left(\Psi_{2}+\rho\rho^{\prime}-\sigma\sigma^{\prime}+\mathrm{c.c.}\right)\pi_{A},$

where c.c. denotes the complex conjugation of the expression in the brackets
above. Ultimately, the 2-dimensional Sen and Levi-Civita connections on the
boundary are related in the following way,

(17)
$\not{\mathcal{D}}_{AB}\phi_{C}=\not{\nabla}_{AB}\phi_{C}+\tfrac{1}{4}A_{AB}\gamma_{C}{}^{D}\phi_{D}+Q_{AB}{}^{L}{}_{C}\phi_{L}.$

The connection $\not{\nabla}_{AB}$ is real, i.e.

$\widehat{\not{\nabla}_{AB}\pi_{C}}=-\not{\nabla}_{AB}\widehat{\pi}_{C},$

while,

$\widehat{\not{\mathcal{D}}_{AB}\pi_{C}}=-\not{\mathcal{D}}_{AB}\widehat{\pi}_{C}+\left(Q_{ABC}{}^{D}+\widehat{Q}_{ABC}{}^{D}\right)\widehat{\pi}_{D}+\tfrac{1}{2}\gamma_{C}^{\phantom{C}D}\widehat{\pi}_{D}A_{AB}.$

In the following we will also require an expression for the Hermitian
conjugate of the connection $\not{D}_{AB}$. A direct computation yields

(18)
$\widehat{\not{D}_{AB}\pi_{C}}=-\not{D}_{AB}\widehat{\pi}_{C}+\tfrac{1}{2}\gamma_{C}{}^{D}\widehat{\pi}_{D}A_{AB}.$

In [24], Szabados proposed the following definition of quasilocal angular
momentum associated with $\partial\mathcal{S}$,

(19) $O[N]\equiv-\frac{1}{2\kappa}\oint_{\partial\mathcal{S}}N^{c}A_{c}dS,$

where $N^{a}$ is a divergence-free vector on $\partial\mathcal{S}$,
$A_{c}=\rho_{a}\Pi^{f}{}_{c}\nabla_{f}\tau^{a}$ is the connection 1-form on
the normal bundle of $\partial\mathcal{S}$ and $\kappa=8\pi G$ the
gravitational coupling constant. After inspecting the definition (14) of a
spinor $A_{AB}$ we immediately see that it is in fact a space-spinor
counterpart of $A_{b}$ in the expression above.

In the sequel we will also make use of a Hodge decomposition on
$\partial\mathcal{S}$. Specifically, given any 1-form $V_{a}$ on
$\partial\mathcal{S}$ there exist two functions $f$ and $f^{\prime}$ such that

$V_{a}=\epsilon_{a}{}^{b}\not{\nabla}_{b}f+\not{\nabla}_{a}f^{\prime},$

where $\epsilon_{ab}$ is the $2$-dimensional Riemannian volume form of the
boundary $\partial\mathcal{S}$.

###  2.3. Conformal rescaling of the 2-dimensional Dirac operator

An action of a 2-dimensional (Levi-Civita) _Dirac operator_ on a spinor
$\pi_{A}$ is given by $\not{\nabla}_{A}{}^{B}\pi_{B}$. The purpose of this
subsection is to explore its properties under conformal rescalings of the
inner boundary metric $\sigma_{ab}$. Indeed, according to the uniformization
theorem for compact Riemannian surfaces (_cf._ [12, Theorem 4.4.1]),
$\sigma_{ab}$ is conformal to the spherical metric, i.e.

$\sigma_{ab}=\Omega^{2}\widetilde{\sigma}_{ab},$

where $\Omega:\mathbb{S}^{2}\rightarrow\mathbb{R}$ is a non-negative smooth
function and $\widetilde{\sigma}_{ab}$ a round 2-sphere metric. Given any
holonomic basis $\\{{\partial}_{{\bm{A}}{\bm{B}}}\\}$ on
$\partial\mathcal{S}$, the covariant derivative of $\pi_{A}$ can be expressed
as
$\not{\nabla}_{AB}\pi_{C}=\partial_{{\bm{A}}{\bm{B}}}(\pi_{\bm{C}})-\Gamma_{{\bm{A}}{\bm{B}}}{}^{\bm{Q}}{}_{\bm{C}}\pi_{\bm{Q}}$,
where
$\Gamma_{{\bm{A}}{\bm{B}}}{}^{\bm{Q}}{}_{\bm{C}}\equiv\frac{1}{2}\Gamma_{{\bm{A}}{\bm{B}}}{}^{{\bm{Q}}{\bm{L}}}{}_{{\bm{C}}{\bm{L}}}$
are the spin coefficients associated to the Christoffel symbols of
$\sigma_{ab}$. Since $\sigma_{ABCD}=\Omega^{2}\widetilde{\sigma}_{ABCD}$, one
arrives at

(20)
$\not{\nabla}_{AB}\pi_{C}=\widetilde{\not{\nabla}}_{AB}\pi_{C}-\tfrac{1}{2}\left(\partial_{AB}(\log{\Omega})\widetilde{\sigma}^{EF}{}_{CF}+\partial_{CF}(\log{\Omega})\widetilde{\sigma}^{EF}{}_{AB}-\partial^{EF}(\log{\Omega})\widetilde{\sigma}_{ABCF}\right)\pi_{E},$

where $\widetilde{\not{\nabla}}_{AB}$ is a space-spinor counterpart of the
Levi-Civita connection on $\mathbb{S}^{2}$. Contracting the second and third
indices in the above with $\widetilde{\epsilon}^{AB}=\Omega^{-1}\epsilon^{AB}$
gives

$\not{\nabla}_{A}{}^{B}\pi_{B}=\Omega^{-1}\Big{[}\widetilde{\not{\nabla}}_{A}{}^{B}\pi_{B}-\frac{1}{2}(\partial_{A}{}^{B}(\log{\Omega})\widetilde{\sigma}^{EF}{}_{BF}+\partial_{BF}(\log{\Omega})\widetilde{\sigma}^{EF}{}_{A}{}^{B}-\partial^{EF}(\log{\Omega})\widetilde{\sigma}_{A}{}^{B}{}_{BF})\pi_{E}\Big{]}.$

However, because
$\widetilde{\sigma}_{ABCD}=\frac{1}{2}(\widetilde{\epsilon}_{AC}\widetilde{\epsilon}_{BD}+\widetilde{\gamma}_{AC}\widetilde{\gamma}_{BD})$
and $\widetilde{\gamma}^{AB}e_{AB}(\log\Omega)=0$ ($\widetilde{\gamma}_{AB}$
is orthogonal to $\partial\mathcal{S}$) the equation above reduces to

(21)
$\not{\nabla}_{A}{}^{B}\pi_{B}=\Omega^{-1}\widetilde{\not{\nabla}}_{A}{}^{B}\pi_{B}.$

This can be used to establish a following fact: if $\widetilde{\pi}_{A}$ is a
Dirac eigenspinor on $\mathbb{S}^{2}$ with an eigenvalue
$\lambda\in\mathbb{R}$, i.e.

$\widetilde{\not{\nabla}}_{A}{}^{B}\widetilde{\pi}_{B}=i\lambda\widetilde{\pi}_{A},$

then

(22) $\not{\nabla}_{A}{}^{B}\pi_{B}=i\tfrac{\lambda}{\Omega}\pi_{A},$

where the spin basis transforms in a following way,
$\widetilde{o}_{A}=\tfrac{1}{\sqrt{\Omega}}o_{A},\
\widetilde{\iota}_{A}=\tfrac{1}{\sqrt{\Omega}}\iota_{A}$ and
$\widetilde{\pi}_{0}=\sqrt{\Omega}\pi_{0},\
\widetilde{\pi}_{1}=\sqrt{\Omega}\pi_{1}$.

##  3\. Approximate twistor equation

###  3.1. Setup

Let $\mathfrak{S}_{1}$, $\mathfrak{S}_{3}$ be the spaces of symmetric valence
1 and 3 spinors over the hypersurface $\mathcal{S}$. The (overdetermined)
spatial twistor operator can be defined as follows,

$\mathbf{T}:\mathfrak{S}_{1}\rightarrow\mathfrak{S}_{3},\qquad\mathbf{T}(\kappa)_{ABC}\equiv\mathcal{D}_{(AB}\kappa_{C)},$

and is a space-spinor counterpart of the twistor operator
$\nabla_{A^{\prime}(A}\kappa_{B)}$ (see [5] for more details). The formal
adjoint of $\mathbf{T}$ is given by

$\mathbf{T}^{*}:\mathfrak{S}_{3}\rightarrow\mathfrak{S}_{1},\qquad\mathbf{T}^{*}(\zeta)_{A}\equiv\mathcal{D}^{BC}\zeta_{ABC}-\Omega_{A}{}^{BCD}\zeta_{BCD},$

and allows to define the _approximate twistor operator_
$\mathbf{L}\equiv\mathbf{T}^{*}\circ\mathbf{T}:\mathfrak{S}_{1}\rightarrow\mathfrak{S}_{1}$,

(23)
$\mathbf{L}(\kappa)_{A}\equiv\mathcal{D}^{BC}\mathcal{D}_{(AB}\kappa_{C)}-\Omega_{A}{}^{BCD}\mathcal{D}_{BC}\kappa_{D},$

which is formally self-adjoint —i.e. $\mathbf{L}^{*}=\mathbf{L}$.

Let $\kappa_{A}$ be a solution of the approximate twistor equation
$\mathbf{L}(\kappa)_{A}=0$. The spinors

$\xi_{A}\equiv\tfrac{2}{3}\mathcal{D}_{A}{}^{Q}\kappa_{Q},\quad\xi_{ABC}\equiv\mathcal{D}_{(AB}\kappa_{C)}$

encode independent components of $\mathcal{D}_{AB}\kappa_{C}$. Moreover, one
has that

$\mathbf{L}(\widehat{\xi})_{A}=0.$

Given the set of asymptotically Cartesian coordinates $(x^{\alpha})$ on
$\mathcal{S}$, the position spinor can be defined as follows,

$x_{\mathbf{A}\mathbf{B}}\equiv\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}x^{1}+\mbox{i}x^{2}&-x^{3}\\\
-x^{3}&-x^{1}+\mbox{i}x^{2}\end{array}\right).$

We will consider a solution of the approximate twistor equation with an
asymptotic behaviour of the form

(24)
$\kappa_{\mathbf{A}}=\bigg{(}1+\frac{m}{r}\bigg{)}x_{\mathbf{A}\mathbf{B}}o^{\mathbf{B}}+o_{\infty}(r^{-1/2}).$

A direct computation shows that

(25a)
$\displaystyle\xi_{\mathbf{A}}=\bigg{(}1-\frac{m}{r}\bigg{)}o_{\mathbf{A}}+o_{\infty}(r^{-3/2}),$
(25b)
$\displaystyle\xi_{\mathbf{ABC}}=-\frac{3m}{2r^{3}}x_{(\mathbf{AB}}o_{\mathbf{C})}+o_{\infty}(r^{-5/2}).$

As a consequence of the above asymptotic expansion of $\kappa_{A}$, one
arrives at the following inequality relating the ADM mass of $\mathcal{S}$ and
an integral of concomitants of the spinor $\kappa_{A}$, provided that the
inner boundary $\partial\mathcal{S}$ is a MOTS [13],

(26) $4\pi
m\geq\oint_{\partial\mathcal{S}}n_{AB}\zeta_{C}\widehat{\mathcal{D}^{(AB}\zeta^{C)}}\mbox{d}S,$

where $n_{AB}$ is the outer directed (i.e. towards $r=\infty$) unit normal on
$\partial\mathcal{S}$ as a surface of $\mathcal{S}$ and
$\zeta_{A}\equiv\widehat{\xi}_{A}$. In the sequel we will use a boundary
condition for $\kappa_{A}$ to refine the inequality (26).

###  3.2. A boundary value problem for the approximate twistor equation

Let

(27)
$\mathcal{D}_{A}{}^{Q}\kappa_{Q}=-\frac{3}{2}\widehat{\phi}_{A}\quad\mathrm{on}\quad\partial\mathcal{S},$

where $\phi_{A}$ is a smooth spinorial field. The approximate twistor equation
together with (27) satisfy the _Lopatinskij-Shapiro_ compatibility conditions
(see eg. [7, 28]). This implies that the associated boundary value problem is
elliptic. Moreover, the decay conditions (4a) and (4b) for the first and
second fundamental forms of the initial data make the approximate twistor
operator $\mathbf{L}$ asymptotically homogeneous. In the sequel we will make
use of an operator $\mathbf{B}$, defined in the in the following way,

$\mathbf{B}:\mathfrak{S}_{1}\rightarrow\mathfrak{S}_{1},\qquad\mathbf{B}(\kappa)_{A}\equiv-\sqrt{2}\gamma_{A}{}^{P}\xi_{P}=-\frac{2\sqrt{2}}{3}\gamma_{A}{}^{P}\mathcal{D}^{Q}{}_{P}\kappa_{Q}.$

The equation (27) now becomes

$\mathbf{B}(\kappa)_{A}{}|_{\partial\mathcal{S}}=\sqrt{2}\gamma_{A}{}^{P}\widehat{\phi}_{P},$

and the associated boundary value problem is

(28)
$\mathbf{L}(\kappa)_{A}=0,\quad\mathbf{B}(\kappa)_{A}{}|_{\partial\mathcal{S}}=\sqrt{2}\gamma_{A}{}^{P}\widehat{\phi}_{P}.$

To discuss the solvability of (28) one has to look at the adjoint operators
$\mathbf{L}^{*}$ and $\mathbf{B}^{*}$. A similar computation as in [13] (in
this case the extrinsic geometry of the boundary is non-trivial) leads to the
following:

######  Proposition 1.

If $\partial\mathcal{S}$ is a MOTS, then the boundary value problem

$\mathbf{L}(\kappa)_{A}=0,\qquad\mathbf{B}(\kappa)_{A}|_{\partial\mathcal{S}}=\sqrt{2}\gamma_{A}{}^{P}\hat{\phi}_{P},$

with a smooth spinorial field $\phi_{A}$ over $\partial\mathcal{S}$ admits a
unique solution of the form

(29) $\kappa_{A}=\mathring{\kappa}_{A}+\theta_{A},\qquad\theta_{A}\in
H^{2}_{-1/2},$

with $\mathring{\kappa}_{A}$ given by the leading term in (24) and where
$H^{s}_{\beta}$ with $s\in\mathbb{Z}^{+}$ and $\beta\in\mathbb{R}$ denotes the
weighted $L^{2}$ Sobolev spaces.

###  3.3. Inequality with the connection 1-form on the normal bundle of the
boundary.

The boundary condition (27) allows to simplify the inequality (26) to the
following form

(30) $4\pi
m\geq\sqrt{2}\oint_{\partial\mathcal{S}}\widehat{\phi}^{A}\gamma_{A}{}^{B}\not{\mathcal{D}}_{BC}\phi^{C}dS,$

where $\phi_{A}$ is a free data in the boundary value problem (28). All
quantities in the integral are now intrinsic to the boundary. A relation (17)
between the 2-dimensional Sen and Levi-Civita connections implies that

$\widehat{\phi}^{A}\gamma_{A}{}^{B}\not{\mathcal{D}}_{BC}\phi^{C}=\widehat{\phi}^{A}\gamma_{A}{}^{B}(\not{\nabla}_{BC}\phi^{C})+\tfrac{1}{4}\widehat{\phi}^{A}\gamma_{A}{}^{B}A_{BC}\,\gamma^{CL}\phi_{L}+\widehat{\phi}^{A}\gamma_{A}{}^{B}Q_{BC}{}^{LC}\phi_{L}.$

Using (12), the GHP expression for $Q_{ABCD}$, it is easy to see that
$Q_{AB}{}^{CB}=\rho^{\prime}\iota_{A}o^{C}$ on a MOTS. This can be combined
with the fact that $\gamma_{A}{}^{B}A_{B}{}^{C}\gamma_{C}{}^{D}=-A_{A}{}^{D}$
to yield

(31) $\begin{split}4\pi
m\geq&-\sqrt{2}\oint_{\partial\mathcal{S}}\widehat{\phi}^{A}\gamma_{A}{}^{B}(\not{\nabla}_{B}{}^{C}\phi_{C})\,dS+\sqrt{2}\oint_{\partial\mathcal{S}}\rho^{\prime}|\phi_{0}|^{2}dS-\frac{1}{2\sqrt{2}}\oint_{\partial\mathcal{S}}A_{AB}\widehat{\phi}^{A}\phi^{B}dS.\end{split}$

It should be noted that since $A_{AB}$ is intrinsic to $\partial\mathcal{S}$
the last term can be expressed as

(32)
$-\frac{1}{2\sqrt{2}}\oint_{\partial\mathcal{S}}\widehat{\phi}^{A}\phi^{B}\sigma_{AB}{}^{CD}A_{CD}dS,$

where $\sigma_{ABCD}$ $=\frac{1}{2}(\epsilon_{AC}\epsilon_{BD}$
$+\gamma_{AC}\gamma_{BD})$ is the spinorial counterpart of the 2-dimensional
metric.

##  4\. Mass-quasilocal angular momentum inequality

In this section we present the main result of this article – the mass-
quasilocal angular momentum inequality for the initial data
$(\mathcal{S},h_{ab},K_{ab})$. It is based on a simplification of (31) under
suitable choice of the boundary spinor $\phi_{A}$. A natural condition for
$\phi_{A}$ arises after inspecting the first term on the right-hand side of
(31) – its 2-dimensional Dirac derivative should be controlled. Indeed, we
will proceed with a following choice,

(33)
$\displaystyle\not{\nabla}_{A}{}^{B}\phi_{B}=i\frac{\lambda}{\Omega}\phi_{A},$

where $\Omega$ is a conformal factor relating a metric $\sigma_{ab}$ on
$\partial\mathcal{S}$ with that on a round sphere $\mathbb{S}^{2}$. It can be
showed that with a suitable choice of the conformal rescaling of the spin
basis the equation (33) corresponds to a Dirac eigenproblem on
$\mathbb{S}^{2}$ (see Subsection 2.3 for more details).

The inequality (31) can now be simplified with the use of (33) to the
following form,

(34) $4\pi
m\geq\sqrt{2}\oint_{\partial\mathcal{S}}\rho^{\prime}|\phi_{0}|^{2}dS-\frac{1}{2\sqrt{2}}\oint_{\partial\mathcal{S}}\widehat{\phi}^{A}\phi^{B}\sigma_{AB}{}^{CD}A_{CD}\,dS,$

where the reality of the ADM mass $m$ has been used to eliminate a (purely
imaginary) term with the eigenvalue $\lambda$, i.e.

(35)
$\lambda\oint_{\partial\mathcal{S}}\left(|\phi_{0}|^{2}-|\phi_{1}|^{2}\right)\Omega^{-1}dS=0.$

To make a connection between the second term on the right-hand side of (34)
and the quasilocal angular momentum (19) we will introduce a spinor $N^{AB}$,
defined as follows,

(36)
$N^{AB}\equiv\sigma^{ABCD}\phi_{C}\widehat{\phi}_{D}=\phi_{0}\overline{\phi}_{1}\iota^{A}\iota^{B}-\phi_{1}\overline{\phi}_{0}o^{A}o^{B}.$

One can verify that $N^{AB}$ is real, i.e. $\widehat{N}^{AB}=-N^{AB}$, so it
corresponds to a real 3-vector. Moreover, $\gamma^{AB}N_{AB}=0$ and

(37) $\displaystyle\not{\nabla}_{a}N^{a}$ $\displaystyle=$
$\displaystyle\not{\nabla}_{AB}\left(\sigma^{ABCD}\phi_{C}\widehat{\phi}_{D}\right)=\not{\nabla}_{AB}\left(\phi^{A}\widehat{\phi}^{B}\right)$
$\displaystyle=$
$\displaystyle-(\not{\nabla}_{B}{}^{A}\phi_{A})\widehat{\phi}^{B}+\phi^{A}(\widehat{\not{\nabla}_{A}{}^{B}\phi_{B}})=0,$

where (33) has been used in the last equality. Hence, $N^{a}$ is intrinsic to
the boundary $\partial\mathcal{S}$ and $\not{\nabla}$-divergence-free, so we
can identify it with $N^{a}$ generating the quasilocal angular momentum (19).
With this choice the inequality (34) yields

(38) $4\pi
m\geq\sqrt{2}\oint_{\partial\mathcal{S}}\rho^{\prime}|\phi_{0}|^{2}dS+\frac{\kappa}{\sqrt{2}}O\left[\sigma^{ABCD}\phi_{C}\widehat{\phi}_{D}\right].$

In the remainder of this section we will simplify (38) and express it in terms
of integrals over a round sphere $\mathbb{S}^{2}$ and the eigenspinor of the
$\mathbb{S}^{2}$-Dirac operator.

The Hodge decomposition can be applied to the connection 1-form on the normal
bundle of the boundary $A_{b}$ to yield
$A_{b}=\epsilon_{b}{}^{c}\not{\nabla}_{c}U+\not{\nabla}_{b}U^{\prime}$, where
$U$ is a rotation potential. This allows to simplify the quasilocal angular
momentum term from (38), i.e.

(39)
$\oint_{\partial\mathcal{S}}N^{a}A_{a}dS=\oint_{\partial\mathcal{S}}U\epsilon^{ab}\not{\nabla}_{a}N_{b}dS.$

The spinorial counterpart of the volume element $\epsilon_{ab}$ of
$\partial\mathcal{S}$ is
$\epsilon_{ABCD}=\frac{i}{2}\left(\epsilon_{AC}\gamma_{BD}+\epsilon_{BD}\gamma_{AC}\right)$,
and

(40)
$\epsilon^{ab}\not{\nabla}_{a}N_{b}=\epsilon^{ABCD}\not{\nabla}_{AB}N_{CD}=2\lambda\Omega\gamma^{AB}\phi_{A}\widehat{\phi}_{B}=2\lambda\Omega(|\phi_{0}|^{2}-|\phi_{1}|^{2}),$

where the definition (36) has been taken into account. Inserting this
expression into (39) and using $dS=\Omega^{2}d\mathbb{S}^{2}$ yields

(41)
$\begin{split}\oint_{\partial\mathcal{S}}N^{a}A_{a}dS=2\lambda\oint_{\mathbb{S}^{2}}U\left(|\widetilde{\phi}_{0}|^{2}-|\widetilde{\phi}_{1}|^{2}\right)d\mathbb{S}^{2},\end{split}$

where $\widetilde{\phi}_{\bm{A}}=\sqrt{\Omega}\,\phi_{\bm{A}}$. The relation
between the volume elements of $\partial\mathcal{S}$ and $\mathbb{S}^{2}$ can
also be utilized to write the first term on the right-hand side of (38) in
terms of an integral over $\mathbb{S}^{2}$. Ultimately,

(42) $4\pi
m\geq\sqrt{2}\oint\displaylimits_{\mathbb{S}^{2}}\rho^{\prime}|\widetilde{\phi}_{0}|^{2}\Omega
d\mathbb{S}^{2}+\frac{\kappa}{\sqrt{2}}O\left[\widetilde{\phi},U\right],$

where

$O\left[\widetilde{\phi},U\right]\equiv-\frac{\lambda}{\kappa}\oint\displaylimits_{\mathbb{S}^{2}}U\left(|\widetilde{\phi}_{0}|^{2}-|\widetilde{\phi}_{1}|^{2}\right)d\mathbb{S}^{2},$

i.e. the quasilocal angular momentum term can now be written only in terms of
the geometry of $\mathbb{S}^{2}$ and the rotation potential $U$.

##  5\. Axisymmetric inner boundary and the Dirac-Killing system.

A natural assumption associated with the existence of a well-defined angular
momentum is that the initial data is axisymmetric, i.e. there exists 1-form
$v_{a}$ such that

$\nabla_{(a}v_{b)}=0\quad\mathrm{on}\quad\mathcal{S}.$

If the inner boundary $\partial\mathcal{S}$ is invariant under the action of
the 1-parameter group of isometries generated by $v_{a}$, then
$v_{a}=\Pi_{a}{}^{b}v_{b}$ ($v_{a}$ is intrinsic to $\partial\mathcal{S}$) and
the projection of the Killing equation gives

$\not{\nabla}_{(a}v_{b)}=0\implies\not{\nabla}_{a}v^{a}=0.$

This suggests that a natural choice for the vector $N^{a}$ defining the
quasilocal angular momentum (19) is that it arises as a solution to the
boundary Killing equation, i.e. $\not{\nabla}_{(a}N_{b)}=0$. However, $N^{a}$
has already been constructed from a spinor $\phi_{A}$ satisfying a first-order
Dirac-type equation (33) on $\partial\mathcal{S}$. Hence, a natural question
arises —can such $N_{a}$ be also a Killing vector of the boundary? We will
show that this cannot be the case on a generic $\partial\mathcal{S}$.

In the sequel we will use an adapted system of coordinates $(\psi,\varphi)$ on
the boundary, such that its metric $\sigma_{ab}$ can be written in the
following form,

$\sigma=-R^{2}\left(\tfrac{1}{F^{2}}d\psi\otimes d\psi+F^{2}d\varphi\otimes
d\varphi\right),\quad\psi\in[\psi_{0},\psi_{1}],\quad\varphi\in[0,2\pi],$

where $F=F(\psi)$, $R$ is a constant and the axisymmetric Killing vector is
now proportional to $\partial_{\varphi}$. To avoid the conical singularities
on the poles we will assume that $F(\psi_{0})=F(\psi_{1})=0$. The NP operators
$\delta$ and $\overline{\delta}$ reduce to

$\delta=\tfrac{1}{\sqrt{2}R}\left(F\partial_{\psi}+\tfrac{i}{F}\partial_{\varphi}\right),\qquad\overline{\delta}=\tfrac{1}{\sqrt{2}R}\left(F\partial_{\psi}-\tfrac{i}{F}\partial_{\varphi}\right),$

in this setting. Moreover,
$\alpha-\overline{\beta}=\\-\frac{1}{\sqrt{2}R}\partial_{\psi}F$ (see [6] for
details).

A straightforward computation yields

$\displaystyle\not{\nabla}_{AB}N_{CD}+\not{\nabla}_{CD}N_{AB}$
$\displaystyle=$ $\displaystyle
2\eth^{\prime}\left(\overline{\phi}_{0}\phi_{1}\right)o_{A}o_{B}o_{C}o_{D}-\left(\eth\left(\overline{\phi}_{0}\phi_{1}\right)+\eth^{\prime}\left(\phi_{0}\overline{\phi}_{1}\right)\right)o_{A}o_{B}\iota_{C}\iota_{D}$
$\displaystyle-\left(\eth\left(\overline{\phi}_{0}\phi_{1}\right)+\eth^{\prime}\left(\phi_{0}\overline{\phi}_{1}\right)\right)\iota_{A}\iota_{B}o_{C}o_{D}+2\eth\left(\phi_{0}\overline{\phi}_{1}\right)\iota_{A}\iota_{B}\iota_{C}\iota_{D},$

where we have used a decomposition of vector $N^{a}$ in terms of a spinor
$\phi_{A}$ in accordance with (36). Ultimately, the condition
$\not{\nabla}_{(a}N_{b)}=0$ implies that

(43) $\phi_{0}\overline{\phi}_{1}=icF,\quad c\in\mathbb{R}.$

Additionally, the spinor $\phi_{A}$ satisfies a first-order Dirac-type
equation (33), which can now be written as

(44)
$\begin{split}F\partial_{\psi}\phi_{1}+\tfrac{\phi_{1}}{2}\partial_{\psi}F-\tfrac{i\sqrt{2}\lambda
R}{\Omega}\phi_{0}=0,\\\
F\partial_{\psi}\phi_{0}+\tfrac{\phi_{0}}{2}\partial_{\psi}F-\tfrac{i\sqrt{2}\lambda
R}{\Omega}\phi_{1}=0.\end{split}$

After multiplying the first equation by $\overline{\phi}_{1}$ and the second
by $\overline{\phi}_{0}$ and performing some manipulations one arrives at

$\begin{split}\partial_{\psi}\left(F|\phi_{1}|^{2}\right)+\tfrac{2\sqrt{2}\lambda
RcF}{\Omega}=0,\\\
\partial_{\psi}\left(F|\phi_{0}|^{2}\right)-\tfrac{2\sqrt{2}\lambda
RcF}{\Omega}=0,\end{split}$

where (43) has been used. Hence

(45) $|\phi_{0}|^{2}=\frac{c_{0}}{F}+\frac{2\sqrt{2}\lambda
Rc}{F}\int\displaylimits_{\psi_{0}}^{\psi}\frac{F}{\Omega}dz,\quad|\phi_{1}|^{2}=\frac{c_{1}}{F}-\frac{2\sqrt{2}\lambda
Rc}{F}\int\displaylimits_{\psi_{0}}^{\psi}\frac{F}{\Omega}dz,$

for some $c_{0},c_{1}\in\mathbb{R}$. On the other hand, one can apply
$F\partial_{\psi}$ to both sides of (43) and use (44) to get

$\sqrt{2}\lambda R\left(|\phi_{1}|^{2}-|\phi_{0}|^{2}\right)=2c\Omega
F\partial_{\psi}F.$

After using (45) we obtain the following compatibility condition for $F$,

(46) $\sqrt{2}\lambda R\bigg{(}\frac{c_{1}-c_{0}}{F}-\frac{4\sqrt{2}\lambda
Rc}{F}\int\displaylimits_{\psi_{0}}^{\psi}\frac{F}{\Omega}dz\bigg{)}=2c\Omega
F\partial_{\psi}F.$

Multiplying the above relation by $F$ and differentiating with respect to
$\psi$ we arrive at

(47) $c\left[\partial_{\psi}\left(\Omega
F^{2}\partial_{\psi}F\right)+4R^{2}\lambda^{2}\Omega^{-1}F\right]=0,$

where $F(\psi_{0})=0$ has been used. This equation implies that the metric of
the inner boundary (via the function $F$) depends on the choice of the
eigenvalue $\lambda$. This cannot be the case, as the former arises as part of
the fixed geometric data associated with the initial hypersurface and the
latter from the first-order Dirac-type equation, which is an auxiliary
condition used to simplify the mass inequality. Hence, the only way to solve
(47) is to assume that $c=0$. In this case $\phi_{A}=0$ and the right-hand
side of the mass-quasilocal angular momentum inequality (38) vanishes.

##  6\. Conclusions

We have obtained a new bound for the ADM mass of of _asymptotically
Schwarzschildean_ initial data in terms of the future inner null expansion of
the inner MOTS boundary and its quasilocal angular momentum. Our approach
bears similarities to the one presented in [13], but we extend it here to
allow for boundaries with nontrivial extrinsic geometry. An expression for
quasilocal angular momentum (in the sense of Szabados [24]) has been recovered
in the bound for the ADM mass by assuming a specific type of boundary
condition for the approximate twistor equation– a spinor $\phi_{A}$ solving a
first-order Dirac-type equation. The strategy developed in this work could
also be applied to obtain Penrose-type inequalities with different type of
asymptotics (e.g. asymptotically hyperbolic) —as long as the concept of
quasilocal angular momentum is well-defined.

Acknowledgments. This project began during JK’s visit to the School of
Mathematical Sciences at Queen Mary University of London. He would like to
express his gratitude to the school for its hospitality and acknowledge the
support provided by the Polish National Science Centre through the MINIATURA
project No. 2021/05/X/ST2/01151.

## References

  * [1] A. Alaee and H. Kunduri.  A Penrose-Type Inequality with Angular Momenta for Black Holes with 3-Sphere Horizon Topology  J. Geom. Anal., 33, 231, 2023.
  * [2] P. Anglada.  Penrose-like inequality with angular momentum for minimal surfaces.  Class. Quantum Grav., 35(4):045018, 2018.
  * [3] P. Anglada.  Comments on Penrose inequality with angular momentum for outermost apparent horizons.  Class. Quantum Grav., 37.6:065023, 2020.
  * [4] A. Ashtekar.  Lectures on non-perturbative canonical gravity.  World Scientific, 1991.
  * [5] T. Bäckdahl and J. A. Valiente Kroon.  Approximate twistors and positive mass.  Class. Quantum Grav., 28:075010, 2011.
  * [6] M.J. Cole, I. Rácz, J. A. Valiente Kroon  Killing spinor data on distorted black hole horizons and the uniqueness of stationary vacuum black holes.  Class. Quantum Grav., 35(20), 2018.
  * [7] S. Dain.  Elliptic systems.  In J. Frauendiener, D. Giulini, and V. Perlick, editors, Analytical and Numerical approaches to General Relativity, volume 692 of Lect. Notes. Phys., page 117. Springer Verlag, 2006.
  * [8] S. Dain.  Geometric inequalities for black holes.  Gen. Rel. Grav., 46:1-23, 2014.
  * [9] S. Dain, M. E Gabach-Clement  Geometrical inequalities bounding angular momentum and charges in General Relativity.  Living Reviews in Relativity, 21(5), 2018.
  * [10] G. W. Gibbons, S. W. Hawking, G. T. Horowitz, and M. J. Perry.  Positivite mass theorems for black holes.  Comm. Math. Phys., 88:295, 1983.
  * [11] J. Jaracz and M. Khuri.  Bekenstein bounds, Penrose inequalities and black hole formation.  Physical Review D, 97, 124026, 2018.
  * [12] J. Jost.  Compact Riemann surfaces.  Springer-Verlag, 2006.
  * [13] J. Kopiński, J. A. Valiente Kroon  New spinorial approach to mass inequalities for black holes in general relativity.  Phys. Rev. D, 103, 2021.
  * [14] J. Kopiński, J. Tafel.  The Penrose inequality for perturbations of the Schwarzschild initial data.  Class. Quantum Grav., 37 015012, 2020.
  * [15] J. Kopiński, J. Tafel.  The Penrose inequality for nonmaximal perturbations of the Schwarzschild initial data.  Class. Quantum Grav., 37 105006, 2020.
  * [16] W. Kulczycki, E. Malec.  Extensions of the Penrose inequality with angular momentum.  Physical Review D, 103.6:064025, 2021
  * [17] M. Ludvigsen, J. A. G Vickers.  A simple proof of the positivity of the Bondi mass.  J. Phys. A: Math. Gen., 15.2, L67-L70, 1982.
  * [18] M. Mars.  Present status of the Penrose inequality  Class. Quantum Grav., 26, 193001, 2009.
  * [19] R. Penrose.  Naked singularities.  Ann. N. Y. Acad. Sci., 224, 1973.
  * [20] R. Penrose and W. Rindler.  Spinors and space-time. Volume 1. Two-spinor calculus and relativistic fields.  Cambridge University Press, 1984.
  * [21] R. Schoen and S. T. Yau.  Proof of the positive mass theorem. II.  Comm. Math. Phys., 79:231-260, 1981.
  * [22] P. Sommers.  Space spinors.  J. Math. Phys., 21:2567, 1980.
  * [23] L. B. Szabados.  Two-dimensional Sen connections in general relativity.  Class. Quantum Grav., 11:1833, 1994.
  * [24] L. B. Szabados.  On a class of $2$-surface observables in general relativity.  Class. Quantum Grav., 23, 2291-2302, 2006.
  * [25] L. B. Szabados.  Total angular momentum from dirac eigenspinors.  Class. Quantum Grav., 25:025007, 2008.
  * [26] J. A. Valiente Kroon.  Conformal Methods in General Relativity.  Cambridge University Press, 2016.
  * [27] E. Witten.  A new proof of the positive energy theorem.  Comm. Math. Phys., 80:381, 1981.
  * [28] J. Wloka, B. Rowley, and B. Lawruk.  Boundary value problems for elliptic systems.  Cambridge University Press, 1995.