\section{Introduction}
\label{sec: intro}

Higher spin (HS) symmetries play a central role in the consistency of HS theories
\cite{Fradkin:1986ka,Fradkin:1986qy,Vasiliev:1990en,Vasiliev:2003ev,Bekaert:2005vh},
and serve as important guideline in the understanding of different AdS/CFT models.
In the recent work \cite{Joung:2014qya} the authors of the current article have studied several aspects --- such as minimal coadjoint orbit, minimal representation and Joseph ideal --- related to symmetries
of massless HS fields of symmetric index type, and found a convenient formula for computing trace in HS algebra.
In this article, we continue our previous work and clarify analogous aspects
of a more general type of HS symmetries ---
global symmetries of a theory involving partially-massless (PM) HS fields
\cite{Deser:2001us}.

PM fields are special spectra which exist (or stay irreducible) only in
the backgrounds with non-vanishing constant curvature. More precisely, they are unitary
irreducible in dS background, whereas in AdS they are irreducible but not unitary.
For a given spin $s$, there are $s$ different PM fields labelled by their \emph{depth} $t=0,1,\ldots,s-1$\,,
where $t=0$ corresponds to the massless field.
They are described by gauge potentials that transform with gauge parameters
of lower rank compared to massless fields of the same spin:
\be
	\delta\,\varphi^{\sst (t)}{}_{\mu_{1}\cdots\mu_{s}}
	=\nabla_{(\mu_{1}}\cdots\nabla_{\mu_{t+1}\phantom{|}}\varepsilon^{\sst (t)}{}_{\mu_{t+2}\cdots \mu_{s})}
	+\mathcal O(\Lambda)\,,
\ee
where $\varphi^{\sst (t)}$ is a PM field of depth $t$, $\nabla$ denotes (A)dS covariant derivative
and $\Lambda$ is the cosmological constant.
The propagating degrees of freedom of the
spin $s$ PM field with depth $t>0$ are more than those of the massless spin $s$ but less than those of the massive one.
In the flat space limit, PM field of depth $t$ is decomposed into a collection of massless fields of spin \mt{s, s-1,\ldots, s-t}\,\cite{Zinoviev:2001dt}.
The Killing tensor of a PM field of spin $s$ and depth $t$ are given by $O(d+2)$ Young diagram with two rows of respective length $s-1$ and $s-1-t$ \cite{Skvortsov:2006at}.
Therefore, algebras involving such generators can be interpreted
as global symmetry of a theory involving PM fields.

In order to construct a physically consistent algebra
based on these generators, we need to require
following two conditions:
\begin{itemize}
\item The isometry algebra $\mathfrak{so}_{d+2}$\footnote{In this article, we do not consider any issue related to the reality structure.} is a subalgebra.
This implies that the corresponding HS theory contains gravity sector.

\item All the generators of different spins transform covariantly under isometries.
This is to say that in the corresponding theory, HS fields couple to gravity in a diffeomorphism invariant manner, which is tantamount to the requirement of the equivalence principle to hold.
\end{itemize}
A priori, finding out an algebra satisfying the above conditions is highly non-trivial task
in a bottom-up approach,  as we have to solve for the Lie brackets among
different PM generators (labelled by infinitely many spins and depths) satisfying Jacobi identity\footnote{These kind of construction was discussed recently in \cite{Boulanger:2013zza}.}. Fortunately,
a series of algebras involving PM HS generators and
satisfying the above conditions have been already found in the literature.
The first discussions are from mathematics literature,
the simplest case \cite{Eastwood:2006}
and the generalizations \cite{Gover:2009, Michel:2011},
where certain PM HS algebras are defined as the maximal symmetries of higher order Laplace operators
\be
	\Box^{\ell}\,\phi=0\,,
	\label{ho single}
\ee
generalizing the definition of the massless HS algebra ($\ell=1$)\,.
In physics, the unfolded formulation of the corresponding PM HS spectra
has been proposed in \cite{Bekaert:2005vh}
and  explored in \cite{Alkalaev:2007bq,Alkalaev:2011zv}.
In \cite{Bekaert:2012vt,Bekaert:2013zya,Basile:2014wua}, the Flato-Fronsdal theorem was generalized to higher-order singletons and arguments were provided for a conjecture of AdS/CFT duality, relating a Vasiliev-type theory with PM symmetric tensor fields to the $O(N)$ model at a multicritical isotropic Lifshitz point.

In this paper, we revisit the PM algebras
including the aforementioned ones
from various points of view.
We first define the PM algebras as particular cosets of
universal enveloping algebra (UEA) of the isometry algebra.
This requires to identify proper ideals of the UEA
and such ideals depend on a parameter $\lambda$\,.
For generic value of $\lambda$\,, the algebra
is spanned by PM generators of any even depths.
When $\lambda$ takes an integer value $\ell$\,,
it truncates into the algebra associated with \eqref{ho single}
allowing only the generators of depth smaller than $2\,\ell$\,.
We also show how the same algebras can be constructed
making use of oscillators and Howe duality.
The oscillator description allows us to perform more concrete calculations
and we derive the explicit form of trace  and bilinear form for these PM algebras.
The explicit expression of bilinear form reveals
a new structure of PM algebras.
For (half-)integer values of $\lambda$ greater than $d/2$\,, with (odd)even $d$,
the algebras truncate into
\emph{finite dimensional} ones.

The paper is organized as follows. In Section \ref{2}, we provide generalities on the global symmetries of free PM HS fields and the vector space of PM HS algebra. In Section \ref{3}, we define PM HS algebras as cosets of
the UEA of (A)dS isometry algebra.
In Section \ref{4}, we present Howe duality and oscillator realization of PM HS algebras.
In Section \ref{5}, the trace for the PM HS algebras are defined
and the bilinear forms are calculated.
In Section \ref{6}, we show
the finite-dimensional truncation of algebras for special values of $\lambda$.
In Section \ref{7}, the relations of the PM HS algebras to
conformal higher spin theory are discussed.
Appendix \ref{A} contains computational details.

\section{Partially-Massless Fields and Their Global Symmetries}
\label{2}

The global-symmetry generators, or Killing tensors, corresponding to PM fields can be
identified by analyzing the relevant Killing equations.
They were first identified in the frame-like formulation in \cite{Skvortsov:2006at}.
Below, we shall rederive the PM Killing tensors in the metric-like ambient formulation.
As it has been shown in \cite{Alkalaev:2011zv,Joung:2012rv},
the Fierz system of the PM field of spin $s$ and depth $t$
has the following gauge equivalence relation,
\be
	\delta_E\,\Phi^{\sst (t)}(X,U)=(U\cdot \partial_{X})^{t+1}\,E^{\sst (t)}(X,U)\,,
\ee
where PM gauge parameters $E^{\sst (t)}$ are traceless tensors:
\be
	\partial_{U}^{2}\,E(X,U)=0\,.
	\label{traceless}
\ee
As usual, being ambient-space fields, $\Phi^{\sst (t)}$ and $E^{\sst (t)}$
should satisfy the tangentiality and homogeneity conditions:
\ba
	&X\cdot\partial_{U}\,\Phi^{\sst (t)}(X,U)=0\,,
	\qquad
	&\left(X\cdot\partial_{X}-U\cdot\partial_{U}+2+t\right)\,\Phi^{\sst (t)}(X,U)=0\,,\nn
	&X\cdot\partial_{U}\,E^{\sst (t)}(X,U)=0\,,
	\qquad
	&\left(X\cdot\partial_{X}-U\cdot\partial_{U}-t\right)\,E^{\sst (t)}(X,U)=0\,,
\ea
in order to be reduced to (A)dS Fierz system.
The global symmetries are given by parameters, that satisfy the Killing equation,
\mt{\delta_E\, \Phi^{\sst (t)}=0}\,, which combined with the tangentiality and homogeneity conditions
defines the following system,\footnote{One
can equally describe PM fields in the Stueckelberg formulation where the gauge
transformation has the standard form,
\be
	\delta_E\,\Phi^{\sst (t)}(X,U)=U\cdot \partial_{X}\,E^{\sst (t)}(X,U)\,,
\ee
whereas the tangentiality condition is modified to
\be
	(X\cdot\partial_{U})^{t+1}\,E^{\sst (t)}(X,U)=0\,.
\ee
Therefore, in the Stueckelberg formulation, we get conditions
equivalent to \eqref{PM Killing} but $X$ and $U$ exchanged.}
\be
	X\cdot\partial_{U}\,E^{\sst (t)}=0\,,
	\quad
	(U\cdot\partial_{X})^{t+1}\,E^{\sst (t)}=0\,,
	\quad
	(X\cdot\partial_{X}-U\cdot\partial_{U}-t)\,E^{\sst (t)}=0\,.
	\label{PM Killing}
\ee
It is obvious to see from here that the PM Killing tensors of depth $t$ correspond to the
$(t+1)$-dimensional representation of $\mathfrak{sp}_{2}$\,:
\ba
	&[h,e]=2\,e\,,\qquad [h,f]=-2\,f\,,\qquad  [e,f]=h\,,&\nn
	&e=U\cdot\partial_X\,,
	\qquad f=X\cdot\partial_U\,,\qquad
	h=X\cdot\partial_X-U\cdot\partial_U\,.&
\ea
We shall also see the relevance of this $\mathfrak{sp}_{2}$ representation later in
the other construction.
The conditions \eqref{PM Killing} and \eqref{traceless} are solved by
\be
	E^{\sst (t)}(X,U)=\sum^{\infty}_{r=0}
	\frac1{r!\,(r-t)!}\,X^{\sst M_{1}}\cdots X^{\sst M_{r}}\,
	U^{\sst N_{1}}\cdots U^{\sst N_{r-t}}\,
	M^{\sst (r,t)}_{\sst M_{1}\cdots M_{r},N_{1}\cdots N_{r-t}}\,,
\ee
with the tensor $M^{\sst (r,t)}$ taking values in the $\{r,r-t\}$ Young diagram of $O(d+2)$\,:
\be
	M_{\sst M_{1}\cdots M_{r},N_{1}\cdots N_{r-t}} \sim\,
	\parbox{145pt}{\begin{tikzpicture}
	\draw (0,0) rectangle (5,0.4);
	\draw (0,0) -- (0,-0.4) -- (3,-0.4) -- (3,0);
	\draw [dotted] (5,0) -- (5,-0.4) -- (3,-0.4);
	\draw [very thin, gray] (0.4,0.4) -- (0.4,-0.4);
	\draw [very thin, gray] (4.6,0.4) -- (4.6,0);
	\draw [very thin, gray] (2.6,-0.4) -- (2.6,0);
	\node at (2.5,0.2){${\st r}$};
	\node at (1.5,-0.2){$\st r-t$};
	\node at (4,-0.2){$\st t$};
	\end{tikzpicture}}\,.
	\label{PM gen}
\ee
Therefore, any theory involving a PM field of spin $s$ and depth $t$
should admit global symmetry
containing generator $M^{\sst (s-1,t)}$\,.

\section{Coset Construction from Universal Enveloping Algebra}
\label{3}

Similarly to the massless HS algebra,
the algebras involving PM HS generators can be approached from the universal enveloping algebra (UEA)
of $\mathfrak{so}_{d+2}$\,.
One of the advantages of UEA construction is that the two physical consistency conditions mentioned in Introduction are automatically satisfied.
Let us remind that the massless HS algebra can be obtained as the coset,
\be
	hs(\mathfrak{so}_{d+2})=\mathcal U\ \Big/
	\left({\tiny \yng(2)}\oplus\,{\tiny{\yng(1,1,1,1)}}\ \right),
\ee
where $\mathcal U$  denotes the UEA of $\mathfrak{so}_{d+2}$,
while (---) denotes the Joseph ideal, generated by the following elements in
$\mathfrak{so}_{d+2}\odot\mathfrak{so}_{d+2}$\,:
\be
	J_{ab}:= M_{(a}{}^{c}{}\odot M_{b)c}-\frac{\eta_{ab}}{d+2}\,
	M^{cd}\odot M_{cd}\sim {\tiny\yng(2)}\,,
	\qquad
	J_{abcd} := M_{[ab}\odot M_{cd]}
	\sim {\tiny\yng(1,1,1,1)}\ .
\ee
This procedure determines the eigenvalues of all Casimir operators of $\mathfrak{so}_{d+2}$.
In particular, the quadratic Casimir is given by
\be
	C_{2}:=\frac12\,M_{ab}\odot M^{ba}=-\frac{(d+2)(d-2)}4\,.
\ee
Let us now consider deformations of the above construction
where we take a quotient with the ideal
generated by only one of $J_{ab}$
 and $J_{abcd}$\,.
 In such cases, the quadratic Casimir $C_{2}$ remains arbitrary while the other Casimirs
are fixed as functions of $C_{2}$ as shown in \cite{Iazeolla:2008ix,Boulanger:2011se}.
Hence, one can take a further quotient with $C_{2}-\nu$\,.
In this way, we have two ideals labelled by a continuous parameter $\nu$ as
\be
	\mathcal I(\nu)=
	\left( \, {\tiny{\yng(1,1,1,1)}}\ \oplus (C_{2}-\nu) \right),
	\qquad
	\mathcal J(\nu)=
	\Big( \,{\tiny{\yng(2)}} \oplus (C_{2}-\nu) \Big)\,.
	\label{ideals}
\ee
Since these ideals contain fewer elements
than the Joseph ideal,
the corresponding coset algebras contain more generators than the original one. These one-parameter
families of algebras have already been investigated in
\cite{Michel:2011} and \cite{Boulanger:2011se}, respectively.
In the case of $\mathcal I(\nu)$\,, we get the global symmetry algebra of
a theory involving PM HS fields,
whereas for $\mathcal J(\nu)$\,, we get the symmetry of
a theory involving mixed symmetry HS fields.
Leaving aside the discussion about the mixed symmetry case,
in the current article we shall consider only the ideal $\mathcal I(\nu)$ and corresponding PM HS algebras.

Let us resume the discussion with the coset algebra,
\be\label{Anu}
 	\mathcal A_{\lambda}=\mathcal U/\,
	\mathcal I(\nu_{\lambda})\,,
	\qquad
	\nu_{\lambda}=-\frac{(d-2\,\lambda)(d+2\,\lambda)}{4}\,,
\ee
where we have parametrized $\nu$ in terms of $\lambda$
for later convenience.
The ideal generated by $J_{abcd}$ contains all $GL(d+2)$-tensors
in $\mathcal U$ having more than two rows, hence
the algebra $\mathcal A_{\lambda}$ is spanned by the generators
which have the symmetry of two-row $GL(d+2)$ Young diagrams.
When decomposed into the traceless $O(d+2)$ tensors, it is
given by
\be
	\mathcal A_{\lambda}\simeq \bigoplus_{p=0}^{\infty}\,
	K_{p}\,,
	\label{PMdecomposition}
\ee
where $K_{p}$ are the vector space of
the Killing tensors associated with the PM field
of any spin and depth \mt{t=2p}\,:
\be
	K_{p}=\bigoplus_{r=0}^{\infty}\
	\parbox{145pt}{\begin{tikzpicture}
	\draw (0,0) rectangle (5,0.4);
	\draw (0,0) -- (0,-0.4) -- (3,-0.4) -- (3,0);
	\draw [dotted] (5,0) -- (5,-0.4) -- (3,-0.4);
	\draw [very thin, gray] (0.4,0.4) -- (0.4,-0.4);
	\draw [very thin, gray] (4.6,0.4) -- (4.6,0);
	\draw [very thin, gray] (2.6,-0.4) -- (2.6,0);
	\node at (2.5,0.2){${\st r+2p}$};
	\node at (1.5,-0.2){$\st r$};
	\node at (4,-0.2){$\st 2p$};
	\end{tikzpicture}}\,.
\ee
Hence,  the generators of $\mathcal A_{\lambda}$
are not only the massless Killing tensors $K_{0}$
but also the PM ones $K_{p\ge1}$ corresponding to the depth $2p$
\cite{Alkalaev:2007bq} (see also \cite{Alkalaev:2014nsa} for related discussions).

As has been shown in \cite{Michel:2011}\,,
when $\lambda$ takes integer values $\ell
=1,2,\ldots,$
the algebra $\mathcal A_{\ell}$ develops the ideals $\mathfrak{q}_{\ell}$
which consist of PM Killing tensors of even depth not smaller than $2\,\ell$\,:
\be
	\mathfrak{q}_{\ell} \simeq
	\bigoplus_{p=\ell}^{\infty} K_{p}\,.
\ee
As a result, the algebra is decomposed into
\be
	\mathcal A_{\ell}=\mathfrak{p}_{\ell}
	\oplus
	\mathfrak{q}_{\ell}\,,
\label{Aell}
\ee
where $\mathfrak{p}_{\ell}$ is the coset algebra:
\be
	\mathfrak{p}_{\ell}=\mathcal A_{\ell}/\mathfrak{q}_{\ell}\simeq
	\bigoplus_{p=0}^{\ell-1} K_{p}\,,
	\label{vect p}
\ee
containing generators of even depth smaller than $2\,\ell$\,.
Therefore, we obtain two algebras $\mathfrak{q}_{\ell}$ and $\mathfrak{p}_{\ell}$ for $\lambda=\ell$\,:
the former algebra contains generators of even depths not smaller than $2\,\ell$\,\footnote{If realized field theoretically, $\mathfrak{q}_{\ell}$ should correspond to  global symmetry of an exotic interacting theory involving
PM fields of even depths not smaller than $2\,\ell$\,.},
whereas the latter one is associated with the theory
involving PM fields of even depths smaller than $2\,\ell$\,.
The usual massless HS algebra corresponds to $\mathfrak{p}_1$\,.

Let us focus now on the coset algebra $\mathfrak{p}_{\ell}$
which contains the usual HS symmetry in the $\ell=1$ case.
In fact, it can be obtained directly from $\mathcal U$ as the coset,
\be
	\mathfrak{p}_{\ell}=\mathcal U/\,\mathcal I_{\ell}\,,
	\label{coset p}
\ee
where the ideal $\mathcal I_{\ell}$ is
given by \cite{Gover:2009}
\be
 	\mathcal I_{\ell}=\left(
	\parbox{74pt}{\begin{tikzpicture}
	\draw (0,0) rectangle (2.5,0.26);
	\draw [very thin, gray] (0.26,0) -- (0.26,0.26);
	\draw [very thin, gray] (2.24,0) -- (2.24,0.26);
	\node at (1.25,0.14) {$\st 2\,\ell$};
	\end{tikzpicture}}
	\oplus\,{\tiny{\yng(1,1,1,1)}}\ \right).
	\label{I ell}
\ee
The vector space \eqref{vect p} of the $\mathfrak{p}_{\ell}$ can be
straightforwardly identified using this ideal,
or equivalently from the corresponding coadjoint orbit defined by
\be
	W^{(a_{1}}{}_{b_{1}}\,W^{c_{1}}{}_{d_{1}}\cdots\,
	W^{a_{\ell}}{}_{b_{\ell}}\,W^{c_{\ell})}{}_{d_{\ell}}\,
	\eta^{b_{1}d_{1}}\,\cdots \eta^{b_{\ell}d_{\ell}}=0\,,
	\qquad
	W^{a[b}\,W^{cd]}=0\,.
	\label{pm orbit}
\ee
Any polynomials of $W^{ab}$ subject to the above conditions
are dual to the generators of $\mathfrak p_{\ell}$
so give rise to the vector space \eqref{vect p}\,.
The values of the quadratic Casimir \eqref{Anu}
can be also calculated from the ideal $\mathcal I_{\ell}$ by
inquiring the consistency of the ideal.

The algebra $\mathfrak{p}_{\ell}$ corresponds in fact to the symmetry algebra of the $\ell$-th power of
Laplacian operator, or in other words the symmetry of the action,
\be\label{boxell}
	S[\phi]=\int d^{d}x\,\phi\,\Box^{\ell}\,\phi\,.
\ee
This symmetry has been studied for $\ell=1,2$ in
\cite{Eastwood:2002su, Eastwood:2006}
and for arbitrary positive integers $\ell$ in  \cite{Gover:2009}.
The aspect of holographic correspondence has been
also explored recently in \cite{Bekaert:2013zya}.

For generic values of $\lambda$\,, the algebra $\mathcal{A}_\lambda$
is generated by the Killing tensors of all even depth PM HS fields.
Hence, if there exists a theory based on this field content, then $\mathcal{A}_\lambda$ may well
correspond to its global symmetry algebra.
However, the viability of such a theory is not obvious:
since the value of the quadratic Casimir used for the quotienting does not correspond any more to
that of a free scalar field on the boundary, it is not clear what might be the CFT dual.
The ideal $\mathfrak{q}_\ell$ appearing for the integer value of $\lambda=\ell$ corresponds to an (higher-order) on-shell condition
for the boundary scalar,
but $\mathcal{A}_{\lambda}$ does not admit such an ideal, hence the boundary scalar seems to be an off-shell one with a fixed conformal weight. It will be interesting to examine the relevant Flato-Fronsdal theorem
to see whether this  CFT might be dual to the putative PM HS theory, but it is beyond the scope of the current paper.\footnote{We are indebted to the anonymous referee for pointing out this subtle point.}

\section{Howe Duality and Oscillator Realization}
\label{4}
The coset construction from UEA provides a concise definition of the PM symmetries,
but for explicit calculations we may consider another description of PM algebra.
Similarly to massless case, PM algebra $\mathcal{A}_{\lambda}$ can be
defined using a reductive dual pair correspondence (introduced in HS context in \cite{Vasiliev:2003ev}), aka Howe duality, as
\be
	\big(\,Sp(2)\,,\,O(d+2)\,\big)\subset Sp(2(d+2))\,.
\ee
The starting point is again the metaplectic representation of $\mathfrak{sp}_{2(d+2)}$
which can be realized by the star product algebra with the product,\footnote{As
 in \cite{Joung:2014qya}, we denote $\mathfrak{sp}_2$ indices by greek letters and $\mathfrak{so}_{d+2}$ indices by roman letters. The $\mathfrak{sp}_2$ indices are rised and lowered using antisymmetric invariant tensors $\epsilon_{\alpha\beta}$ and $\epsilon^{\alpha\beta}$, while $\mathfrak{so}_{d+2}$ indices are rised and lowered by symmetric metric tensor of $\mathfrak{so}_{d+2}$ (which can be taken as Kroneker delta, since we work with the complex algebra, not a particular real form) and its inverse.}
\be
	(f\star g)(y)
	=\exp\Big(\frac12\,\epsilon_{\alpha\beta}\,\partial_{y_{\alpha}}\!\cdot \partial_{z_{\beta}}\Big)\,
	f(y)\,g(z)\,\Big|_{z=y}\,.
\ee
The mutual stabilizers $\mathfrak{sp}_2$ and $\mathfrak{so}_{d+2}$ correspond to
\be
	K_{\alpha\beta}=y_{\alpha}\!\cdot y_{\beta}\,,
	\qquad
	M_{ab}=y_{\alpha a}\,y^{\alpha}{}_{b}\,,
	\label{M rep}
\ee
and the off-shell HS algebra  $\widetilde{hs}(\mathfrak{so}_{d+2})$
is defined as the centralizer of $\mathfrak{sp}_2$\,:
any element $f(y)$ in $\widetilde{hs}(\mathfrak{so}_{d+2})$ satisfies
\be
	[\,f(y)\ \overset{\star},\  K_{\alpha\beta}\,]=0\,.
\ee
By solving this condition, one can show that $\widetilde{hs}(\mathfrak{so}_{d+2})$
is spanned by polynomials of $M_{ab}$\,.

The usual (on-shell) massless HS algebra is
the coset of $\widetilde{hs}(\mathfrak{so}_{d+2})$
by the $\mathfrak{sp}_2$-triviality relation, $K_{\alpha\beta}\sim0$\,.
This relation actually defines the trivial representation of   $\mathfrak{sp}_2$\,,
which is dual to the singleton representation of $\mathfrak{so}_{d+2}$\,.
In the following, we shall consider generalization
of the last step and define the PM algebras $\mathcal A_{\lambda}$
and $\mathfrak p_{\ell}$\,.

\paragraph{PM algebra $\bm{\mathcal A_{\lambda}}$}
In order to obtain $\mathcal A_{\lambda}$ from $\widetilde{hs}(\mathfrak{so}_{d+2})$\,,
we can quotient the latter with
the $\mathfrak{sp}_{2}$ Casimir relation,
\be
	C_{2}(\mathfrak{sp}_2)=
	\frac12\,K_{\alpha\beta}\star K^{\alpha\beta}
	\sim (1-\lambda)(1+\lambda)\,,
	\label{soCasimir}
\ee
which also fixes the quadratic Casimir of $\mathfrak{so}_{d+2}$
in terms of $\lambda$ to \eqref{Anu}
since in $\widetilde{hs}(\mathfrak{so}_{d+2})$ the two Casimirs are related by
\be
	K_{\alpha\beta}\star K^{\alpha\beta}
	+ M_{ab}\star M^{ba}=-\frac{(d+2)(d-2)}{2}\,.
\ee
In the case of $\lambda=\ell$\,,
the relation \eqref{soCasimir}
gives an indecomposable representation
consisting of the $(\ell+1)$-dimensional representation
and the infinite-dimensional representation with lowest (or heighest) weight
$\ell+1$ (or $-(\ell+1)$)\,.
Hence, the algebra $\mathcal A_{\ell}$
becomes the semi-direct-sum of  the coset algebra $\mathfrak{p}_{\ell}$
and the ideal algebra $\mathfrak q_{\ell}$ as in \eqref{Aell}.

\paragraph{PM algebra $\bm{\mathfrak{p}_{\ell}}$}

In order to obtain directly $\mathfrak{p}_{\ell}$ from $\widetilde{hs}(\mathfrak{so}_{d+2})$\,,
we can quotient the latter with the following equivalence relations
\be\label{ell rep}
	K_{++}\sim 0\,,
	\quad
	K_{--}^{\, \ell} \sim 0\,,
	\quad
	\left(K_{+-}+\ell-1\right) \sim 0\,,
\ee
which define the $(\ell+1)$-dimensional representation
of $\mathfrak{sp}_2$\,.
This relations have asymmetric form with respect to $\mathfrak{sp}_2$ indices, but
one can show in fact that any other relations defining the $(\ell+1)$-dimensional representation
such as,
\be
	K_{++}^{\,1+k}\sim 0\,,
	\quad
	K_{--}^{\,\ell-k} \sim 0\,,
	\quad
	\left(K_{+-}+\ell-k-1\right) \sim 0\,,
	\qquad[\,k=0,1,\ldots,\ell\,]\,,
	\label{PM equiv}
\ee
give the same algebra $\mathfrak{p}_\ell$.
Instead of \eqref{PM equiv},
one may also use the equivalence relation\footnote{The same $\mathfrak{sp}_{2}$ tensors have been used in decomposing generators of off-shell HS algebra into PM ones in \cite{Alkalaev:2007bq}.}
\be
K_{(\alpha_1\alpha_2}\star K_{\phantom{(}\alpha_3\alpha_4}\star
\cdots\star K_{\alpha_{2\ell-1}\alpha_{2\ell})}\sim 0\,,\label{PM covariant condition}
\ee
which generalizes the relation $K_{\alpha\beta}\sim 0$ of the massless HS algebra
to PM cases in a manifestly $\mathfrak{sp}_2$ covariant manner.
However, this relation is weaker than \eqref{PM equiv}
since the relations $K_{(\alpha_{1}\alpha_{2}}\star\cdots\star K_{\alpha_{2n-1}\alpha_{2n})}\sim0$
with $n=0,1,\ldots,\ell-1$ also imply \eqref{PM covariant condition}.
Therefore, only complemented by \eqref{soCasimir},
the relation
\eqref{PM covariant condition}
becomes equivalent to \eqref{PM equiv}
so defines the PM algebra $\mathfrak{p}_{\ell}$\,. The condition \eqref{soCasimir} implies multiplicity one for each generator of given spin and depth.
In the following we shall use the relations \eqref{soCasimir}
and \eqref{PM covariant condition}
to derive the trace of PM algebras.
Instead, in Appendix \ref{A}, we make use of relations \eqref{ell rep} for alternative derivation of the trace.

\subsection*{Decomposition of traceful generators into PM ones}

For concreteness, let us look at
the generators of $\mathcal A_{\lambda}$
(and $\mathfrak{p}_{\ell}$) more closely
and show that they indeed correspond to PM Killing tensors given in \eqref{PM gen}.
We begin with $\widetilde{hs}(\mathfrak{so}_{d+2})$ spanned by
\be
	\tilde M^{\sst (r)}_{a_{1}\cdots a_{r},b_{1}\cdots b_{r}}
	=y_{a_{1}\alpha_{1}}\,y_{b_{1}}{}^{\alpha_{1}}\,\cdots\,y_{a_{r}\alpha_{r}}\,y_{b_{r}}{}^{\alpha_{r}}\,,
	\label{tf tensor}
\ee
which are traceful tensors corresponding to $GL(d+2)$ $\{r,r\}$ Young diagrams.
The traceful tensor $\tilde M^{\sst (r)}$ \eqref{tf tensor} can be decomposed into traceless tensors $M^{\sst (r,2t)}$
with $t=0,\ldots, [r/2]$\,,
\be
	M^{\sst (r,2t)}_{a_{1}\cdots a_{r},b_{2t+1}\cdots b_{r}}
	=\tilde M^{\sst (r)}_{a_{1}\cdots a_{r},b_{1}\cdots b_{r}}\,
	\eta^{b_{1}b_{2}}\cdots \eta^{b_{2t-1}b_{2t}}
	-({\rm traces})\,,
\ee
which can be interpreted as PM generators of spin $r+1$ and depth $2t$\,.
These tensors can be conveniently handled by contracting them with
\be
	\tilde W^{ab}=\tilde w^{a}_{\alpha}\,\tilde w^{\alpha b}\,.
	\label{W para}
\ee
Then, the traceless generators of $\widetilde{hs}(\mathfrak{so}_{d+2})$ given by
\be
	M^{\sst (r,2t)}(\tilde W)
	= M^{\sst (r,2t)}_{a_{1}\cdots a_{r},b_{2t+1}\cdots b_{r}}\,
	\eta_{b_{1}b_{2}}\cdots \eta_{b_{2t-1}b_{2t}}\,
	\tilde W^{a_{1}b_{1}}\cdots \tilde W^{a_{r}b_{r}}\,,
	\label{M W}
\ee
can be expressed solely through  traces of $M_{ab}$'s and $\tilde W^{ab}$'s,
\be
	\la M^{p}\,\tilde W^{q}\,M^{r}\,\cdots\,\tilde W^{s}\ra\,,
	\label{traces}
\ee
as they do not have any free index.
Using the following identity based on \eqref{M rep} and \eqref{W para},
\be
	A\,B\,A=\frac12\,A\,\la A\,B\ra\,,
	\qquad [\,A,B=M\ {\rm or}\ \tilde W\,],
	\label{M W id}
\ee
one can show that any trace of type \eqref{traces} will
be broken down to a function of
\be
	\la M^{2}\ra \la \tilde W^{2}\ra\,,
	\qquad
	\la M\,\tilde W\ra\,,\qquad \la (M^{2}\,\tilde W^{2})^{p}\ra\,.
\ee
Here again, one can show that the last object with  $p\ge2$
is not independent and can be expressed in terms of the rest.
Hence, $M^{\sst (r,2t)}(\tilde W)$ in \eqref{M W} are polynomials of $\la M^{2}\ra \la \tilde W^{2}\ra$\,,
$\la M\,\tilde W\ra$ and $\la M^{2}\,\tilde W^{2}\ra$\,.
Rewriting them as $\star$ polynomials,
we can replace $\la M^{2}\ra$ by $2(1-\lambda^{2})$ in $\mathcal A_{\lambda}$\,.
In this way, any element of $\mathcal A_{\lambda}$\,, namely PM generators,
can be written as $\star$ polynomials of  $\la M\,\tilde W\ra$ and $\la M^{2}\,\tilde W^{2}\ra$
with $\la \tilde W^{2}\ra$-dependent coefficients:
\ba
M^{\sst (r,2t)}(\tilde W)
&\sim&\la M^{2}\,\tilde W^{2}\ra^{\star t}\star
\la M\,\tilde W\ra^{\star (r-2t)}+\,\la \tilde W^{2}\ra\,\Big(
a\,\la M^{2}\,\tilde W^{2}\ra^{\star (t-1)}\star\la M\,\tilde W\ra^{\star (r-2t-2)}\nn
&&+\,b\,\la M^{2}\,\tilde W^{2}\ra^{\star(t-2)}\star
\la M\,\tilde W\ra^{\star (r-2t-4)}+\cdots \Big)+\cdots,
\ea
where $a$ and $b$ are some constants.
Here, one can notice that
the $\star$ polynomial corresponding to $M^{\sst (r,2t)}$
has maximal orders $r-2t$ and $t$ in $\la M\,\tilde W\ra$ and $\la M^{2}\,\tilde W^{2}\ra$\,, respectively.

 In $\mathcal A_{\lambda}$\,,
 any $\star$ powers of  $\la M^{2}\,\tilde W^{2}\ra$ can appear
 so PM generators of any even depth are present there.
However in $\mathfrak{p}_{\ell}$\,,
one can show using \eqref{PM covariant condition} that
the $\ell$-th $\star$ power of $\la M^{2}\,\tilde W^{2}\ra$
reduces to lower order ones:
\be
	 \la M^{2}\,\tilde W^{2}\ra^{\star \ell}\sim
	  a\,\la \tilde W^{2}\ra\,\la M^{2}\,\tilde W^{2}\ra^{\star (\ell-1)}+\cdots+
	 \la \tilde W^{2}\ra^{\ell}\,,
\ee
with some constant $a$\,.
Therefore, the generators of depth not smaller than $2\,\ell$ become
linear combinations of generators of lower depth.
In order to remove this redundancy,
one can impose on $\tilde W^{ab}$ the orbit condition \eqref{pm orbit}, which is equivalent to
\be
	W^{ab}=w^{a\alpha}\,w^{b}{}_{\alpha}\,,\qquad
w_{(\alpha_1}\!\cdot w_{\alpha_2}\,\cdots\,w_{\alpha_{2\ell-1}}\!\cdot w_{\alpha_{2\ell})}=0\,.
    \label{orbit W para}
\ee
Here, $W^{ab}$ denotes $\tilde W^{ab}$ satisfying the above condition.
With this, any higher than $\ell-1$ $\star$-powers of $\la M^{2}\,W^{2}\ra$ vanish in $\mathfrak{p}_\ell$\,:
\be
	 \la M^{2}\,W^{2}\ra^{\star n}=0
	 \qquad
	 [\,n=\ell,\ell+1,\ldots\,]\,.
\ee
It will be convenient to require another condition  on $W^{ab}$\,,
\be
	\la W^{2}\ra=0
	\qquad \Leftrightarrow
	\qquad
	w_\alpha{}\cdot w_{\gamma}\,w^{\gamma}\cdot w_\beta=0\,,
	\label{W2=0}
\ee
to hold. Note that this requirement does not project any component of
generators in \eqref{M W} so we do not loose any generality by imposing it.
To summarize, with \eqref{orbit W para} and \eqref{W2=0},
any elements of $\widetilde{hs}(\mathfrak{so}_{d+2})$
are projected to the generators of $\mathcal{A}_\lambda$ (or $\mathfrak{p}_\ell$) in a
simple manner as
\be
    M^{\sst (r,2t)}(W)=\left\{\begin{array}{ccc}
    &\la M^{2}\,W^{2}\ra^{\star t}\star
    \la M\,W\ra^{\star (r-2t)}\,,\qquad &t=0,1,2\ldots\\
    & 0\,, \qquad &t=\ell,\ell+1,\ldots
    \end{array}\right.,
\ee
where the second case only holds for $\mathfrak{p}_\ell$\,.

\section{Trace and Bilinear Form}
\label{5}
Analogously to the case of massless HS algebra,
it will be useful to define a trace in PM HS algebras,
from which the bi-/tri-linear forms and the structure constants can be derived.
We define the trace as the coefficient of identity, similarly to the
massless HS algebra,
\be
{\rm Tr}\left[c_0+c_a T^a\right]=c_0\,,
\label{Tr Def}
\ee
where $T^a$ denotes a generator of PM HS algebra.
The expression of trace for any element in the algebra
can be obtained in principle by only using
the definition \eqref{Tr Def} as in the case of massless HS algebra \cite{Joung:2014qya},
and  the trace formula for PM HS algebra $\mathfrak{p}_{\ell}$
is obtained in that way for a few lower $\ell$'s in Appendix \ref{A}.

In the following, we will derive the bilinear form of
$\mathcal A_{\lambda}$, that should manifest the appearance of the ideal $\mathfrak{q}_{\ell}$
and the corresponding coset $\mathfrak{p}_{\ell}$
when $\lambda=\ell$\,.

\subsection*{Non-Gaussian projector}

In the massless case, it turns out \cite{Joung:2014qya} that the trace obtained from the definition \eqref{Tr Def}
admits a simple representation,
 \be
	{\rm Tr}_{\mathfrak{p}_{1}}[f]=(\Delta\star f)(0)\,,
	\label{trace}
\ee
where $\Delta(y)$ is the \emph{projector} introduced in \cite{Vasiliev:2001wa,Vasiliev:2004cm}\footnote{Strictly speaking it is not a projector, since $\Delta\star\Delta$ is not well defined as shown in \cite{Vasiliev:2004cm}. See also recent discussion in the context of $(A)dS_2$ algebras in \cite{Alkalaev:2014qpa}.}.
The function $\Delta(y)$ giving the trace is not unique as shown in \cite{Joung:2014qya}
even though the definition \eqref{Tr Def} is unambiguous.

Generalizing this massless result to PM one, we assume that the
trace of the PM algebra $\mathcal A_{\lambda}$ admits also an analogous expression,
\be
	{\rm Tr}_{\mathcal A_{\lambda}}[f]
	=(\Delta_{\lambda}\star f)(0)\,,
	\label{trace lambda}
\ee
with the conditions
\be
	[\,\Delta_{\lambda}\,\overset{\star},\,K_{\alpha\beta}\,]=0=
	[\,\Delta_{\lambda}\,\overset{\star},\, M_{ab}\,]\,.\label{commutation condition}
\ee
The above conditions are satisfied by the ansatz,
\be
\Delta_{\lambda}=\Delta_{\lambda}(z)\,,\qquad z=K_{\alpha\beta}\,K^{\alpha\beta}\,.
\ee
In the massless case, we also impose the condition that
$K_{\alpha\beta}\star \Delta_{1}(z)=0$\,.
In the case of $\mathcal A_{\lambda}$\,, we impose instead
\ba\label{l-1condition}
\left(
C_2(\mathfrak{sp}_2)-(1-\lambda)(1+\lambda)\right)
\star\Delta_{\lambda}(z)=0\,.\label{c2condition}
\ea
In order to determine such $\Delta_{\lambda}$\,, we assume
\be
	\Delta_{\lambda}(z)=\sum_{n=1}^{\infty}\,d_{n}(\lambda)\,\Pi_{n}(z)\,,
	\label{projector def}
\ee
with $\Pi_{n}(z)$ satisfying
\be\label{covariant projector definition}
K_{(\alpha_1\alpha_2}\star K_{\alpha_3\alpha_4}\star\cdots\star K_{\alpha_{2n-1}\alpha_{2n})}
\star \Pi_{n}(z)=0\,.
\ee
We shall show below that \eqref{c2condition}
fixes uniquely the coefficients $d_{n}(\lambda)$ and
determines the projector.
Since
the condition \eqref{PM covariant condition}
should hold for $\lambda=\ell$\,,
the coefficients in \eqref{projector def} will satisfy $d_{n\ge\ell}(\ell)=0$\,.
The identification of $\Pi_{n}(z)$ is in parallel with $\Delta_{1}$ in \cite{Vasiliev:2004cm},
and the equation \eqref{covariant projector definition} for $\Pi_{n}$
takes the following form,
\be\label{Delta equation}
\left(2\,z\,\partial_z^2+(d+1)\,\partial_z+1\right)^{n}\Pi_{n}(z)=0\,.
\ee
Clearly $\Pi_{m}(z)$ with $m<n$ also solves the equation for $\Pi_{n}$
but since we are considering a linear combination of them,
without a loss of generality
we can take a particular solution of \eqref{Delta equation},
\be\label{Delta cov}
\Pi_{n}(z)=\int_{-1}^{1}ds\,(1-s^2)^{\frac{d}2-n}\,e^{i\,s\,\sqrt{2\,z}}\,,
\ee
for which $\left(2\,z\,\partial_z^2+(d+1)\,\partial_z+1\right)^{n-1}\Pi_{n}(z)\neq 0$.
With the above $\Pi_{n}(z)$'s,
the condition \eqref{c2condition}
fixes the coefficients $d_{n}(\lambda)$ in the formula \eqref{projector def} as
\be
d_n(\lambda)=c\,\frac{(1-\lambda)_{n-1}\,(1+\lambda)_{n-1}}{(\frac32)_{n-1}\,(n-1)!}\,
\frac1{d-2n}\,,
\label{cl}
\ee
up to an overall factor $c$ which can in turn be fixed from the
normalization condition ${\rm Tr}(1)=1$
(or, equivalently $\Delta_{\lambda}(0)=1$).
Notice that the coefficient $d_{n}(\lambda)$ is ill-defined for $d=2n$\,.
This problem is in fact an artifact of non-Gaussian projector.
In terms of Gaussian projector, this factor will disappear
and in the end all formulas are well-defined.
Details of computation for \eqref{cl} can be found in Appendix \ref{A}\,.
Remark also that
the dependence of $\lambda$ and $d$ in \eqref{cl}
is such that there exist
special values of $\lambda$ where
the series 	\eqref{projector def}
 acquire special properties.
These special cases deserve more careful analysis since
they may imply the appearance of an ideal in the algebra $\mathcal A_{\lambda}$\,.
For concrete treatment, we shall explicitly calculate the
bilinear form using the trace formula.
However, as in the massless case, the
non-Gaussian dependence on oscillators in the projector
 complicates the computations.
Fortunately,
we can find an alternative projector which gives
the same trace but admits a treatable expression.

\subsection*{Gaussian projector}

The integral form \eqref{Delta cov} of $\Pi_{n}$
is simple but not very useful for
actual computations as it involves
\be
	\exp\left(i\,s\sqrt{\frac{y_{\alpha}\cdot y_{\beta}\,y^{\alpha}\cdot y^{\beta}}2}\right).
\ee
However, using the same method as the one described in the Appendix of \citep{Joung:2014qya}, one can show that the expression
\be\label{Projector1}
P_{n}=
\int_0^1dx\, x^{\frac12}\, (1-x)^{\frac{d-2}2-n}\, e^{-2\,\sqrt{x}\,y_+\cdot\, y_-}\,
\ee
is equivalent to $\Pi_{n}$\, in the following sense:
\be
	(\Pi_{n}\star f)(0)=(d-2n)\,(P_{n}\star f)(0)\,.
\ee
Here, one can see that the problematic factor $d-2n$ cancels out
when working with the Gaussian projector.
Since $\Delta_\lambda$ will be used only within the setting of \eqref{trace lambda}, we can replace $\Pi_n$ by $(d-2n) P_n$ in the formula of the projector.
Evaluating the sum in \eqref{projector def}, we obtain
new projector as
\be\label{Projector}
	D_{\lambda}=N_{\lambda} \int_0^1 dx\, x^{\frac12}\, (1-x)^{\frac{d-4}2}\,
	{}_2F_1\left(1+\lambda\,,\, 1-\lambda\,;\, \frac32\,;\, \frac1{1-x}\right)
	e^{-2\,\sqrt{x}\,y_+\cdot\,y_-}\,,
\ee
with
\be
N_{\lambda}=\frac{(-1)^{\lambda-1}\,\Gamma(d+1)}{2^{d-1}\,
\Gamma(\frac{d}2-\lambda)\,\Gamma(\frac{d}2+\lambda)}\,.
\label{NG projector}
\ee
Let us emphasize again that this projector is equivalent to $\Delta_\lambda$
since
\be
    {\rm Tr}(f)=\left(\Delta_\lambda \star f\right)(0)=
    \left(D_\lambda \star f\right)(0)\,,
    \qquad \forall\,f\in \widetilde{hs}(\mathfrak{so}_{d+2})\,.
    \label{trace D}
\ee
The apparent advantage of this projector is that it involves
a Gaussian function of oscillators and therefore is convenient for actual calculations.
Hence, we will use the expression \eqref{Projector}
in the computation of invariant bilinear form of the PM HS algebra.
The trace defined with \eqref{Projector}
can be equally obtained starting from the definition \eqref{Tr Def}
as we have demonstrated for $\mathfrak{p}_{\ell}$
in Appendix \ref{A}.

\subsection*{Bilinear form}

Analogously to massless case \cite{Joung:2014qya}, we consider now the invariant bilinear form of the algebra given as a trace from star-product of generating functions of all PM generators:
\be
B(W_1,W_2)={\rm Tr}(M(W_1)\star M(W_2))
={\rm Tr}\left(e^{y_+\cdot W_1\cdot y_-}\star e^{y_+\cdot W_2\cdot y_-}\right).
\label{bilin}
\ee
We use \eqref{trace D} for the trace,
hence consider first
\be
\label{starrho}
\frac1{G^{\sst (2)}(\rho,W)}=e^{\rho\, y_{+}\cdot\,y_{-}}\star e^{y_+\cdot W_1\cdot y_-}\star e^{y_+\cdot W_2\cdot y_-}\,\big|_{y_{\alpha}=0}\,,
\ee
which  can be evaluated using the
star product composition formula derived in \cite{Didenko:2003aa}
and used in \citep{Joung:2014qya} as
\ba\label{G2}
G^{\sst (2)}(\rho,W)=\frac{\det_{\sst N\times N}\!\left[\frac12\,\frac{1+\rho}{1-\rho}
\,\frac{2+W_1}{2-W_1}\,\frac{2+W_2}{2-W_2}+\frac12\right]}
{\det_{\sst N\times N}\!\left[\frac12\,\frac{1+\rho}{1-\rho}+\frac12\right]
\det_{\sst N\times N}\!\left[\frac12\,\frac{2+W_1}{2-W_1}+\frac12\right]
\det_{\sst N\times N}\!\left[\frac12\,\frac{2+W_2}{2-W_2}+\frac12\right]}\,.
\ea
This expression involves determinants of
$N\times N$ matrices, which can be simplified to
\be
G^{\sst (2)}(\rho,W)
=\det{\!}_{\sst N\times N}\!\left[1+\rho\,\frac{W_1+W_2}2+\frac{W_1\,W_2}4\right].
\ee
Further simplification can be made by
defining the \emph{dual} matrices $V_{ij}$ as
\be
	(V_{ij})_{\alpha\beta}=(w_i)_{\alpha}{}^{A}\,(w_{j})_{\beta A}
	\qquad
	\left[\,(W_{i})^{AB}=(w_{i})^{\alpha A} (w_{i})_{\alpha}{}^{B}\,\right],
\ee
and this enables us to rewrite (using Sylvester's determinant theorem) the $N\times N$ determinant
as $4\times 4$ one as
\be
G^{(2)}(\rho,W)=\det{}_{\sst 4\times 4}\!
\begin{bmatrix}
  1+\frac{\rho}2\,V_{11}+\frac14\,V_{12}\,V_{21}  &\quad \frac{\rho}2\,V_{12}+\frac14\,V_{12}\,V_{22} \\
    \frac{\rho}2\,V_{21}       					  & 1+\frac{\rho}2\,V_{22}
\end{bmatrix}.
\ee
Using the parameterization \eqref{orbit W para} and
the condition \eqref{W2=0},
the above reduces to
\be
G^{\sst (2)}(\rho,W)=
\det{\!}_{\sst 2\times 2}\left[1+\frac{1-\rho^2}4\,V_{12}\,V_{21}
+\frac{\rho}2\,V_{11}-\frac{\rho\,(1-\rho^2)}8\,
V_{12}\,V_{22}\,V_{21}\right],
\ee
which can be straightforwardly evaluated to give
\be\label{G2PM}
G^{\sst (2)}(\rho,W)=
\left(1+\frac{1-\rho^2}8\,\la W_1\,W_2 \ra\right)^2
+\frac{\rho^2\,(1-\rho^2)}{16}\,\la W_1^2\,W_2^2 \ra\,.
\ee
Finally, using \eqref{Projector}, \eqref{starrho} and \eqref{G2PM},
the invariant bilinear form can be expressed as
\ba\label{Bilinear}
B(W_1,W_2)
=N_{\lambda}\,\int_0^1 dx\ \frac{x^{\frac12}\,
(1-x)^{\frac{d-4}2}\, {}_2F_1(1+\lambda,\ 1-\lambda;\ \tfrac32;\ \tfrac1{1-x})}
{\left(1+\frac{1-x}8\,\la W_1\,W_2\ra\right)^2
+\frac{x\,(1-x)}{16}\,\la W_1^2\,W_2^2\ra}\,.
\ea
Let us make a few observations on the bilinear form \eqref{Bilinear}.
First of all, it is manifestly symmetric in $\lambda \to -\lambda$, which is a straightforward consequence of the symmetry of the quadratic Casimir.
Moreover, when $\lambda=\ell$
the bilinear form become degenerate for the generators
$M^{\sst (r,2t)}$ with $t\ge \ell$\,.
This degeneracy is less manifest since
$\la W_{1}\,W_{2}\ra^{r-2t}\,\la W_{1}^{2}\,W_{2}^{2}\ra^{t}$
does not exactly correspond to the contribution of $M^{\sst (r,2t)}$
--- in other words, the bilinear form is not diagonal.
Taking the simplest example of $\lambda=\ell=1$ case,
the $W^{2}$-order term of the bilinear form is given by
\ba
	&&
	a\left(\la W_{1}\,W_{2}\ra^{2}
	-\frac{4}{d}\,\la W_{1}^{2}\,W_{2}^{2}\ra\right)
	=W_{1}^{a_{1}b_{1}}\,W_{1}^{a_{2}b_{2}}\,
	W_{2}^{c_{1}d_{1}}\,W_{2}^{c_{2}d_{2}}\times\nn
	&&\qquad \times
	\left[{\rm Tr}\left(M^{\sst (2,0)}_{a_{1}a_{2},b_{1}b_{2}}\,
	M^{\sst (2,0)}_{c_{1}c_{2},d_{1}d_{2}}\right)
	+{\rm Tr}\left(M^{\sst (2,2)}_{a_{1}a_{2}}\,M^{\sst (2,2)}_{c_{1}c_{2}}\right)
	\eta_{b_{1}b_{2}}\,\eta_{d_{1}d_{2}}\right],
\ea
with some constant $a$\,.
By decomposing the LHS of the equation into traceless
tensors, one can verify that the PM generators of the second depth
have vanishing bilinear from:
\be
	{\rm Tr}\left(M^{\sst (2,2)}_{a_{1}a_{2}}\,M^{\sst (2,2)}_{c_{1}c_{2}}\right)=0\,,
\ee
hence form  an ideal.
Likewise, the appearance of the ideal $\mathfrak{q}_{\ell}$ and the
corresponding coset algebra $\mathfrak{p}_{\ell}$ in $\mathcal A_{\ell}$
can be checked with diagonalization. And as we discussed before,
 the ideal part --- having vanishing bilinear form --- can be
 conveniently discarded by contracting with $W^{ab}$ satisfying
 the orbit condition \eqref{pm orbit}.
 We shall keep this condition in the following section.

\section{Finite dimensional truncations}
\label{6}
The normalization factor $N_{\lambda}$ \eqref{NG projector} of the bilinear form \eqref{Bilinear}
contains
\be
	\frac1{\Gamma(\frac d2-\lambda)}\,,
\ee
which vanishes for $\lambda-\frac d2=0,1,\ldots$\,.
Remind that $N_{\lambda}$ is fixed with the condition
${\rm Tr}(1)=1$ which means the norm of the identity
is one. When $N_{\lambda}$ itself
vanishes, the only way to keep the normalization of the identity
is that the integral \eqref{Bilinear} should diverge
for the identity part.
In this case, if there exists generators which give
finite integral, then their norm will vanish
and  the algebra develops a new ideal corresponding
to such generators.

In order to study these new ideals and corresponding coset algebras,
it is useful to use the transformation of hypergeometric function,
\be
{}_2F_1(1+\lambda,\,1-\lambda;\,\tfrac32;\,\tfrac1{1-x})=
(\tfrac{x}{x-1})^{\lambda-1}{}_2F_1(\tfrac12-\lambda,\,1-\lambda;\,\tfrac32;\,\tfrac1{x})\,,
\ee
to rewrite the bilinear form as
\ba\label{Bilinear odd d}
B(W_1,W_2)
=N_{\lambda}\,(-1)^{\lambda-1}\int_0^1 dx\,
\frac{x^{\lambda-\frac12}\ (1-x)^{\frac{d-2}2-\lambda}\,
{}_2F_1\left(\frac12-\lambda\,,\,
1-\lambda\,;\, \frac32\,;\, \frac1{x}\right)}
{\left(1+\frac{1-x}8\,\la W_1\,W_2\ra\right)^2
+\frac{x\,(1-x)}{16}\,\la W_1^2\,W_2^2\ra}\,.
\ea
Focusing on the $M^{\sst (r,2t)}$ part, which corresponds to
\be
	\la W_1\,W_2\ra^{r-2t}\,\la W_1^2\,W_2^2\ra^{t}\,,
	\label{Wr Wt}
\ee
the integral is proportional, up to finite coefficient, to
\ba
&&\int_0^1 dx\, x^{\lambda-\frac12+t}\,
(1-x)^{\frac{d-2}2-\lambda+r-t}\,
{}_2F_1\left(\tfrac12-\lambda\,,\, 1-\lambda\,;\, \tfrac32\,;\, \tfrac1{x}\right)\nn
&&=\,\Gamma\left(\frac{d}2-\lambda+r-t\right)
\sum_{n=0}^{\infty}\,\frac{(\frac12-\lambda)_{n}\,
(1-\lambda)_{n}\,\Gamma(\lambda+t+\frac12-n)}
{n!\,(\frac32)_{n}\,\Gamma(\frac{d+1}2+r-n)}\,.
\label{integral}
\ea
This integral diverges when the first factor $\Gamma(\frac d2-\lambda+r-t)$ diverges
while the summation part is always finite.
Hence, for a fixed value of $\lambda$ with vanishing $N_{\lambda}$\,:
\be
\lambda=\frac d2+k\,,\qquad k=0,1,\ldots
\label{trunc lambda}
\ee
the algebra develops an infinite-dimensional ideal corresponding to
the generators $M^{\sst (r,2t)}$ with $r-t>k$\,. The coset algebra,
denoted henceforth by $\mathfrak{f}_{k}$\,,
 is finite dimensional one with the generators $M^{\sst (r,2t)}$ satisfying
\be
	r-t\le k\,.
\ee
The Young diagram of $M^{\sst (r,2t)}$ has $r$
and $r-2t$ boxes in the first and second row respectively,
so  the total number  is $2(r-t)$\,.
Therefore, the algebra $\mathfrak{f}_{k}$ consists
of the generators whose Young diagram contains no more than $2k$ boxes.
This situation is very similar to the one with massless HS algebras in three dimensions\footnote{In fact, the three dimensional HS algebra $hs[\lambda]$ can be interpreted as the partially-massless HS algebra in $(A)dS_2$, and allows for finite-dimensional truncations, interpreted as partially-massless algebras \cite{Alkalaev:2014qpa}.}, and five dimensions \cite{Joung:2014qya,Manvelyan:2013oua}  where the corresponding HS algebras  are the symmetries of symmetric tensor representations of $\mathfrak{sl}_4$ or, equivalently, (anti-)self-dual (spin-)tensors of $\mathfrak{so}_6$ described by rectangular three-row Young diagrams.
For example, $\mathfrak{f}_{k}$ with $k=3$ have (we omit the identity generator, corresponding to Young diagram with no boxes)
\ba
&\left\{\ \small\yng(1,1)\ ,\ \small\yng(2)\ \right\},\qquad
\left\{\ \small\yng(2,2)\ ,\ \small\yng(3,1)\ , \ \small\yng(4)\ \right\},&\nn
&\left\{\ \small\yng(3,3)\ ,\ \small\yng(4,2)\ ,
\ \small\yng(5,1)\ ,\ \small\yng(6)\ \right\},&
\ea
where we have organized the generators in terms of the number of boxes.
 Regrouping the same set according to spins, we get
\ba
\left\{\ \small\yng(1,1)\ \right\},\qquad
\left\{\ \small\yng(2,2)\ ,\ \small\yng(2)\ \right\},\qquad
\left\{\ \small\yng(3,3)\ ,\ \small\yng(3,1)\ \right\},\qquad \nn
\left\{\ \small\yng(4,2)\ ,\ \small\yng(4)\ \right\},\qquad
\left\{\ \small\yng(5,1)\ \right\},\qquad
\Big\{\ \small\yng(6)\ \Big\}\,.
\ea
The dimension of $\mathfrak{f}_{3}$ is
$\left(\tfrac{(d+1)(d+2)(d+6)}6\right)^{2}$\,.
For generic integer value of $k$\,,
the dimension of $\mathfrak{f}_{k}$ is given again by  a perfect square $M_{k}^{2}$ with
\be
	M_{k}=\frac{(d+1)_{k-1}(d+2k)}{k!}\,.
\ee
This suggests that the algebra $\mathfrak{f}_{k}$ is isomorphic
to $\mathfrak{gl}_{M_{k}}$\,,
an endomorphism of a $\mathfrak{so}_{d+2}$ representation
with dimension $M_{k}$\,.
It turns out that this representation corresponds to
the finite-dimensional one-row $\mathfrak{so}_{d+2}$ Young diagram
with length $k$\,:
\be
	\begin{tikzpicture}
	\draw (0,0) rectangle (4,0.4);
	\draw [very thin, gray] (0.4,0) -- (0.4,0.4);
	\draw [very thin, gray] (3.6,0) -- (3.6,0.4);
	\node at (2,0.2) {$k$};
	\end{tikzpicture}\,.
	\label{fin k rep}
\ee
This can be verified by comparing the quadratic Casimir of \eqref{fin k rep}
--- given by $k(k+d)$ ---
with that of \eqref{trunc lambda}\,: the latter gives
\be
	C_{2}(so(d+2))=-\frac{(d-2\,\lambda)(d+2\,\lambda)}{4}
	=k(k+d)\,,
\ee
hence they coincide as expected.
This proves that the finite-dimensional coset algebra $\mathfrak{f}_{k}$
corresponds to the symmetries of the finite-dimensional representation
\eqref{fin k rep}.
The latter representation can be also understood from
the boundary point of view:
they should correspond to the solution (sub)space of higher-order Laplace equation for scalar field $\Box^{d/2+k}\phi=0$\,.
Since any higher-order Laplace equation, $\Box^{\ell}\phi=0$\,,
is equivalent to the set of equations in the $(d+2)$-dimensional ambient space:
\be
    \partial_X^2\,\Phi(X)=0\,,\qquad
    \left(X\cdot \partial_X+\frac{d}2-\ell\right)\,\Phi(X)=0\,,
    \label{amb sys}
\ee
one can study its solutions using the above system (see e.g. \citep{Bekaert:2013zya})\footnote{In \eqref{amb sys},
we can consider the transformation
\be
	\Phi(X)\quad\to\quad (X^{2})^{\ell-d}\,\Phi(X)\,,
\ee
to flip the sign of $\ell$\,. Higher order singletons
are usually described with the positive sign --- that is, with homogeneity
$-(d+2\ell)/2$ --- together with $(X^{2})^{\ell}\,\Phi(X)=0$\,. Equivalent description is given by equation $(\partial_X^2)^{\ell}\Phi(X)=0$ supplemented with condition $(X^2)\,\Phi(X)=0\,$. These conditions can be solved to get higher order conformal wave operators (GJMS operators) in arbitrary background \cite{Manvelyan:2007tk}.}.
For $\ell=\frac d2+k$\,, the ambient field $\Phi$ has homogeneity
degree $k$\,, any $k$-th order polynomials in $X^a$ dual to
traceless tensors \eqref{fin k rep} become solutions of the system \eqref{amb sys}.

Let us conclude this section with a short summary.
In this section,  using the bilinear form \eqref{Bilinear odd d}, we have identified the finite dimensional PM HS algebras  corresponding to the symmetry of the symmetric tensor representation \eqref{fin k rep}  of (A)dS algebra.
Retrospectively,  these algebras could be identified with no reference to the bilinear form but
only scanning possible finite dimensional representations of (A)dS algebra whose tensor square
satisfy the condition of PM algebras.

\section{PM Algebras and CHS Symmetries}
\label{7}

Conformal Higher Spin (CHS) theory is an interacting theory of conformal higher-spin fields
 ---
that is, defined by Fradkin-Tseytlin free action \cite{Fradkin:1985am} ---
of all spins \cite{Segal:2002gd,Bekaert:2010ky} and its symmetry coincides
with the symmetry of massless
HS theory in one higher dimensions.
It has been shown that $(d+1)$-dimensional free conformal spin-$s$  field
can be decomposed, around $(A)dS_{d+1}$\,, into the set of spin-$s$ PM fields
with all depths \cite{Joung:2012qy,Metsaev:2014iwa,Nutma:2014pua}:
\be
	{\rm CHS}_{s}=\bigoplus_{t=0}^{s+\frac{d-5}2} {\rm PM}_{(s,t)}\,,
\ee
where ${\rm PM}_{(s,t)}$ with $t\ge s$ are massive fields.
Therefore, CHS symmetry itself can be regarded as the symmetry of PM fields of any depths.
Then a natural question arises: whether the even depth PM symmetries
discussed in this paper ($\mathcal{A}_{\lambda}, \mathfrak{p}_{\ell}$ or $\mathfrak{f}_{k}$)
can be embedded as a subalgebra in CHS symmetry.
First of all, it is rather straightforward to see that the even depth PM part of CHS symmetry form a subalgebra.
CHS symmetry in $(A)dS_{d+1}$ is isomorphic to  the massless HS algebra
in $(A)dS_{d+2}$\,, so can be realized by $d+2$ sets of oscillators, $Y_{\alpha A}=(y_{\alpha a},z_{\alpha})$\,. Now considering the map,
\be
	\rho\,:\ \left(\,y_{\alpha a}\ ,\, z_\alpha\,\right)
	\quad \to \quad \left(\,y_{\alpha a}\ ,\,  -z_\alpha\,\right),
\ee
which is an automorphism
 for $hs(\mathfrak{so}_{d+3})$
the  $\rho$ invariant space
of $hs(\mathfrak{so}_{d+3})$ forms a subalgebra
and such space is generated by even depth PM generators.

The next question is to what this symmetry corresponds.
One of the simplest way to answer this question is
by examining the Howe duality in the oscillator construction.
The quotient relation of massless HS algebra in $AdS_{d+2}$
can be translated in $AdS_{d+1}$ to
\be
	K_{\alpha\beta}+z_{\alpha}\,z_{\beta}\sim 0\,,
\ee
where $K_{\alpha\beta}$ is the $\mathfrak{sp}_{2}$ generators dual
to $\mathfrak{so}_{d+2}$ --- not $\mathfrak{so}_{d+3}$\,.
The above relation simply means that the dual $\mathfrak{sp}_{2}$
carries the representation realized by $z_{\alpha}\,z_{\beta}$\,.
By calculating the quadratic Casimir, we get
\be
	C_{2}(\mathfrak{sp}_{2})\sim \frac12\,z_{\alpha}\,z_{\beta}\star z^{\alpha}\,z^{\beta}
	=\frac{3}{4}\,,
\ee
and can recognize that it coincides with the $C_2$ of  $\mathcal A_{\frac12}$\,.
This proves that $\mathcal A_{\frac12}$ is the even depth subalgebra of CHS symmetry.

The existence of even depth PM subalgebra $\mathcal A_{\frac12}$ inside of CHS symmetry has an interesting implication towards a unitary truncation of CHS theory.
CHS theory is considered to be non-unitary since its linearized spectrum is described by
a higher derivative action.
The non-unitarity of the latter become clear when it is decomposed into the actions of PM fields around $(A)dS_{d+1}$:
\be
	S^{\sst\rm CHS}_{s}=\sum_{t=0}^{s+\frac{d-5}2} (-1)^{t}\,
	S^{\rm\sst PM}_{(s,t)}\,,
	\label{CHS PM}
\ee
where even depth PM fields and odd depth PM fields have relatively negative sign for kinetic terms.
One may wonder whether the above CHS action
can be truncated into a unitary one by simply selecting
positive sign part.
Although it is not certain whether this truncation might be consistent,
one can immediately identify one necessary condition for
the consistency: the symmetry of the truncated
spectrum should form a subalgebra of the original symmetry.
For example,
the spectrum of Conformal Gravity (CG)
decomposes into massless spin two and PM spin two,
and their kinetic terms have relatively negative sign.
The symmetry of CG --- that is the conformal symmetry ---
contains the symmetry of Einstein Gravity --- that is the isometry ---
as a subalgebra,
and the action of CG
written as \eqref{CHS PM}
can be consistently truncated to Einstein Hilbert one
\cite{Metsaev:2007fq,Maldacena:2011mk,Deser:2012qg}.
Now coming back to CHS theory, if CHS action can be
truncated into a theory of even depth PM fields,
we can escape from the non-unitarity problem
at least at the linear level.
The fact that $\mathcal A_{\frac12}$
corresponds to a subalgebra of CHS symmetry suggests
that such a truncation
might be viable.

Another interesting question is that whether other even depth PM symmetries we have encountered in this paper can
be still considered as a truncation of (a variant of) CHS theory.
Although we do not have a proof, we do not see how the other PM algebras can be embedded in CHS symmetry.
Moreover, we do not see even the massless HS symmetry,
that is $\mathfrak{p}_{1}$\,, can be so:
interpreted differently, we do not see
how usual HS symmetry can be
embedded inside of the same symmetry in one higher dimensions.
Considering the embedding of $\mathfrak{p}_{\ell}$ within conformal-like symmetry,
there is a chance that $\mathfrak{p}_{\ell}$ is supplemented with odd depth generators
to form a new symmetry.
Corresponding theories may coinside with Weyl-like theories of partially-massless fields of both even and odd depth lower than a given number that were considered in \cite{Joung:2012qy}. Conversely, Weyl-like theories of \cite{Joung:2012qy} may have even-depth unitary truncations that make use of the algebras $\mathfrak{p}_{\ell}$.

\acknowledgments

We thank S. Rey  for useful discussions.
The work of EJ was supported by Basic Science Research
Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A6A3A04056670)
and by the Russian Science Foundation grant 14-42-00047 associated with Lebedev Institute.
The work of KM was supported by the BK21 Plus Program funded by the Ministry of Education (MOE, Korea) and National Research Foundation of Korea (NRF).


\section{Introduction}

The High Altitude Water Cherenkov (HAWC) Gamma-ray Observatory is a newly commissioned instrument on volcano Sierra Negra in the state of Puebla, Mexico~\cite{bib:hawcsensi}.
Here we describe the design and capabilities of a high-level (i.e., event-level) analysis framework.
The framework is a set of software objects that
provides particular functionality for statistical analysis of event-level data and that facilitates the development of analysis programs by end users.
Event-level data are air shower events observed in the HAWC detector that have been reconstructed (i.e., a direction
and other shower observables have been estimated)
and have passed some low-level selection triggers.
Our framework facilitates statistical analysis using the well-known likelihood formalism \cite{Neyman}.
The framework supports a model-centric approach for a user to design an analysis program.
That is, a user of the framework mainly concentrates on the physics model or detector model that she wants to test.
The software interfaces to the physics model, the detector response, the binned data, and the calculation of likelihood values are standardized.
The actual details of the physics models, the details of the detector response, and the binning of the data are defined by the user through these standardized interfaces.

Some noteworthy design features of the framework are:
\begin{itemize}
  \item The physics model definition has no restriction with regard to energy spectrum, source spatial morphology, or number of sources.
  \item The physics model can be arbitrarily parameterized and any number of parameters can be optimized with a maximum likelihood routine.
  \item The detector response model can be parameterized enabling details of the detector response to be estimated directly from data
using a well-known gamma-ray source (e.g., the Crab Nebula).
  \item The data is binned on the sky using Healpix bins~\cite{bib:healpix}. The data can be binned in an arbitrary number of other dimensions (e.g., event quality, energy, etc.) that are defined by the user. The binned data (the size of which is typically $~10^8$ bins, neglecting the time dimension) is efficiently stored on disk and handled in memory.
  \item A routine to calculate a profile maximum likelihood is supported.
  \item The framework is fully compatible with the Multi-mission Maximum Likelihood (3ML) architecture \cite{Vianello_icrc} which enables the
incorporation of multiple instruments (e.g., HAWC and Fermi/LAT) in a joint-likelihood fit.
\end{itemize}

In the following sections we elucidate the design features of the
HAWC event-level analysis framework, also called \textbf{Li}kelihood \textbf{F}itting \textbf{F}ramework (LiFF).

\section{Likelihood Formalism for Statistical Significance}

The main purpose of the framework  is to facilitate the calculation of the likelihood of a particular physics model and detector response model given an observed set of data.
Currently, the framework is designed with the minimum time window as one transit (one day)~\footnote{We currently have other software routines for dealing with shorter
time windows. We plan to expand the capabilities of this framework to shorter time windows in the future.}.
To substantially reduce the computational workload, we bin the data -- we calculate a binned likelihood value.
HAWC event-level data is binned in a multi-dimensional array of event attributes. For a typical HAWC analysis, the array dimensions include the event direction
as a set of two sky coordinates (e.g., right ascension and declination)
and event quality and/or an energy metric.

The calculation of the logarithm of the likelihood is computationally convenient:
$$ \ln \mathcal{L} (\theta; \mathbf{N}_{obs}) = \sum_{i = 0}^{All Bins} \ln(f((N_{obs})_i | \theta )),$$
where $\theta$ is a set of model parameters, $(N_{obs})_i$ is the number of events in the ith bin, $\mathbf{N}_{obs}$ is the collection of all $(N_{obs})_i$,
and $f$ is the probability of observing $(N_{obs})_i$ given $\theta$.
To estimate the optimum set of parameters $\theta_0$, we find the $\theta$ that maximizes the value of $\ln \mathcal{L}$,
or equivalently we minimize the value of $-\ln \mathcal{L}$.
Estimates of the parameter uncertainties and the correlations between parameters can be calculated with well known methods \cite{UncertaintyEstimation}.

To compare hypotheses, the Likelihood Ratio test is used. A test statistic is defined as
$$ TS = 2 \ln \frac{ \mathcal{L} (\mbox{Alternative Hypothesis} ; \mathbf{N}_{obs})}  { \mathcal{L} (\mbox{Null Hypothesis}; \mathbf{N}_{obs})}.$$
When the null-hypothesis is true and in the case of nested models, $TS$ is approximately distributed as the $\chi^2$ distribution
with the number of degrees of freedom equal to the difference in the number of free parameters between the hypotheses (Wilks' theorem \cite{bib:Wilks}).
For example, the usual source discovery test has the null hypotheses as ``there is no source'' and the alternative hypothesis as ``there is a source with a flux normalization $X$''.
In this case, the difference between the number of free parameters is one (the flux normalization). When the null hypothesis is true, $TS$ is distributed as $\chi^2(DOF = 1)$,
and $\sqrt{TS}$ can be interpreted as the significance of the over-abundance of events in units of "Gaussian sigma".
Using the $TS$ formalism enables comparison of arbitrary hypotheses
with the full rigor of modern frequentist statistics. Note that non-nested models cannot use the Wilks' theorem approximation but there
are other methods available \cite{Non-nested}.

When the expected number of signal events in a bin is small (e.g., either due to a weak source, or a short observation time, or both),
Poisson probability functions must be used in the calculation of $\ln \mathcal{L}$ or the calculation
of the $TS$ will be biased and not easily interpreted.
In general, when using Poisson probability functions, $\ln \mathcal{L}$ is nonlinear with regard to the model parameters $\theta$.
Therefore, the minimization of the log-likelihood function is performed with a iterative, gradient search algorithm. We use the Minuit \cite{Minuit} minimization package.
If we assume a Gaussian probability function instead of Poisson, $\ln \mathcal{L}$ is linear with regard to the flux normalization parameter and can be solved algebraically.
This method is known as Gaussian weighting in the literature (e.g., \cite{GaussianWeighting}). Our framework incorporates a method that implements
the Gaussian weighting algorithm. The results of the Gaussian weighting algorithm are typically used as
a first guess for the minimization of the non-linear $\ln \mathcal{L}$ function.

The model parameters $\theta$ can be divided into physics parameters (the parameters of main interest) and nuisance parameters.
An example of a nuisance parameter for HAWC is the background normalization (the amount of hadronic events that pass event quality cuts) in a particular bin.
A profile maximum likelihood takes into account the idea of physics parameters and nuisance parameters.
Two nested loops are implemented where the nuisance parameters are estimated in the inner loop and the physics parameters
are estimated in the outer loop.
In this way, the outer loop can incorporate likelihood values for a particular physics model from several different instruments with the nuisance parameters
associated with each instrument effectively decoupled. This idea is part of the 3ML architecture, which is described in detail in \cite{Vianello_icrc}.
The 3ML idea can be thought of as standardizing the statistical language that is used both within a collaboration
and between collaborations -- it has wide applicability.

\section{Software Architecture}

In Figure 1 we show a diagram of the framework architecture and how it is used in the analysis of HAWC data.
The framework is composed of software objects grouped into five modules: 1.) Physics Model, 2.) Detector Response Model, 3.) Binned Data, 4.) Likelihood Calculation, and 5.) Minimizer.
The dependencies of the framework are \textsc{Root} \cite{Root} and HEALPix.
In the following sub-sections we briefly overview the highlights of the framework components.

\begin{figure}
\includegraphics[width=1.0\textwidth]{FrameworkDiagram1.png}
\caption{A diagram of the framework architecture. The code writting by the user is typically not complex. For example, software objects are instantiated and initialized,
then a \textit{GetLikelihood} method is called.}
\label{fig1}
\end{figure}

\subsection{Binned Event-level Data}

Event-level data is binned on the sky using HEALPix pixels, which are equal area pixels. A set of HEALPix pixels that covers the whole sky is called a HEALPix Map.
Event-level data are also binned in other dimensions.
For example, we currently typically use 10 event-quality bins. The events in a particular event-quality bin have a similar angular resolution and gamma-hadron resolution.
That is, the detector response is a function of the event-quality. Our current metric for event-quality is the number of hit photomultiplier tubes in the array,
which is strongly correlated to the energy deposited in the array.
One HEALPix map with event counts is stored for each event-quality bin.

An arbitrary number of bin dimensions beyond the two sky coordinates can be defined and implemented by the user.
Each additional bin, in a dimension other than the sky coordinates, is realized as an additional HEALPix Map.
For example, the data can be binned in both event-quality and in a metric which estimates the energy of the primary particle.
This energy metric is the energy deposited on the array adjusted for core location on the array and zenith angle.
Binning in an energy metric is important for detailed spectral analysis.
If we had 10 event quality bins and 10 energy bins, we would have a total of 100 HEALPix Maps. In this way, a HEALPix map is a "higher dimension bin".

In a typical analysis using our framework, only the bins that have a signal expectation according to the Physics Model and Detector Response Model
are summed over in the calculation of $\ln \mathcal{L}$.
For a Physics Model that is a point-source, for example, only sky pixels in a region of typically 5$^{\circ}$ radius around the hypothesized source location
have to be analyzed, since the point spread function of HAWC vanishes for larger angular separations.
Thus, we typically only keep a small fraction of the total number of sky-pixels in memory at any one time.
This significantly lowers our memory requirements (important if jobs are to be run in parallel) as well as our disk access time.
To handle partial sky regions in memory, we developed the \textit{SkyMap} class that inherits from the HEALPix base class.
The \textit{SkyMap} class stores a limited number of sky pixels where the boundaries of the pixel-set can be arbitrary polygonal or circular shapes.

As a fast and compressed file storage format for HEALPix maps, we implemented a class called \textit{MapTree} that depends on \textsc{Root} object classes.
Each pixel value is stored in a \textit{TTree} object with a single branch,
where the entry index number is identical to the HEALPix pixel index, and the \textit{TTree} is stored as a \textsc{Root} \textit{TFile}.
We have found that with this on-disk format, writing and reading times for a full set of sky pixels is faster than or comparable to
the FITS~\cite{bib:fits} on-disk format. More importantly, the writing and reading of partial sky regions are substantially faster.

Our time bin, the length of time for which data is collected to generate a single set of HEALPix maps, is an arbitrary time period.
The maximum length of the time bin is typically determined by the duration of a stable detector configuration, relevant in particular during the construction phase. To facilitate time-dependent analyzes like in~\cite{bib:agn}, the current minimum time interval is one sidereal day.

\subsection{Physics Model}

The physics model is implemented as a collection of one or more sources. We define two types of sources: 1.) A point source, and 2.) An extended source.
The point source has a location and a spectrum.
A general spectrum is defined by a table of differential flux values, binned in arbitrary energy intervals. For many analyzes, the spectrum can also be described by a simple function, for example a power law, in which case the spectrum is stored as a \textit{TF1} object \cite{Root}.
The extended source has a location, a spectrum, and a spatial morphology.
The spatial morphology interface is still under development but will eventually allow arbitrary tabular or functional inputs describing flux intensity and spectra as a function of location.

This implementation gives the user a tremendous amount of freedom in setting the physics model.
Note that each source can have an independent spectrum and that point sources and extended sources can be implemented in the same physics model.
Some examples of physics models that can be implemented and tested are
a closely spaced double source, a soft-spectrum source embedded in a hard-spectrum extended source, and the diffuse emission from the galactic plane.

\subsection{Detector Response Model}

If the time bin is an integer number of sidereal days, then the detector response is a function of the declination of the source but not the right ascension.
We define the detector response as a function of the declination and the source spectrum.
Given a source declination and arbitrary source spectrum
the code within the Detector Response Module calculates the expected point spread functions
(how events are spread on the sky due to detection limitations)
for each of the higher-dimension bins and the number of events in each of the higher-dimension bins.
This information is then used by code within the Likelihood Calculation Module to determine the expected number of signal events in each bin for a particular
physics model.

The description of the detector response is realized as a collection of histograms (\textit{TH1} objects
~\cite{Root}). These histograms are generated from simulation using
a particular spectrum. If the source spectrum called for in the physics model is different from the spectrum used to generate the original detector response, the detector
response is modified by re-weighting the energy distribution histograms for each higher dimension bin. This is a substantial computational short cut and is accurate as long
as the spectrum called for in the physics model is not too different from the spectrum used to generate the original detector response.

There is functionality such that the detector response can be parameterized and these parameters can be free in the minimization of $- \ln \mathcal{L}$.
This allows for the determination of the detector response directly from the data.
For example, if we assume that the spatial extent of the Crab Nebula is much smaller than our point spread function,
we can determine the point spread function that is most consistent with the collection of events with directions near the Crab~\cite{bib:icrccrab}.

\subsection{Likelihood Calculation}

To calculate $\ln \mathcal{L}$ we need to know the expected number of signal events and the expected number of background events for each bin.
The expected number of signal events in each relevant bin (bins that have zero expected source events are not included) is calculated
using the physics model parameters set by the user and the relevant detector response information.
The expected number of background events (i.e., hadron events that have passed our event quality selections) are normally not estimated from simulation but rather directly from the data.
The background is either estimated in a prior step that occurs when the data is first binned \cite{Andy'sContribution}, or the
background is estimated as a nuisance parameter as part of the minimization of $-\ln \mathcal{L}$.

Then to calculate the log likelihood using Poisson probability:
$$ \ln \mathcal{L} = \sum^{All Bins} (N_{obs} \ln N_{exp} - N_{exp} - \ln(N_{obs} !)).$$
Note that calculating $(\ln (x!))$ is numerically efficient (e.g., using Stirling's approximation).

A minimization of $-\ln \mathcal{L}$ can be performed with an arbitrary number of free parameters.
Two nested minimization loops are supported. This enables the calculation of a profile likelihood where nuisance parameters (e.g., the
expected number of background events in each bin) are estimated in the inner loop and the parameters associated with the physics model
are estimated in the outer loop. The calculation of such a profile likelihood is key to incorporating multiple instruments (e.g., HAWC and Fermi/LAT) in a joint-Likelihood fit.

\section{Conclusions}

The HAWC Collaboration has designed and coded a software framework for the statistical analysis of event-level data. The likelihood formalism is used.
With the addition of a small amount of code, a user can efficiently explore any number of complex hypotheses.
The framework is model-centric in that the user mainly focuses on the design of the physical model that she wants to confront with data.
The complexities of handling the binned data, the detector response, the calculation of $\ln \mathcal{L}$, and the minimizations are handled by
standardized software interfaces.

The framework was designed along side and is fully compatible with the 3ML architecture.
The 3ML architecture allows for the incorporation of multiple instruments in a joint-likelihood fit and is a powerful addition to the broader field
of astronomy and cosmic-ray physics.

\section*{Acknowledgments}
\footnotesize{
We acknowledge the support from: the US National Science Foundation (NSF);
the US Department of Energy Office of High-Energy Physics;
the Laboratory Directed Research and Development (LDRD) program of
Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnolog\'{\i}a (CONACyT),
Mexico (grants 260378, 55155, 105666, 122331, 132197, 167281);
Red de F\'{\i}sica de Altas Energ\'{\i}as, Mexico;
DGAPA-UNAM (grants IG100414-3, IN108713,  IN121309, IN115409, IN113612);
VIEP-BUAP (grant 161-EXC-2011);
the University of Wisconsin Alumni Research Foundation;
the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory;
the Luc Binette Foundation UNAM Postdoctoral Fellowship program.
}