text:
\section{Introduction}
Weak gravitational lensing by large-scale structure of the light from
distant galaxies is a powerful probe for constraining the cosmological
parameters and distinguishing between competing models of the Universe
\citep{Blandfordetal1991,Seitzetal1994,JainSeljak1997,Kaiser1998,Schneideretal1998,ZhangPetal2007}.
The first measurements of the correlations in the shapes of distant
background galaxies date back more than two decades
\citep{Baconetal2000,Kaiseretal2000,VanWaerbekeetal2000,Wittmanetal2000}. Since
then, cosmic shear observations have become ever more precise as the
coupling of techological advancements and algorithmic developments
have enabled us to conduct unprecedented deep optical imaging surveys
of the cosmos \mbox{KiDS}\footnote{{\tt kids.strw.leidenuniv.nl}},
DES\footnote{{\tt www.darkenergysurvey.org}} and HSC\footnote{{\tt
  hsc.mtk.nao.ac.jp/ssp/}}, with current state-of-the-art surveys now
mapping thousands of square degrees
\citep{Hildebrandtetal2017,Troxeletal2018,HSC2018,Hikageetal2019,
  Asgarietal2021}. By the end of the decade planned experiments like
\emph{Euclid}\footnote{{\tt www.cosmos.esa.int/web/euclid}} and the Rubin
Telescope\footnote{{\tt www.lsst.org}} \citep{Euclid2011,LSST2009}
will map volumes close to the entire physical volume of our observable
Universe. In order to make optimal use of these rich data sets we will
need to push forward our understanding and modelling of various
physical and measurement effects. In particular: accurate modelling of
the nonlinear evolution of large-scale structure, including the
baryonic physics effects; accurate modelling and correction of the
point-spread function of the telescope; correcting the bias in the
weak lensing shape estimation algorithms; and accounting for the
intrinsic alignments, to name but a few of the main systematics
\citep[see][for a more detailed discussion of these effects]{Schneider2006p3,Masseyetal2013,TroxelIshak2015}.

If the underlying matter density field were a Gaussian random field,
then all of the information in a weak lensing survey would be
contained in the shear two-point correlation function.  However,
physical effects like: the nonlinear growth of structure
\citep{Bernardeauetal2002}, the mapping between cosmic shear and
galaxy ellipticities \citep{Miralda-Escude1991}, and lensing beyond
the Born approximation
\citep{Hilbertetal2009,PrattenLewis2016,Fabbianetal2018}, all
introduce non-Gaussianity in the maps. Furthermore, the nonlinear
evolution also induces correlations in the convergence power spectrum
multipoles, which grow stronger on small scales. This means that the
information content of the second order statistics becomes saturated
after a given multipole
\citep{Satoetal2011,Hilbertetal2012,Kayoetal2013,Marianetal2013,Byunetal2017}. Thus
in order to capture all of the cosmological information available in
lensing surveys one must look to the higher order statistics of the
shear field
\citep{Schneideretal1998,Bernardeauetal2002,SchneiderLombardi2003}. Furthermore,
owing to the different ways in which the cosmological parameters and
nuisance parameters depend on the higher-order statistics, the
inclusion of such measurements brings with it the further virtue of
being able to break parameter degeneracies, e.g. by combining second
and third order statistics
\citep{KilbingerSchneider2005,Sembolonietal2011,Fuetal2014}, or by
incorporating the information found in the statistical properties of
the peaks in the shear field \citep{Marianetal2013,Kacprzaketal2016}.

A powerful method to disentangle systematic effects from cosmic shear
signals is the E/B decomposition
\citep{Crittendenetal2001,Schneideretal2002a}. At leading order, pure
weak lensing signals are sourced by a scalar lensing potential, which
means that their deflection fields are curl free.  Equivalently, the
ring-averaged cross component of the shear is expected to be zero (the
B mode), while the tangential one contains all the lensing signal (the
E mode). Thus B modes enable a robust test for the presence of
systematic errors. One method to take advantage of this E/B
decomposition is the so-called `aperture mass statistics'
\citep{Kaiser1995,Schneider1996,Schneideretal1998}. `Aperture mass'
(${\mathcal M}_{\rm ap}$) and `Map-Cross' $(\mathcal{M}_{\times})$ are obtained by
convolving the tangential and cross shear with an isotropic filter
function. Therefore by construction they are E/B-decomposed. Taking
the second moment leads to the variance of aperture mass, the third to
the skewness, the fourth to the kurtosis, etc.

The standard approach for measuring the aperture mass statistics in
data utilises the fact that, for the flat sky, any $n$-point moment
can be expressed in terms of integrals over the $n$-point shear
correlation functions, modulo a kernel function
\citep{Schneideretal2002a,Jarvisetal2004}. The reason for adopting
this strategy stems from the fact that the correlation functions can
reliably be estimated in the presence of a nontrivial survey
mask. However, for these estimators to be accurate and E/B decomposed,
one requires three conditions to be satisfied: (i) the $\xi_+$/$\xi_-$
correlations need to be measured down to zero separation; (ii) they
also need to be measured up to a maximum angular scale, set by the
exact form of the aperture mass filter and its angular scale; (iii)
the angular bins must be sufficiently fine for the discretisation of
the integrals to be reliable
\citep{KilbingerSchneider2005,Fuetal2014}. Owing to galaxy image
blending, signal-to-noise issues and the finite size of the survey,
the lower bound is never possible and the upper bound means that
biases can occur due to edge effects. In addition, while the mean
estimate is unbiased, the covariance matrix does require one to
carefully account for the mask
\citep{Schneideretal2002b,Friedrichetal2016}. More recent developments
that also make use of the shear correlation functions, while
circumventing the issues of E/B leakage on small scales are the ring
statistics and COSEBIs
\citep{Schneideretal2007,Schneideretal2010}. While those approaches
can in principle be extended to higher order statistics, the
estimation of the $n$-point correlation functions turns out to be
notoriously time consuming \citep{Schneideretal2005a,
  Jarvisetal2003}. Further methods to extract non-Gaussian information
from the aperture mass look at its probability density function as a
whole \citep{Bernardeauetal2000, Munshietal2004, Barthelemyetal2020}
or at the distribution of its signal-to-noise peaks
\citep{Marianetal2012, Heydenreichetal2020, Martinetetal2021}.

In \citet{Porthetal2020} we took a different approach and explored a
computationally efficient (accelerated) implementation of the original
direct estimator of the aperture mass dispersion
\citep{Schneider1998}. Rather than measuring the correlation functions
of the shear polar, in this formulation one instead directly measures
cumulants of $\mathcal{M}_{\rm ap}$ on a set of apertures and then
uses an optimised weighting scheme to average the estimates, along
with a restriction on the types of apertures that are acceptable. The
present work extends our previous investigation in a number of
important ways. First, we construct accelerated direct estimators for
the higher-order aperture mass moments, including the skewness,
kurtosis, etc. Second, we also develop further the multiscale
aperture moments \citep{Jarvisetal2003,Schneideretal2005a}. These two
improvements enable us to better trace the full, harmonic mode,
configuration dependence of the convergence polyspectra.

This paper is organised as follows: In \Sect{sec:Theory} we introduce
key concepts of weak lensing, define the aperture mass and show how
its connected cumulants are related to the convergence polyspectra. In \Sect{sec:Estimators} we
revisit the direct estimators for higher order aperture mass measures
and construct suitable bases, in which each statistic can be computed
in linear time complexity. After investigating the variance of the
direct estimators, we give details of our updated algorithm used to
perform the measurements. In \Sect{sec:Measurements} we
empirically verify the linear scaling and the measurements of our implementation of the direct estimator on Gaussian mocks. In \Sect{sec:DetectionSignificance} we then apply the estimator to the SLICS simulation suite in
order to assess up to which order one can expect to extract
information from the aperture mass statistics on a KiDS-1000 like
survey. Finally, in \Sect{sec:Conclusions} we summarise our
findings, conclude and discuss future work.

\section{Higher order aperture mass measures for cosmic shear}\label{sec:Theory}

\subsection{Weak gravitational lensing and aperture mass}\label{ssec:apmass}

In this paper we are mainly concerned with the weak lensing of distant
background (source) galaxy shapes by the intervening large-scale
structure
\citep[for detailed reviews of the topic
  see][]{BartelmannSchneider2001,Schneider2006p1,Schneider2006p3,Dodelson2003,Dodelson2017,Kilbinger2015,Mandelbaum2018}.
The two fundamental quantities describing this mapping from true to
observed galaxy images are the convergence $\kappa$ and the complex
shear $\gamma=\gamma_1+i\gamma_2$, which, assuming a metric theory of
gravity, are all derived from an underlying scalar lensing
potential. In a galaxy survey the effective convergence at angular
position $\bm \theta$ and radial comoving distance $\chi$ can be
connected to the density contrast $\delta(\chi\bm\theta,\chi)$ through:
\begin{align}
  \kappa({\bm\theta}) & =  \frac{3}{2}\Omega_{\rm m,0}\left(\frac{H_0}{c}\right)^2
  \int_0^{\chi_{\rm H}} {\rm d}\chi' \frac{\chi'}{a(\chi')} g(\chi') \delta(\chi'{\bm\theta},\chi') \ ,\label{eq:kappaeff}
\end{align}
where $\Omega_{\rm m,0}$ is the total matter density, $H_0$ denotes
the Hubble constant, $a$ is the scale factor, $c$ is the speed of
light, $\chi_H$ is the comoving distance to the horizon and $g(\chi)$
is a weight function related to the normalized redshift distribution
${\rm d} n(z)/{\rm d} z$ of the source galaxies as
\begin{align}
  g(\chi') &\equiv \int_{z(\chi')}^{z_{\rm H}} {\rm d} z \
   \frac{{\rm d} n(z)}{{\rm d} z}
  \frac{\left[\chi(z)-\chi'\right]}{\chi(z)} \ . \label{eq:gchi}
\end{align}
Aperture mass was developed by \citet{Schneider1996} as a technique to
estimate projected mass overdensities enclosed within a circular
region:
\begin{equation}
{\mathcal M}_{\rm ap}({\bm\theta}_0;{\vartheta}) = \int_{{\mathbb{R}^2}} {{\rm d}^2\bm\theta}_1 \kappa({\bm\theta}_1)
U(\left|{\bm\theta}_1-{\bm\theta}_0\right|;{\vartheta}) \ , \label{eq:MapKappa}
\end{equation}
where $U$ is a compensated filter function.  In the flat sky limit the
(cross) aperture mass can be expressed in terms of a related
circularly symmetric filter function $Q(U)$ and the complex shear
field $\gamma$ in its E/B-decomposed basis:
\begin{equation}
\mathcal{M}_{{\rm ap}/\times} ({\bm\theta}_0;{\vartheta}) \equiv \int_{{\mathbb{R}^2}} {{\rm d}^2\bm\theta}_1 \gamma_{{\rm t}/\times}({\bm\theta}_1;{\bm\theta}_0)
Q(|{\bm\theta}_{1}-{\bm\theta}_0|;{\vartheta}) \ , \label{eq:MapShear}
\end{equation}
where the tangential and cross components of the shear field at
position ${\bm\theta}+{\bm\theta}_0$ with respect to the aperture center
${\bm\theta}_0$ are defined as \citep{BartelmannSchneider2001}:
\begin{align}
\gamma_{\rm t}({\bm\theta};{\bm\theta}_0) + i \gamma_{\times}({\bm\theta};{\bm\theta}_0) & \equiv  - \gamma(\bm\theta+{\bm\theta}_0){\rm e}^{-2i\phi}  \label{eq:gt1} \ ,
\end{align}
in which $\phi$ denotes the polar angle associated with the vector
${\bm\theta}$.  In the absence of systematic errors (B-modes) in the
lensing data, map-cross should vanish \citep{Schneideretal2002a}.

For this work we will make use of the polynomial filter function
introduced by \citet{Schneideretal1998}:
\begin{equation}
Q(\theta;{\vartheta}) = \frac{6}{\pi {\vartheta}^2}
\left(\frac{\theta}{{\vartheta}}\right)^2
\left[1-\left(\frac{\theta}{{\vartheta}}\right)^2\right]
\mathcal{H}({\vartheta} - \theta) \label{eq:schneiderQ} \ ,
\end{equation}
where ${\vartheta}$ is the characteristic scale of the filter and
${\mathcal H}(x)$ is the Heaviside function, which guarantees
that the filter function has compact support.

\subsection{A hierarchy of aperture mass measures}

One may construct moments of the aperture mass field, and this gives
rise to the so called aperture mass statistics. At the two-point level
this gives us the variance $\MapStatEns{2}_{\rm c}(\vartheta_1)$ and
at the three-point, the skewness $\MapStatEns{3}_{\rm
  c}(\vartheta_1)$, etc., where the subscript $\rm c$ stands for the
connected cumulant obtained from the moments
\citep{ScoccimarroFrieman1996a}.  Owing to the fact that the aperture
mass is a convolution of the convergence field with a filter function,
it is possible to rewrite these moments in terms of their Fourier
space counterparts, that is the convergence spectra. For example for
the variance and skewness we have:
\begin{align}\label{eq:Map2Kappa}
  \MapStatEns{2}_{\rm c}(\vartheta) & = \int \dellnorm{1} C_{\kappa,2}(\ellvec{1})
  \ \left|\widetilde{U}_{\vartheta}(\ellvec{1})\right|^2    \ ; \\
  %
  \label{eq:Map3Kappa}
  \MapStatEns{3}_{\rm c}(\vartheta)
  &=
  \int \dellnorm{1} \cdots \int \dellnorm{3} \ (2\pi)^2\delta^D\left(\sum_{i=1}^3\ellvec{i}\right)
  \nonumber \\ &\hspace{-0.5cm} \times
  C_\kappa(\ellvec{1}, \cdots, \ellvec{3}) \
  \widetilde{U}_{\vartheta}(\ellvec{1})\widetilde{U}_{\vartheta}(\ellvec{2}) \widetilde{U}_{\vartheta}(\ellvec{3}) \  \ ,
\end{align}
where $\widetilde{U}_{\vartheta_i}$ denotes the Fourier transform of
the aperture mass filter function $U(\theta; \vartheta_i)$ and
$C_\kappa(\ellvec{1})$ denotes the convergence power spectrum, and
$C_\kappa(\ellvec{1}, \ellvec{2},\ellvec{3})$ the convergence
bispectrum. These spectra can formally be defined:
\begin{align}
    \left\langle \tilde{\kappa}(\ellvec{1}) \tilde{\kappa}(\ellvec{2}) \right\rangle_{\rm c}
    & =
    (2\pi)^3 \delta^D\left(\ellvec{1}+\ellvec{2}\right)  \
    C_\kappa(\ellvec{1}) \ ; \\
    \left\langle \tilde{\kappa}(\ellvec{1}) \tilde{\kappa}(\ellvec{2})\tilde{\kappa}(\ellvec{3}) \right\rangle_{\rm c}
    & =
    (2\pi)^3 \delta^D\left(\sum_{i=1}^3\ellvec{i}\right)  \
    C_\kappa(\ellvec{1},\ellvec{2},\ellvec{3})\ .
\end{align}
This of course can be generalised to $n$-point aperture mass moments:
\begin{align}\label{eq:MapnKappa}
    \MapStatEns{n}_{\rm c}(\vartheta)
    &=
    \int \dellnorm{1} \cdots \int \dellnorm{n} \ (2\pi)^2\delta^D\left(\sum_{i=1}^n\ellvec{i}\right)
    \nonumber \\ &\hspace{-0.5cm} \times
    C_\kappa(\ellvec{1}, \cdots, \ellvec{n}) \ \widetilde{U}_{\vartheta}(\ellvec{1}) \cdots
    \widetilde{U}_{\vartheta}(\ellvec{n}) \  \ ,
\end{align}
where the $n$-point convergence spectrum is defined:
\begin{align}
  \left\langle \tilde{\kappa}(\ellvec{1}) \dots
  \tilde{\kappa}(\ellvec{n}) \right\rangle_{\rm c}
    & =
    (2\pi)^3 \delta^D\left(\sum_{i=1}^n\ellvec{i}\right)  \
    C_\kappa(\ellvec{1},\dots,\ellvec{n})\ .
\end{align}

It is worth noting that due to the fact that $\widetilde{U}$ is a
sharply peaked filter function in Fourier space, the aperture mass
moment on a given scale only carries information about a specific range
of wavemodes $\vec{\ell}$ from the underlying polyspectrum. In order to
extract more of the information that is available one needs to compute
\Eqn{eq:MapnKappa} for a large set of aperture radii
\citep{Schneideretal2005a}.

\subsection{Multiscale aperture mass moments and their correlators}\label{sec:TheoryMultiscale}

Even if one considers a wide range of aperture radii there will be
certain wavemode configurations of the polyspectra that are
suppressed when compared with other configurations. This may result in
a loss of sensitivity to certain physical effects that are only
manifest in the higher-order polyspectra, such as those induced by
modifications of gravity or primordial non-Gaussianities. In order to
combat this one can further generalise the aperture mass moments in
several ways. First, if we choose different scales for the aperture
mass filter function, then we get the multiscale aperture mass moments.
For the $n$-point multiscale aperture mass moment this can be written:
\begin{align}\label{eq:MapnKappaMS}
    \MapStatEns{n}_{\rm c}(\vartheta_1, \cdots, \vartheta_n)
    &=
    \int \dellnorm{1} \cdots \int \dellnorm{n} \ (2\pi)^2\delta^D\left(\sum_{i=1}^n\ellvec{i}\right)
    \nonumber \\ &\hspace{-0.5cm} \times
    C_\kappa(\ellvec{1}, \cdots, \ellvec{n}) \
    \widetilde{U}_{\vartheta_1}(\ellvec{1}) \cdots \widetilde{U}_{\vartheta_n}(\ellvec{n}) \  \ .
\end{align}

Second, if we correlate a set of apertures at different spatial
positions in the sky, then one can define the multiscale aperture
mass moment correlators \citep{Szapudietal1997,Munshietal2003}.  There
are two special cases where this approach can be applied, the first is
the case where the separation of the aperture is directed
perpendicular to the line of sight. The second case is where the
apertures are placed along the same line of sight, but where different
tomographic bins of source galaxies are used to estimate the aperture
mass. The former case measures the correlation of the cumulants on the
same redshift slice, but at different angular positions. The latter
case corresponds to correlating aperture measures in different surveys
with overlapping footprints, or between photometric redshift bins
within the same survey. As the aperture mass filter carries most of its weight in a compact region surrounding the aperture center one expects the signal to fall off rapidly for
aperture separations that exceed beyond a few times the aperture
radius. Generalizing the result of \citep{Schneideretal1998} we can
formally write this as follows:
\begin{align}\label{eq:MapStatCorr}
    \left\langle
    \mathcal{M}_{\text{ap}}^n
    \mathcal{M}_{\text{ap}}^m
    \right\rangle_{\rm c}(\vartheta_1, \cdots, \vartheta_n, \vartheta'_1,
    \cdots, \vartheta'_m; \overrightarrow{\Delta})
    &=
    \nonumber \\ &\hspace{-5cm}
    \int \dellnorm{1} \cdots \dellnorm{n+m} \ (2\pi)^2 \delta^D
    \left(\sum_{j=1}^{n+m}\vec{\ell_j}\right) \ C_\kappa(\vec{\ell}_1, \cdots, \vec{\ell}_{n+m}) \
    \nonumber \\ &\hspace{-4.5cm} \times \hspace{.2cm}
    \tilde{U}_{\vartheta_1}(\vec{\ell}_1) \cdots \tilde{U}_{\vartheta'_m}(\vec{\ell}_{n+m})
    \ e^{i\overrightarrow{\Delta} \sum_{j=1}^m\vec{\ell}_{n+j}}\ ,
\end{align}
where $\vec{\Delta}$ is a separation vector. Note that for zero
separation we recover the $(m+n)$th cumulant. In addition, we can
assess the impact of the exponential factor by evaluating the two
point cross-correlation coefficients $r_{mn}$, which are defined
in a similar way to those in \citep{Munshietal2005}:
\begin{align}\label{eq:twocrosscorrcoeff}
r_{mn}(\Delta) \equiv \frac{\left\langle X^m X^n \right\rangle_c(\Delta)}{\left\langle X^{m+n} \right\rangle_c}  \ ,
\end{align}
where for our case $X^m=\mathcal{M}_{\text{ap}}^m$.  In this work,
however, we do not consider the cosmological information contained in
\Eqn{eq:MapStatCorr}, but instead use it to assess how fast the
$r_{mn}$ converge to unity - this can be seen as a proxy for how
densely apertures need to be sampled within a survey footprint to
retrieve all available signal.

\section{Estimators for higher order aperture mass statistics}\label{sec:Estimators}

\subsection{Direct estimators for the aperture mass moments and
  their evaluation in linear order time}

In this subsection we concern ourselves with estimators for higher
order aperture mass statistics that mimic the original theoretical
expressions \Eqn{eq:MapnKappa} more closely. At first, let us
investigate the special case of all the radii being equal.

Consider an aperture of angular radius $\vartheta$, centred on the
position ${\bm\theta}_0$. The aperture contains $N$
galaxies\footnote{Strictly speaking, we select galaxies within the
support of the $Q$ filter function of that aperture. For the filter
functions we use in this work the support is always concentric around
the aperture center and linearly scaling with aperture
radius. Therefore we will continue referring to $N$ as the number of
galaxies per aperture.} with positions ${\bm\theta}_i$, complex
ellipticities $e_i$ and weights $w_i$. Then, for a single aperture,
one can write down an estimator for the $n$th order aperture mass
statistic \Eqn{eq:MapnKappa} as \citep{Schneideretal1998,
  Munshietal2003}
\begin{align}\label{eq:MapnEst}
  \MapStatEst{n} = (\pi \vartheta ^2)^n\frac{\sum_{(i_1, ..., i_n)^N} w_{i_1}Q_{i_1}e_{t,i_1}
    \cdots w_{i_n}Q_{i_n}e_{t,i_n}}{\sum_{(i_1, ..., i_n)^N}w_{i_1} \cdots w_{i_n}} \ ,
\end{align}
where we defined the shorthand notation
\begin{equation} \label{eq:ShorthandSum}
\sum_{(i_1,\dots,i_n)^N} \equiv \sum_{i_1=1}^N \sum_{i_2\neq i_1}^N
\dots \sum_{i_n\neq i_{n-1}\neq\dots\neq i_{1}}^N\ .
\end{equation}
In certain cases we might use further abbreviations, meaning that $(i_1,...,i_n)^N \equiv (i_1,...,i_n) \equiv \ \neq$. On applying the above estimator to the case of $n=2$, one can easily
show that that this estimator is unbiased after averaging over the
intrinsic ellipticity distribution, the galaxy positions within the
aperture, and finally over cosmological ensembles
\citep{Schneideretal1998, Porthetal2020}.

If we were to apply the above estimator given by \Eqn{eq:MapnEst} to
determine the hierarchy of aperture mass moments, then this naive
implementation would appear to result in an estimator that requires of
the order $N^n$ operations to compute. However, following our earlier
work \citep{Porthetal2020}, one can complete the sums to transform the
estimators into sums and products of linear order terms. In
Appendix~\ref{app:estimators1} we explicitly show, using elementary
means, how one can compute the skewness ($\MapStatEst{3}$) and
kurtosis ($\MapStatEst{4}$) using linear sums. The results for
second, third and fourth orders are:
\begin{align}
\widehat{M_{\rm ap}} &= \MsEst{1} \ ; \\
\MapStatEst{2} &= \frac{\MsEst{1}^{2} - \MsEst{2} }{1-\SEst{2}}  \ ;\\
\MapStatEst{3} &= \frac{\MsEst{1}^{3} - 3\MsEst{2}\MsEst{1} + 2\MsEst{3}}{1-3\SEst{2}+2\SEst{3}} \ ;\\
\MapStatEst{4} &= \frac{\MsEst{1}^{4} - 6\MsEst{2}\MsEst{1}^{2} + 3\MsEst{2}^{2} +8\MsEst{3}\MsEst{1} -6\MsEst{4}}
           {1-6\SEst{2}+3\SEstPow{2}{2}+8\SEst{3}-6\SEst{4}}\ ,
\end{align}
where we have introduced two additional quantities: $\SEst{m}$ and
$\MsEst{m}$, which are defined:
\begin{align}
  \MsEst{m} &\equiv \SEst{m} \ (\pi \vartheta^2)^m \
  \frac{\sum_{i=1}^N w_i^m Q_i^m e_{t,i}^m}{\sum_{i=1}^N w_i^m} \ ;
  \label{eq:MapnEstPowerSumsA}\\
  \SEst{m} &\equiv \frac{\sum_{i=1}^N w_i^m}{\left(\sum_{i=1}^N w_i\right)^m} \ \ ;
  \label{eq:MapnEstPowerSumsB}
\end{align}

Applying the elementary approach described in
Appendix~\ref{app:estimators1} beyond fourth order rapidly becomes
cumbersome, to say the least. We have therefore developed an analytic
method for generation of the $n$th order estimator decomposed into
linear sums. This follows from noting that the sum in \Eqn{eq:MapnEst}
runs over unequal indices and that one can express any statistic
$\MapStatEst{n}$ as a sum of the power sums
\Eqn{eq:MapnEstPowerSumsA}, where the coefficients preceding each term
are determined with the help of the complete Bell polynomials
$B_n$. Hence, for the general $n$th order estimate one has:
\begin{align}\label{eq:MapnEstBell}
\MapStatEst{n} =
\frac{B_n\left(-\MsEst{1}, -\MsEst{2}, -2\MsEst{3}, ..., -(n-1)!\MsEst{n}\right)}
     {B_n\left(-\SEst{1}, -\SEst{2}, -2\SEst{3}, ..., -(n-1)!\SEst{n}\right)} \ .
\end{align}
For full details of this derivation we refer the reader to
Appendix~\ref{app:estimators2}. Here we only note that each argument
that goes into \Eqn{eq:MapnEstBell} is a single sum over the galaxies
in the aperture and is therefore independent of the order of the
statistic. Using this formalism, we extend our decomposition to 5th
and 6th order, giving us:
\begin{align}
\MapStatEst{5} &= \frac{1}{N_5}
\left[ \MsEst{1}^5 - 10\MsEst{2}\MsEst{1}^{3} + 15\MsEst{2}^{2}\MsEst{1} +
  20\MsEst{3}\MsEst{1}^{2} \right. \nn \\
  & \left. \hspace{0.2cm} -20 \MsEst{3}\MsEst{2} -30\MsEst{4}\MsEst{1} + 24\MsEst{5} \right]\ ; \\
\MapStatEst{6} &= \frac{1}{N_6}
\left[\MsEst{1}^{6}-15\MsEst{2}\MsEst{1}^{4}+45\MsEst{2}^{2}\MsEst{1}^{2}-15\MsEst{2}^{3}
  +40\MsEst{3}^{2}\right. \nn \\
  & \hspace{0.2cm} -90\MsEst{4}\MsEst{1}^{2}+40\MsEst{3}\MsEst{1}^{3}-120\MsEst{3}\MsEst{2}\MsEst{1} \nn \\
  & \hspace{0.2cm} \left. +90\MsEst{4}\MsEst{2}+144\MsEst{5}\MsEst{1}-120\MsEst{6}\right]
  \hspace{0.2cm} \ ,
\end{align}
where
\begin{align}
  N_5 =& 1-10\SEst{2}+15\SEstPow{2}{2}+20\SEst{3}-20\SEst{3}\SEst{2}-30\SEst{4}+24\SEst{5} \ ;\\
  N_6 =&  1-15\SEst{2}+45\SEstPow{2}{2}-15\SEstPow{2}{3}+40\SEst{3}-120\SEst{3}\SEst{2} \nn \\
  & +40\SEstPow{3}{2}-90\SEst{4}+90\SEst{4}\SEst{2}+144\SEst{5}-120\SEst{6}\ .
\end{align}

\subsection{Direct estimators for the multiscale aperture mass moments}

In complete analogy we can write down an unbiased direct estimator for
the full multiscale aperture mass moments of \Eqn{eq:MapnKappaMS}:
\begin{align}\label{eq:MapnREst_h1}
\MapStatEst{n}(\vartheta_1, ..., \vartheta_n) &= (\pi \vartheta_1 ^2) \cdots (\pi \vartheta_n ^2)
\nonumber \\
&\hspace{-.8cm}\times \frac{\sum_{(i_1, ..., i_n)} w_{i_1}Q_{\vartheta_1,i_1}e_{t,i_1}
  \cdots w_{i_n}Q_{\vartheta_n,i_n}e_{t,i_n}}{\sum_{(i_1, ..., i_n)^N}w_{i_1} \cdots w_{i_n}} \ ,
\end{align}
where each index runs through all the galaxies within the aperture of
the largest radius. In this case the power sums of
\Eqns{eq:MapnEstPowerSumsA}{eq:MapnEstPowerSumsB} do not form a
sufficient basis to express these estimators, but we are still able to
write down the estimators from elements within the sets
\begin{equation}
X_n \equiv \left\{ X_{(\mathscr{s}_1, ..., \mathscr{s}_n)}^{(m)} \  \left|  \  \mathscr{s}_i \in \{0,1\} \ ,
\  \sum_{i=1}^n\mathscr{s}_i = m \leq n \ \right. \right\} \ ,
\end{equation}
where $X \in \{ M_{\rm s}, S \}$ and the corresponding elements
constitute of multivariate power sums being defined as
\begin{align}\label{eq:Msvar}
\MsEstPoly{m}{\mathscr{s}_1, ..., \mathscr{s}_n}
&\equiv
\left(\prod_{k=1}^n \left(\pi \vartheta^2_k\right)^{\mathscr{s}_k}\right)
\sum_{i=1}^{N(\mathscr{s})} w_i^m \prod_{j=1}^n \left[ e_{t,i}Q_{\vartheta_j,i} \right]^{\mathscr{s}_j} \ ,
\nonumber \\
\SEstPoly{m}{\mathscr{s}_1, ..., \mathscr{s}_n} &\equiv
\sum_{i=1}^{N(\mathscr{s})} w_i^m  \ ,
\end{align}
where $N(\mathscr{s})$ denotes the number of galaxies within the aperture
of the smallest radius for which $\mathscr{s}_i$ is not zero. Despite the
more complicated looking form compared to the equal radius case these
estimators can also be computed in $\mathcal{O}(N)$ time using the
$\left|\hat{X}_n\right| = 2^n-1$ distinct multivariate power sums
\Eqn{eq:Msvar} and summing over various partitions $P$ of the set
$\{1, \cdots, n \}$:
\begin{align}\label{eq:MapnREst}
  \MapStatEst{n}&(\vartheta_1, ..., \vartheta_n) = \nonumber\\
    &\frac{\sum_{m=1}^n \sum_{\pi \in P_{n,m}} (-1)^m \prod_{i=1}^{m} (n_i-1)! \
    \MsEstPoly{n_i}{\mathscr{s}_1(\pi_i), \cdots, \mathscr{s}_n(\pi_i)}
  }{\sum_{m=1}^n \sum_{\pi \in P_{n,m}} (-1)^m \prod_{i=1}^{m} (n_i-1)! \
    \SEstPoly{n_i}{\mathscr{s}_1(\pi_i), \cdots, \mathscr{s}_n(\pi_i)}} \ .
\end{align}
In this expression the combination of the two outer sums run through
each partition $\pi$ that consists of $m$ blocks and the $\alpha(\pi_i)$
denote the value of the $\alpha$ as evaluated from the $i$th block of
the partition. For a motivation of this equation and explicit
expressions we again refer to \App{app:estimators2}.

\subsection{Estimators applied to a large survey}

In order to estimate any aperture statistics
\begin{equation} \mathfrak{M}\in \left\{\left<M_{\rm ap}^{2}\right>(\vartheta_1),\ \left<M_{\rm
    ap}^{3}\right>(\vartheta_1),\dots \right\}\end{equation}
on a contiguous survey field one can simply place an ensemble of
apertures on the field and compute their weighted means
\begin{align}\label{eq:MapStatEns}
  \hat{\mathfrak{M}} = \frac{\sum_{\rm ap} w_{\rm ap}
    \hat{\mathfrak{M}}_{\rm ap}}{\sum_{\rm ap} w_{\rm ap}} \ ,
\end{align}
where the weights $w_i$ should be chosen to minimize the variance of
the estimator. Owing to the linearity of \Eqn{eq:MapStatEns}, if the
estimator of a single aperture is unbiased, then so is
\Eqn{eq:MapStatEns}. Thus including more apertures will increase the
signal-to-noise of the ensemble estimator.

\subsection{Variance of the direct estimators}\label{sec:DirectCov}

In order to understand how to weight the apertures we need to obtain
expressions for the variance of the moment estimators. On generalizing
the prescriptions outlined in \citep{Schneideretal1998,
  Munshietal2003} to include the shear weights, as in
\citet{Porthetal2020}, one can work out expressions for the variance
of the higher-order direct estimators \Eqn{eq:MapnEst} for a given
aperture. For the explicit derivation of for the variance of the third order statistic see
Appendix D in the online supplementary material. From this analysis we see that the general
expression can be written:
\begin{align}\label{eq:MapStatCov}
\sigma^2\left[\MapStatEst{n}\right]
&=
\sum_{\ell=0}^n  \sum_{m=\ell}^n \frac{\sum_{\neq} w_{i_1}^2\cdots w_{i_m}^2 w_{i_{m+1}}
  \cdots w_{i_n} w_{j_{m+1}}\cdots w_{j_n}}{\left( \sum_{\neq} w_{i_1}\cdots w_{i_n}\right)^2}
\nonumber \\
&\times C(n,\ell,m) M_{g,2}^{\ell}\sigsq{\ell}  \MsMapStatEns{m-\ell}{2(n-m)} \ ,
\end{align}
where the sum over the galaxy weights can again be decomposed as sums
of (bivariate) power sums and the multiplicities are given by
\begin{align}
C(n, \ell, m) & = \frac{1}{\ell!(m-\ell)!} \left(\frac{n!}{(n-m)!}\right)^2 \ .
\end{align}
For a discussion of the origin for the multiplicity factor $C(n, \ell,
m)$, as well as a motivation of \Eqn{eq:MapStatCov} and some of its limits we refer the reader to \App{app:covariance}; in particular we obtain for the shape noise dominated limit
\begin{align}
    \sigma^2\left[\MapStatEst{n}\right]
    \approx
    n!\frac{\sum_{\neq} w_{i_1}^2\cdots w_{i_n}^2}{\left( \sum_{\neq} w_{i_1}\cdots w_{i_n}\right)^2}  \ \left(\frac{\sigma^2_\epsilon}{2}\right)^n \MgStat{n} \  . \label{eq:MapStatCovShotLimit}
\end{align}
The above formula gives the variance per aperture, thus for the
estimator over the full survey field, \Eqn{eq:MapStatEns}, the
variance can be written down as:
\begin{align}\label{eq:MapStatEnsCov}
  \sigma^2\left[\hat{\mathfrak{M}}\right] &= {\rm cov}\left( \frac{\sum_i w_i \hat{\mathfrak{M}}_i}{\sum_i w_i},
  \frac{\sum_j w_j \hat{\mathfrak{M}}_j}{\sum_j w_j} \right)
    \nonumber \\ &=
    S_2 \ \sigma^2\left[{\hat{\mathfrak{M}}_{\rm ap}}\right] +
    \frac{\sigma^2\left[{\hat{\mathfrak{M}}_{\rm ap}}\right]}{\left(\sum_i w_i\right)^2}
    \sum_{i \neq j} w_i w_j \rho (\hat{\mathfrak{M}}_i, \hat{\mathfrak{M}}_j)\ ,
\end{align}
\begin{figure}
    \centering
    \centering{
      \includegraphics[width=0.45\columnwidth]{FIGS/SpatialCrossCorrExampleRho.pdf}
      \includegraphics[width=0.45\columnwidth]{FIGS/SpatialCrossCorrExampleSN.pdf}
      }
      \caption{Example configuration of the correlation coefficient $\rho$ (left) and its effect on the signal contained in a survey field as predicted from \Eqn{eq:MapStatEnsCovb} (right).}
\label{fig:SpatialCorrCoeffExample}
\end{figure}
where in the above we have defined the cross-correlation coefficient
between apertures whose centres are at position $\bm\theta_i$ and
$\bm\theta_j$ to be:
\begin{equation} \label{eq:twocrosscorrcoeff2}
\rho(\mathfrak{M}_i, \mathfrak{M}_j)
\equiv \frac{\left<\mathfrak{M}_i, \mathfrak{M}_j\right>}
       {\sqrt{\left<\mathfrak{M}_i, \mathfrak{M}_i\right>\left<\mathfrak{M}_j, \mathfrak{M}_j\right>}} \ .
\end{equation}
Note that for the case of well separated apertures, the
cross-correlation coefficient will vanish and only the first summand
needs to be taken into account, which for unity weights gives the
familiar $1/N_{\rm ap}$ scaling of the variance. If the apertures are
oversampled, this assumption is no longer valid and the term involving
$\rho$ must be included. Owing to the fact that $\rho$ should only
depend on the relative spatial distance $\Delta$ between the aperture
centres, we can rewrite \eqref{eq:MapStatEnsCov} as a weighted sum over all possible distances between aperture center pairs:
\begin{align}
    \sigma^2\left[\hat{\mathfrak{M}}\right] = \
    &S_2 \ \sigma^2\left[{\hat{\mathfrak{M}}_{\rm {ap}}}\right] +
    \nonumber \\
    &\frac{\sigma^2\left[{\hat{\mathfrak{M}}_{\rm {ap}}}\right]}{\left(\sum_i w_i\right)^2}
    \sum_{b \in \rm bins} \left(\sum_{i,j \in \mathcal{I}_b} w_i w_j \right)
    \rho(\hat{\mathfrak{M}}, \Delta_b)
     \\ \approx
    &\frac{\sigma^2\left[{\hat{\mathfrak{M}}}_{\mathrm{ap}}\right]}{N_{\mathrm{ap}}} + 2\pi \  \frac{\sigma^2\left[{\hat{\mathfrak{M}}}_{\mathrm{ap}}\right]}{A_{\mathrm{survey}}}  \int_{R_{\mathrm{ap}}/\alpha}^\infty {\rm d}\Delta \ \Delta \  \rho (\hat{\mathfrak{M}}, \Delta)
    \ , \label{eq:MapStatEnsCovb}
\end{align}
where in the first step the bins are defined as a partition of the reals, and $\mathcal{I}_b
\equiv \{ i,j | \Delta(i,j) \in b \}$ collects all the aperture center
pairs falling into bin $b$. For the second step we make the approximation that each aperture contains roughly the same signal such that the weights can be set to unity and we furthermore rewrote the expression in a continuous version, which makes the interpretation of the cross term more concise. In particular, we parametrize the lower bound of the integral in terms of the aperture oversampling rate $\alpha \equiv R_{\mathrm{ap}}/\Delta_{\mathrm{min}}$.

In a realistic scenario we expect $\rho$ to to rapidly decrease from unity and then to approach zero for $\Delta \gg R_{\mathrm{ap}}$. An example of such a correlation coefficient is shown in \Fig{fig:SpatialCorrCoeffExample}. Here we explicitly see the importance of taking into account the cross term once there is a substantial overlap between neighbouring apertures. In this example we would infer that measuring the statistics with $\alpha \approx 2$ would be sufficient to extract most of the signal.

\subsection{Implementation and scaling of the direct estimator}

A practical implementation of \Eqn{eq:MapStatEns} consists of three
steps:
\begin{enumerate}
\item Spatially organise the shape catalog to allow for a fast
  assignment of galaxies to apertures.
  %
\item For each aperture of the ensemble compute
  $\hat{\mathfrak{M}}_{\rm ap}$, the associated weight and optionally
  additional systematics (i.e. the coverage fraction $c_k$). Store
  each of these values in an array.
  %
\item Based on some aperture selection and aperture weighting criteria
  $w_{\rm ap}$, update the weights and evaluate the weighted sum.
\end{enumerate}
In what follows we will explore each of these steps in more detail and
for clarity we will denote the number of galaxies in the survey and in
the aperture as $N_{\rm g}$ and $N_{\rm g, ap}$, respectively.

\subsubsection{Assigning galaxies to apertures}

For our implementation we use a spatial hashing data structure. We
start by covering the survey footprint with an equal area mesh of
$N_{\rm pix}$ pixels and create a hash table with the ID of each pixel
as the key and the galaxy IDs as values. The hash function in our case
is the ordinary pixel assignment function. For each aperture we
iterate over the associated galaxies within pixels that partially lie
within the $Q$ filter's support. The construction of the hash table
scales as $\mathcal{O}\left(\max (N_{\rm pix}, N_{\rm g}) \right)$ and the assignment is achieved in $\mathcal{O}\left(\max(N_{\rm pix,ap},
N_{\rm g,ap})\right)$ time per aperture.  We found that when
making a sensible choice of the mesh's coarseness, this data structure
is more stable than a naive KD-tree based implementation as it does not require an additional
range search operation which scales as $\mathcal{O}(\log(N_{\rm g}))$
per aperture and thus becomes a bottleneck for small apertures.

\begin{figure*}
\centering{
    \includegraphics[width=1.\columnwidth]{FIGS/timing_mapgrid_singleR_whash.pdf}
    \includegraphics[width=1.\columnwidth]{FIGS/timing_mapgrid_polyR_whash.pdf}}
  \caption{\small{Computational complexity of the direct estimators
      for equal (left) and unequal (right) aperture radii as a
      function of the number of galaxies. All results are given for
      apertures of radius $10'$ which are oversampled by a factor of
      sixteen $(\alpha=4)$ on a survey field of size $(12
      \rm{deg})^2$. Different colors indicate different orders of the
      evaluated statistics.  The black dashed line indicates the time
      spent in constructing the spatial hash. We see that for equal
      aperture radii the evaluation of higher order statistics
      basically comes for free, while for unequal radii there is a
      constant multiplicative offset based on the relative size of the
      radii and on the order which traces the number of multivariate
      power sums that need to be evaluated. All the scaling were
      obtained when running the estimator on a single CPU
      core. }}\label{fig:DirectEstimatorScaling}
\end{figure*}

\subsubsection{Computing the statistics per aperture}

For the case of all radii being equal we first compute the power sums
in \Eqns{eq:MapnEstPowerSumsA}{eq:MapnEstPowerSumsB} and then
recursively transform them to the corresponding moments via the
recurrence relation \citep{Comtet1974}
\begin{align}\label{eq:BellRecurrence}
B_{n+1}(x_1, \cdots, x_{n+1}) &=
\sum_{i=0}^n \binom{n}{i} B_n(x_1, \cdots, x_n) x_{i+1} \ ,
\end{align}
where $B_0 \equiv 1$. Evaluating each power sum is linear in $N_{\rm
  g, ap}$ and for all practical applications the time taken for
transforming to the $M_{\rm ap}^n$ basis can be neglected.

For the general case we need to compute the relevant multivariable
power sums \Eqn{eq:Msvar} and bring them to the aperture moments basis
by the transformation \Eqn{eq:MapnREst}. In order to
dynamically allocate and evaluate those expressions we use a
combinadic counting scheme to organize the power sum basis whereas the
transformation equation is generated with the help of restricted growth strings \citep{Knuth2005_D}.

\subsubsection{Choice of weights for the averaging}

Following our findings in \citet{Porthetal2020} we employ an inverse
shot noise weighting scheme with an additional hard cutoff $c_{\rm
  lim}$ for the aperture coverage $c_{\rm ap}$, which for second order
statistics was found to lower the mask induced bias while increasing
the signal-to-noise compared to equal weights. The explicit form of
the weights for the $n$th moment can be found from \Eqn{eq:MapStatCovShotLimit}
when neglecting all constant contributions:
\begin{align}\label{eq:InvShotWeight}
\mathcal{W}^{\rm (shot)}_{\rm ap}(c_{\rm lim}) \equiv
\left[\frac{\sum_{(i_1, \cdots, i_n)} w_{i_1}^2 \cdots w_{i_n}^2 }
  {\left(\sum_{(i_1, \cdots, i_n)} w_{i_1} \cdots w_{i_n}\right)^2}\right]^{-1}
\mathcal{H}(c_{\rm ap} - c_{\rm lim}) \ .
\end{align}
Dependent on whether we are dealing with the case of equal or unequal
aperture radii the sums can be decomposed in a similar fashion as
described above and evaluated together with the corresponding
linearised direct estimator. As a further refinement one could also
include the weights and completenesses of the surrounding apertures
weighted by the spatial cross correlation coefficient $\hat{\rho}$ -
this would upweight apertures that are close to a mask as they cover
more unique area.

\section{Results: application to Gaussian mocks}
\label{sec:Measurements}

\subsection{Aperture mass statistics and Gaussian lensing fields}

\begin{figure*}
    \centering
    \centering{
      \includegraphics[width=1.0\columnwidth]{FIGS/gaussmocks_ext_3_Map2VsTheorysigma_0p0.pdf}
      \includegraphics[width=1.0\columnwidth]{FIGS/gaussmocks_ext_3_Map4VsTheorysigma_0p0.pdf}}
    \caption{\small{Comparison of the measured aperture mass moments
        with their theoretical prediction. {\bf Left panel: } The
        upper subpanel shows the aperture mass dispersion as a
        function of the aperture scale. The red line shows the
        theoretical predictions evaluated from the input power
        spectrum and the blue line shows the measurement from the
        mocks. The blue shaded regions show the standard deviation of
        the corresponding measurement across the ensemble. The lower
        subpanel shows the relative error between the measurement and
        the theory, with the line styles as before. {\bf Right panel:
        } Same as left panel, but for the kurtosis of aperture
        mass. }} \label{fig:GaussMock1}
\end{figure*}

In order to validate that our hierarchy of aperture mass moment
estimators are unbiased and do indeed recover correct results, we
first apply them to a set of Gaussian mock lensing simulations. In
this case, the whole moment hierarchy can be written as powers of the
variance of the aperture mass. Hence, this motivates us to define the
scaled aperture mass moments:
\begin{align} \label{eq:GaussCum}
s_n (\vartheta_1)
\equiv
\frac{1}{(n-1)!!}\frac{\MapStatEns{n}(\vartheta_1)}{\left[\MapStatEns{2}(\vartheta_1)\right]^{n/2}}
=
\delta^K_{n,2\mathbb{N}} \ ,
\end{align}
where the final equality is true for a Gaussian field only.

In order to test this we generated 256 Gaussian lensing mocks. The
methodology to create each mock was as follows:
\begin{itemize}
\item We first generate a Gaussian convergence field over a $12
  \times12 \ {\rm square\ deg}$ survey area. The area is tiled by a
  rectangular mesh of $8192^2$ pixels. The variance of the convergence
  is obtained through specifying the convergence power spectrum, and
  we do this for a source distribution similar to that for the \mbox{CFHTLenS}
  \ survey \citep{Fuetal2014}.
\item We next obtain the shear field. This is done by Fourier
  transforming the convergence field and making use of the
  \citet{Kaiseretal1993} approach\footnote{In order to suppress edge
  effects introduced by the FFT we build the pixelated convergence
  field on an a plane having $16$ times the area of the mock.}.
\item We then sample $4\times 10^6$ galaxies into the survey footprint
  and use a multilinear interpolation of the shear field onto each
  galaxy.
\end{itemize}
Note that since we are assessing the accuracy of the estimators only
we choose to set the intrinsic ellipticities of our source galaxies to
zero. On repeating the analysis below when including this term we did
not find a shift of the curves.

\subsection{Computational scaling tests}

Owing to the fact that each step of our algorithm is strictly linear,
we expect a linear relationship between the elapsed time for estimator
evaluation and the number of galaxies, for any given statistic. In
addition, for the equal radius case the order of the statistics should
not strongly impact the evaluation time. However, for the unequal
radius case, this does not necessarily hold true, since the computation
depends on the relative sizes of the apertures as well as on the order
of the statistics to be evaluated.

\Figure{fig:DirectEstimatorScaling} shows the elapsed time of the
direct estimator calculation for a Gaussian mock, where the number of
sampled galaxies in the mock is increased. Focusing on the left panel
first, this shows the case for the standard aperture mass estimators
with equal radii and here we compute all of the moments up to the 20th
order. As expected the computational time for all of the moments
scales linearly with the size of the problem and we also see that
there is no obvious drop in performance for the higher order moments.

The right panel of \Figure{fig:DirectEstimatorScaling} is the same as
the left panel, but now for the case of unequal radii aperture mass
moments, and here we only consider moments up to 6th order. There are
two differences between the equal and non-equal radius case. First, we
can see that there is a much larger multiplicative offset between
adjacent orders for the generalized statistics. This is expected as
the number of basis elements that need to be allocated in that case is
given by $2^n$ compared to the $n$ ones in the equal radius
statistics. We also observe that for the second order statistic the
unequal radius calculation does roughly need four times as long as the
equal radius one. We can explain this offset when noting that for our
example the ratio of the largest and smallest scale was set to
two. With our definition of the oversampling rate as being relative to
the smallest aperture radius this implies that we need to allocate
four times as many galaxies.

Finally, we note that the superior scaling of the direct estimator
compared to traditional estimation methods should not come as a
surprise. Looking back at the original definition \Eqn{eq:MapnKappa}
of the aperture mass one sees that it depends on the positions and
shapes of the galaxies with respect to the aperture origin. In
contrast, when switching to the description of aperture mass
in terms of the shear correlation functions
\citep[i.e.][]{Schneideretal2002a}, the main dependence shifts to the
relative distance and shapes between tuples of galaxies. This change
of reference position makes the evaluation of correlation function
based estimators intrinsically much more complex than a simple
discretization of \Eqn{eq:MapnKappa}.

\begin{figure}
  \centering{
    \includegraphics[width=1.0\columnwidth]{FIGS/gaussmocks_ext_3_GaussConsistencysigma_0p0.pdf}}
  \caption{\small{Scaled $n$th order aperture mass moments $s_n
      (\vartheta)$ (see Eq.\,\ref{eq:GaussCum}), measured in the
      ensemble of 256 Gaussian mocks, as a function of the aperture
      scale, for all moments up to 10th order. The solid lines of
      varying colours show the mean of the measurements. The dotted
      black lines show the Gaussian theoretical expectations. For a
      Gaussian mock, the even order $s_n (\vartheta)$ give unity, and
      the odd ones vanish. }} \label{fig:GaussMock2}
\end{figure}

\begin{figure*}
  \centering{
    \includegraphics[width=1.0\columnwidth]{FIGS/gaussmocks_8192_ext_3_Map4VaryRVsTheorysigma_0p0.pdf}
    \includegraphics[width=1.0\columnwidth]{FIGS/gaussmocks_8192_ext_3_Map6VaryRVsTheorysigma_0p0.pdf}}
  \caption{\small{Multiscale aperture mass moments as a function of
      the scale parameter $\vartheta$, measured in the ensemble of 256
      Gaussian mocks. Line styles are the same as in
      \Fig{fig:GaussMock1}.  {\bf Left panel: } the fourth order
      aperture statistics. In this case, the vector of aperture scales
      was set to $\bm\alpha=(0.5,0.8,1.,2.)$. {\bf Right hand side: }
      Same as left hand side, but this time for the sixth order
      statistic with the vector of aperture scales set to $\bm{\alpha}
      = (0.5,0.7,1.,1.1,1.5,2.)$. }}
    \label{fig:GaussMock3}
\end{figure*}

\subsection{A hierarchy of aperture mass moments}

\Figure{fig:GaussMock1} shows a comparison of the direct estimators
for the second and fourth order aperture mass moments as a function of
angular scale as applied to the 256 Gaussian mocks.  Here we consider
the case where all the aperture radii are equal (recall that for a
Gaussian field all of the odd moments vanish). In both cases the
curves are in very good agreement with the Gaussian theory
predictions, indicated by the solid red lines. We also note that for
increasingly large aperture radii the measured results appear to be
slightly below the theoretical expectation. This discrepancy can be
attributed to finite field effects, as well as border effects being
introduced by the Kaiser-Squires inversion method \citep[see][for a
  discussion]{Piresetal2020}.

\Figure{fig:GaussMock2} presents the measured $s_n$ (see
\Eqn{eq:GaussCum}) for all of the aperture mass moments up to 10th
order as a function of the aperture scale. We see that they are
consistent with the Gaussian theoretical expectations. Note that in
order to obtain this good agreement and circumvent the finite field
effects described above, we used the ensemble mean of the measured
aperture mass variance as the denominator in $s_n$.

\Figure{fig:GaussMock3} displays the fourth and sixth order
multiscale aperture mass statistics as a function of the scale
parameter. Note that there are a number of options for exploring the
configuration dependence of the multiscale aperture mass moments,
here we focus on fixing the ratio of the filter lengths and varying
the overall scale of the configuration $\vec{a}$ with the parameter $\vartheta$,
e.g. for the kurtosis we would have
\begin{equation}
\MapStatEns{4}_{\rm c}(\vec{\vartheta})
\equiv
\MapStatEns{4}_{\rm c}(a_1\vartheta, \cdots,
a_4\vartheta) \
\end{equation}
where the constant $a_i \in \mathbb{R}_+$ specify the configuration. The estimates shown in the figure were obtained
using our generalized estimator \Eqn{eq:MapnREst}. As for the previous
cases, we find good agreement between the measurements and the
Gaussian predictions, which were obtained by making use of
\Eqn{eq:MapnKappaMS} and Wicks theorem for the convergence polyspectra
\citep{Bernardeauetal2002}.

\begin{figure*}
\centering{
    \includegraphics[width=1.\columnwidth]{FIGS/SLICS_KiDS450_SN_10_EqualRadii_2.pdf}
    \includegraphics[width=1.\columnwidth]{FIGS/SLICS_KiDS450_SN_10_EqualRadii_3.pdf}\\
    \includegraphics[width=1.\columnwidth]{FIGS/SLICS_KiDS450_SN_10_EqualRadii_4.pdf}
    \includegraphics[width=1.\columnwidth]{FIGS/SLICS_KiDS450_SN_10_EqualRadii_5.pdf}}
  \caption{\small{Measurement of the aperture mass statistics in the
      SLICS simulation suite for different aperture oversampling rates
      $\alpha$. All measurements were done on an ensemble of 819
      realizations with an angular area of 100 deg$^2$ each, where
      the $n(z)$ follows the KiDS-450 distribution.  The upper part of
      the panels correspond to the mean and rescaled standard deviation from
      the ensembles. The lower panel shows the signal-to-noise for the
      corresponding statistics when rescaled to match a 1000 deg$^2$ survey. For the aperture mass dispersion we
      additionally plot the theoretical prediction as the black
      line. For the fourth and fifth order plots we restrict ourselves
      to the contribution of the connected part of the convergence
      polyspectra. We see that choosing an oversampling parameter of
      $\alpha \gtrsim 4$ recovers most of the
      information.}}\label{fig:SLICSMockEqualSN}
\end{figure*}

\section{Results: Detection significance of higher order moments}\label{sec:DetectionSignificance}

In this section we now turn to the question of the detection
significance of higher order aperture statistics from current and
future surveys.

\begin{figure*}
\centering{
    \includegraphics[width=.8\columnwidth]{FIGS/SLICS_CorrelationMatrixMoments.pdf}
    \includegraphics[width=1.2\columnwidth]{FIGS/SLICS_DetectionSignificanceMoments.pdf}
    }
  \caption{\small{Correlation coefficient matrix (left) and cumulative detection significance (right) for connected moments of the aperture mass statistics. We take into account aperture sizes between $10'$ and $100'$. In the correlation matrix the lower triangle shows the results without shape noise while the upper part includes this term and serves as the basis for the computation of the detection significance.}}
  \label{fig:SLICSMockDetectionSignificanceEqual}
\end{figure*}
\begin{figure*}
\centering{
    \includegraphics[width=1.\columnwidth]{FIGS/SLICS_KiDS450_SN_UnEqualRadii_3.pdf}
    \includegraphics[width=1.\columnwidth]{FIGS/SLICS_KiDS450_SN_UnEqualRadii_5.pdf}
    }
  \caption{\small{Measurements of the unequal radii aperture mass
      statistics of third (left) and fifth (right) order in the SLICS
      simulation suite. Each line corresponds to a different set of
      relative aperture sizes as given in Table
      \ref{tab:MapVaryConfig}.}}
  \label{fig:SLICSMockUnqualSN}
\end{figure*}
\begin{figure*}
\centering{
    \begin{minipage}[b]{.61\columnwidth}
    \includegraphics[width=\linewidth]{FIGS/Hyper2.pdf}\\
    \includegraphics[width=\linewidth]{FIGS/Hyper3.pdf}\\
    \includegraphics[width=\linewidth]{FIGS/Hyper4.pdf}
    \end{minipage}
    \includegraphics[width=1.39\columnwidth]{FIGS/HyperCorrcoef4.pdf}
    }
  \caption{\small{Measurements (left) and correlation matrix (right) of the multiscale aperture mass statistics of second, third and fourth order. For each of those statistics we compute all configurations for ten logarithmically spaced radii between $5'$ and $50'$ in which all the apertures have unequal radii. The black lines indicate the blocks of the (cross-) correlations of different orders.}}
  \label{fig:SLICSMockMultiscaleCorrcoef}
\end{figure*}

\begin{figure}
\centering{
    \includegraphics[width=.8\columnwidth]{FIGS/Indices3D.pdf}
    }
  \caption{\small{Path in which the set of non redundant aperture scales for the third order statistics is traversed. The starting point is the lower corner. The subpath in the $\vartheta_3=50'$ plane corresponds to the full path taken for the second order statistics.}}
  \label{fig:Indices3D}
\end{figure}

\subsection{The SLICS mocks}

In order to answer this question we make use of the
SLICS\footnote{\url{https://slics.roe.ac.uk/}} mocks -- this is a
large suite of lensing mock catalogues generated from a large set of
cosmological $N$-body simulations \citep[for full details
  see][]{Harnois-Derapsetal2018}. Each SLICS mock corresponds to a
survey area of $100 \deg^2$. These are generated from the past light
cone extracted from fully independent gravity-only $N$-body
simulations, which evolve $N=1532^3$ particles within a comoving box
of length $L=505\,\, h^{-1}{\rm Mpc}$. The lensing maps are constructed using the
Born approximation. We adopt the catalogues for which the galaxies are
randomly distributed within the lightcone according to the KiDS-450
source distribution \citep{Hildebrandtetal2017}. The shape noise has
been set to $\sigma_{\epsilon} = 0.29$ per shear component. In order
to mimic a constraining power that is comparable to the KiDS-1000 data while not being too noisy, we rescale the errorbars by a factor of $\sqrt{10}$. This provides us
with effectively $819$ simulated $1000({\rm deg})^2$ surveys with which to perform
our analysis.

When estimating the aperture mass statistics from the SLICS
mocks using the estimator given by \Eqn{eq:MapStatEns}, the
achievable signal-to-noise ratio will depend on the number of sampled
apertures selected. If too few are chosen then our estimate will be
inefficient, on the other hand due to the fact that there are
aperture-to-aperture correlations choosing too many will capture all
of the available information, but ultimately will be computationally
inefficient. We therefore expect that the information will saturate for
a given oversampling rate, and that to sample at a higher rate would
be of little use. To investigate this we proceed as in
\citet{Porthetal2020} and place apertures on a regular grid with
spacing $\Delta$, corresponding to an aperture oversampling rate of
$\alpha \equiv \min\left(\{\vartheta_i\}\right) /2\Delta$.

\subsection{Measurement in the SLICS mocks}

\Figure{fig:SLICSMockEqualSN} shows the detection significance of the
equal radii aperture mass statistics for the second, third, fourth and
fifth order aperture mass statistics as a function of the aperture
scale and for various choices of the oversampling rate. For the second
order statistics we also plot the theoretical prediction of the
aperture mass dispersion evaluated from \Eqn{eq:MapnKappa}, where
the convergence power spectrum was computed with
CCL\footnote{\url{https://github.com/LSSTDESC/CCL}}
\citep{Chisarietal2019} using {\tt Halofit} \citep{Smithetal2003}, but
with the modifications of \citet{Takahashietal2012}, as the matter
power spectrum. While part of the difference between the curves for
small aperture radii could be attributed to uncertainties in {\tt
  Halofit}, our suspicion is that they mostly stem from the limited
particle mass resolution in the SLICS mocks (see Fig. 6 in
\cite{Harnois-Derapsetal2018} for the resulting suppression of the
shear correlation functions for small separations).

Several important points are worth noting from these measurements.
First, we see that for a KiDS-1000 like survey there is sufficient
fidelity to detect the aperture mass statistics up to fourth order\footnote{We find the cumulative detection significance of the fifth order statistics to be at the $2.9\sigma$ level.},
with the signal-to-noise peaking at an aperture size of around
$\vartheta\approx10'$ for all statistics. This is exciting, as this
has never before been achieved with standard correlation function
based estimators, and if correct would represent the first robust
detection of these statistics using these methods. Second, while for
the case of the two-point statistics the signal-to-noise ratio (shown
in lower sub-panels for each plot) falls off slowly for larger
apertures, this ratio approaches zero faster for the connected parts
of the higher order statistics\footnote{Owing to the fact that the
aperture mass has zero mean, the full and connected moments differ
only for even order moments of four or more.}. Third, while an
aperture oversampling rate of $\alpha\approx2$ seems sufficient to
capture all the signal for second order statistics, it misses some
information for subsequent orders where it becomes necessary to use
$\alpha \gtrsim 4$.

\Figure{fig:SLICSMockDetectionSignificanceEqual} displays the correlation structure of the aperture mass cumulants as well as the cumulative detection significance. We only consider measurements with $\vartheta \geq 10'$ as this is where the SLICS mocks do agree reasonably well with the theoretical predictions and and due to the fact that the robust theoretical modelling of those statistics might reach its limits at around those scales. We see that while for shape noise free ellipticity catalogs there are strong correlations for small aperture radii, this is not the case for the realistic mocks in which those scales are still shape noise dominated. We further note large correlations around the diagonal between different orders, where the degree of correlation increases with the order of the cumulants. For the cumulative detection significance we see that that the cumulants beyond third order do not add a substantial amplitude to the cumulative signal-to-noise. This is expected, given the relatively lower signal-to-noise as well as the larger portion of cross-covariances that need to be taken into account. One should note that this type of analysis does not imply that the higher order cumulants are obsolete as they still may add complementary information by breaking cosmological parameter degeneracies.

\subsection{Multiscale aperture mass measurements}

We now shift to the measurement of the multiscale statistics for which there are a number of ways on how to select various aperture scale multiplets. In \Figure{fig:SLICSMockUnqualSN} we focus on a fixed set of aperture propositions and then simply scale them with a single parameter $\vartheta$. The different configurations $\vec{a}$ of aperture radii that we have employed are shown in Table~\ref{tab:MapVaryConfig}. We see that for both, the third and the fifth order moments there does not appear to be a strong decline in detection significance for multiscale apertures compared to the associated moments, even if the relative spread of radii is large.

Another way to select aperture scale multiplets for a statistic of order $n$ is choose a list of $m \geq n$ aperture scales and to compute the statistics for each choice of $n$ elements within that list. For our purposes we choose the subset in which none of the aperture radii are equal, as this speeds up our calculation, see \App{app:ApplicationToMap} for the details. In the left hand side of \Fig{fig:SLICSMockMultiscaleCorrcoef} we show our measurements for the second, third and fourth order connected cumulants of the multiscale aperture mass statistic using ten logarithmically spaced scales between $5'$ and $50'$. The first index of the multiplet corresponds to the selection of the smallest possible aperture scales from which we then start choosing the next lowest radius in the subsequent dimension up until we reach the combination of the largest possible set of aperture radii - for an example of this path for the third order statistics see \Fig{fig:Indices3D}. Recalling that the second order aperture mass statistic is simply a filtered version of the power spectrum, we should not expect the multiscale extension add any information to that order\footnote{For a different form of the $Q$ filter function like the one proposed in \cite{Crittendenetal2001} one can easily work this out analytically, see .i.e. \cite{Schneideretal2005a}.}. For the three statistics we again find a detection significance that is comparable to the equal scale case, meaning that we can extract substantial signal from convergence spectra configurations which are not corresponding to regular polygons. In the right hand side of  \Fig{fig:SLICSMockMultiscaleCorrcoef} we plot the joint correlation coefficient of the multiscale cumulants. On the investigated range of scales we only find a slight to modest correlation between the higher order multiscale statistics and the second order one. It also appears that the higher order cumulants exhibit a stronger auto- and cross correlation. However, this is (at least partially) an artefact of the range and sampling density of the chosen radii.

\begin{table}
    \centering
    \begin{tabular}[h]{|c c|c c|}
    \hline
    \multicolumn{2}{|c|}{Third order} &  \multicolumn{2}{c|}{Fifth order} \\
    Label  & Configuration $\vec{a}$ & Label & Configuration $\vec{a}$ \\ \hline
    $X_1$ & $(1, 1, 1)$ & $X_5$ & $(1, 1, 1, 1, 1)$ \\
    $X_2$ & $(1, 2, 2)$ & $X_6$  & $(1, 1, 2, 2, 2)$ \\
    $X_3$ & $(1, 5, 5)$ & $X_7$  & $(1, 1, 5, 5, 5)$ \\
    $X_4$ & $(1, 3, 5)$ & $X_8$ & $(1, 2, 3, 4, 5)$ \\
    \hline
    \end{tabular}
    \caption{Cofigurations of the aperture radii displayed in \Fig{fig:SLICSMockUnqualSN}}
    \label{tab:MapVaryConfig}
\end{table}

\section{Conclusions and Discussion}\label{sec:Conclusions}

In this paper we have explored an alternative method for estimating
the aperture mass statistics in weak lensing cosmic shear
surveys. This study extended our previous work \citep{Porthetal2020}
in a number of ways: First, we generalized the direct estimator
approach to higher statistics, and showed how to rewrite the standard
estimator as a product of linear order time sums.  Second, we provided
the details of the computation of the variance of these estimators.
Third, we further generalised the aperture mass statistics to include
the multiscale approach. Again, we showed how one could estimate
these using linear order products of power sums. The work can be
summarised as follows:

In \Sect{sec:Theory} we reviewed the background theory of cosmological
weak lensing and showed how the connected cumulants of the aperture
mass statistics are related to the convergence polyspectra.

In \Sect{sec:Estimators} we introduced the direct estimator for
moments of the aperture mass statistics. We then gave expressions for
how the nested sums can be decomposed into a linear combination of
products of (multivariate) power sums that facilitates a linearly
scaling estimation procedure in the number of galaxies within an
aperture. We then generalized this estimator to an ensemble of
overlapping apertures and computed its variance. We argued that the
aperture cross correlation coefficient leads to a substantial
correction to the naive $1/N$ scaling if the apertures are not well
separated, and that it also can be used to assess the degree of
aperture oversampling that is necessary to capture most of the
available information. Finally, we gave a detailed explanation of the
algorithms used for our implementation.

In \Sect{sec:Measurements} we successfully validated our method on
Gaussian mock simulations and furthermore verified the linear scaling.

In \Sect{sec:DetectionSignificance} we turned to the SLICS simulation
suite and assessed the signal-to-noise of the statistics for a $1000$
degree survey following a KiDS-450 like $n(z)$ distribution
function. We found that with these specifications significant
detections of up to fourth order can be expected for the equal and
unequal radii cumulants and that an aperture oversampling rate of at
least four extracts nearly all the signal.

In this paper we have neglected the impact of survey masks on the
measurement process and the possible bias that this could induce, the
exploration of this is sufficient to warrant its own publication and
this is the subject of our associated publication (Porth et al. in
prep.).  Throughout this paper we were mainly concerned with making
the extraction of information from higher order statistics of galaxy
shape catalogs computationally feasible and accurate. However, we
remained agnostic about further challenges that need to be addressed
before applying our methods to real data. For example, one should
investigate the required PSF modelling, shape measurement and shear
bias calibration quality to not introduce substantial biases in the
measurement. Additionally, the range of measurements that can
ultimately be used for obtaining cosmological parameter constraints
will be limited to the scales for which one can theoretically
accurately model those higher order statistics.

\section*{Acknowledgements}
We would like to thank Patrick Simon, Laura Marian, Stefan Hilbert, Peter Schneider, Cora Uhlemann and Gary Bernstein for useful discussions, as well as the anonymous referee for helpful comments. We would like to thank Joachim
Harnois-Deraps for making public the SLICS mock data, which can be
found at http://slics.roe.ac.uk/.
LP acknowledges support from a STFC Research Training Grant (grant
number ST/R505146/1). RES acknowledges support from the STFC (grant
number ST/P000525/1, ST/T000473/1). This work used the DiRAC@Durham
facility managed by the Institute for Computational Cosmology on
behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment
was funded by BEIS capital funding via STFC capital grants
ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham
University and STFC operations grant ST/R000832/1. DiRAC is part of
the National eInfrastructure. This research made use of numpy, a
library used for scientific computing and technical computing and
matplotlib, a Python library for publication quality graphics
\citep{Harris2020, Hunter2007}.

\section*{Data Availability}
The SLICS mock catalogs are available at \url{https://slics.roe.ac.uk/}. Additional data underlying this article will be shared on reasonable request to the corresponding author.

\bibliographystyle{mnras}


text:
\section{Method}

We use the partial reconstruction of the
process $\ee \to \dss D_{s1}^-$. In this analysis 4-momentum
conservation allows us to infer the 4-momentum of the undetected
part. The method used  was described in Ref.~\cite{DstDst} and applied to the
measurement of the $\ee \to D^{(*)+} D^{(*)-}$ cross sections.

Here we reconstruct the process $\ee \to \dss D_{s1}^-$ using two
different tagging procedures. The first one (denoted as the \dso\ tag)
includes the full reconstruction of the \dso\ meson via $\dso\to
\overline{D}{}^* K$ decay and observation of the photon from $\dss \to
\ds \gamma$, while the $D_s^+$ is not reconstructed. The measured
signal yield with the \dso\ tag is proportional to the branching fractions
of the reconstructed $\overline{D}{}^*$ modes. In the second procedure
(denoted as the \dss\ tag) we require full reconstruction of \dss\
through $\dss \to \ds \gamma$ and observation of the kaon from
$D_{s1}^{-}\to \overline{D}{}^{*} K$, but the $\overline{D}{}^{*}$ is
not reconstructed. Since the \ds\ meson is reconstructed in the
channel of interest, $\ds\to K^+ K^- \pi^+$, the signal yield measured
with the \dss\ tag is proportional to this \ds\ branching
fraction. The (efficiency-corrected) ratio of the two measured signal
yields is equal to the ratio of well-known $\overline{D}{}^{*}$
branching fractions and the branching fraction of the \ds.

In order to calculate reconstruction efficiencies and optimize event
selection criteria,  Monte Carlo (MC) signal events are generated and
simulated using a GEANT-based full simulator, including  initial
state radiation (ISR), and assuming no form-factors for \dss\ and \dso\ mesons.

To identify the signal we study the mass recoiling against the
reconstructed particle (or combination of particles) denoted as $X$. This recoil mass is defined  as:
\begin{eqnarray}
\RM(X) \equiv \sqrt{ {(E_{CM} - E_X)}^2 - P_X^2 },
\end{eqnarray}
where $E_X$ and $P_X$ are the center-of-mass (CM) energy and momentum of X,
respectively; $E_{CM}$ is the CM beam energy. We expect a peak in the \RM\ distribution at the
nominal mass of the recoil particle.

The resolution in \RM\ is   $\sim 50\ensuremath{\,{\mathrm{MeV}/c^2}}$ according to the MC, which is not sufficient to separate different final states, {\emph e.g.}  $\ds
D_{s1}^{-}$, $\dss D_{s1}^{-}$ and non-resonant $D_s^+ \overline{D} K$. To
disentangle the  contribution of these final states we use another
kinematic variable, the recoil mass difference \RMD:
\begin{eqnarray}
\RMD(D^{-}_{s1} \gamma) \equiv \RM(D^{-}_{s1})-\RM(D^{-}_{s1}\gamma), \\
\RMD(D^{*+}_s K) \equiv \RM(D^{*+}_s)-\RM(D^{*+}_s K).
\end{eqnarray}
In the {\dso\ tag} procedure the signal events make a narrow peak in
the $\RMD(D^{-}_{s1} \gamma)$ distribution at the nominal $\dss -\ds$ mass
difference with a resolution of $\sigma \sim 5\ensuremath{\,{\mathrm{MeV}/c^2}}$, dominated by the
$\gamma$ energy resolution according to the MC. The  $\RMD(D^{*+}_s K)$
spectrum peaks at the $D_{s1}^{-} - \overline{D}{}^{*}$ mass
difference with a resolution of $\sigma \sim 2\ensuremath{\,{\mathrm{MeV}/c^2}}$.

\section{Selection}

All primary charged tracks are required to be consistent with
originating from the interaction point. Charged kaon candidates are
identified using information from $dE/dx$ measurements in the central
drift chamber, aerogel Cherenkov counters and time-of-flight
system. No identification requirements are applied for pion
candidates. \ks\ candidates are reconstructed from $\pi^+
\pi^-$ pairs with an invariant mass within $15\ensuremath{\,{\mathrm{MeV}/c^2}}$ of the nominal
\ks\ mass. The distance between the two pion tracks at the \ks\ vertex
is required to be smaller than $1\,\mathrm{cm}$, the flight distance
in the plane perpendicular to the beam from the interaction point is
required to be greater than $0.1\,\mathrm{cm}$ and the angle between
the \ks\ momentum direction and decay path in this plane should be
smaller than $0.01\,\mathrm{rad}$. Photons are reconstructed in the
electromagnetic calorimeter as showers with energies above $50 \mev$
that are not associated with charged tracks. $\pi^0$ candidates are
reconstructed by combining pairs of photons with invariant masses
within $15\ensuremath{\,{\mathrm{MeV}/c^2}}$ of the nominal $\pi^0$ mass.

$D^0$ candidates are reconstructed using five decay modes: $K^-
\pi^+$, $K^- K^+$, $K^- \pi^- \pi^+ \pi^+$, $\ks \pi^+\pi^-$ and $K^-
\pi^+ \pi^0$. $D^+$ candidates are reconstructed using the $K^- \pi^+
\pi^+$ channel. \ds\ candidates are reconstructed using the $K^+ K^-
\pi^+$ decay mode. A $\pm 15\ensuremath{\,{\mathrm{MeV}/c^2}}$ mass window is used for all $D$
modes (approximately $2.5\sigma$ in each case) except for the $D^0
\to K^- \pi^- \pi^+ \pi^+$ where $\pm 10\ensuremath{\,{\mathrm{MeV}/c^2}}$ is applied and $D^0\to
K^- \pi^+ \pi^0$ with $\pm 20\ensuremath{\,{\mathrm{MeV}/c^2}}$ mass window.  All \ds, $D^+$ and
$D^0$ candidates are subjected to a mass and vertex constrained fit to
improve their momenta and thus the recoil mass resolution.

$D^*$ candidates are selected via $D^{*+} \to D^0 \pi^+$ and $D^{*0}
\to D^0 \pi^0$ decay modes with $\pm 2\ensuremath{\,{\mathrm{MeV}/c^2}}$ mass window. \dss\
candidates are reconstructed via $\ds \gamma$ channel with $\pm
10\ensuremath{\,{\mathrm{MeV}/c^2}}$ mass window. \dso\ is reconstructed in $\overline{D}{}^{*0}
K^-$ and $D^{*-} \ks$ decay modes.

In the case of multiple candidates, the candidate with the minimum value of
$\chi^2_{tot}$ ($\chi^2_{tot}=\chi^2_{M(D)} + \chi^2_{M(D^*)}$ for \dso\ tag
and $\chi^2_{tot}=\chi^2_{M(\ds)}$ for \dss\ tag, respectively) is chosen. Each
$\chi^2$ is defined as the square of the ratio of the deviation of the
measured mass from the nominal value and the corresponding resolution.

\section{Reconstruction}

\subsection{\dso\ tag}

As the ratio of $\dso\to\dstbkm$  and $\dso\to
\dstmks$ branching fractions is unknown, we perform the analysis for
these two channels separately and average the results.

The \dstbkm\ and \dstmks\ mass spectra with a preselection requirement
on $\RM(\overline{D}{}^*K)~<3~\ensuremath{\,{\mathrm{GeV}/c^2}}$ are shown in
Fig.~\ref{ds1.mds1}. Clean peaks at the nominal \dso\ mass without
substantial combinatorial background are evident in both
distributions. We fit these distributions with a  sum of two Gaussians with  a common mean
 and   the background parametrized with a  threshold
function:
\begin{equation}
Bg(m)=(A\cdot m + B )\sqrt{m-M^{PDG}_{\overline{D}{}^*}-M^{PDG}_{K}}.
\label{ds1.bg}
\end{equation}
The fit yields a Gaussian central value of $M=2535.1\pm
0.4\ensuremath{\,{\mathrm{MeV}/c^2}}$ for \dstbkm\ and $M=2535.0\pm 0.3\ensuremath{\,{\mathrm{MeV}/c^2}}$ for \dstmks\
channels, respectively, which is in  good agreement with both the PDG
and MC values.

\begin{figure}[htb]
\begin{center}
\epsfig{file=ds1.mass.eps,width=0.7\textwidth}
\end{center}
\caption{The distribution of invariant mass a) \dstbkm\ and b)
\dstmks. The solid curves represent the fit described in the text.}
\label {ds1.mds1}
\end{figure}

The signal and sideband regions for \dso\ are defined as $|M(\dso)-2.5354|<0.005\ensuremath{\,{\mathrm{GeV}/c^2}}$ and $0.010\ensuremath{\,{\mathrm{GeV}/c^2}}<|M(\dso)-2.5354|<0.015\ensuremath{\,{\mathrm{GeV}/c^2}}$, respectively. The $\RM(\dso)$ spectra  are
shown in Fig.~\ref{ds1.rm}.  The significant excess in the region from
$1.9$ to $2.2 \ensuremath{\,{\mathrm{GeV}/c^2}}$ includes $\ee \to \ds \dso$ and $\ee \to \dss
\dso$ processes.   The small
excess around $2.3 \ensuremath{\,{\mathrm{GeV}/c^2}}$ could be interpreted as a
contribution from the process $e^+e^-\to  D_{sJ}^+\dso$.  The excess
around $2.5 \ensuremath{\,{\mathrm{GeV}/c^2}}$ corresponds to the excited $D_s^{**+}=D_{s1}^+, ~
D_{s2}^+$ states.

\begin{figure}[htb]
\begin{center}
\epsfig{file=ds1.recoilmass.eps ,width=0.7\textwidth}
\end{center}
\caption{a) $\RM(\dso \to \dstbkm)$ distribution; b) $\RM(\dso \to
\dstmks )$ distribution.  Points with error bars correspond to the
events in the signal region. The \dso\ sideband contribution is shown as hatched
histogram. The  signal windows are indicated by the
vertical lines. }
\label{ds1.rm}
\end{figure}

We define the signal $\RM(\dso)$ window as
$|\RM(\dso)-2.1121|<0.100\ensuremath{\,{\mathrm{GeV}/c^2}}$ (Fig.~\ref{ds1.rm}). As the recoil mass resolution of $\sim50\ensuremath{\,{\mathrm{MeV}/c^2}}$ is
not sufficient to resolve  $\ee \to \ds \dso$ and $\ee \to \dss
\dso$ processes, to disentangle them we
reconstruct the $\gamma$ from $\dss \to  \ds\gamma$ decays. The reconstructed
\dso\ candidates from the signal $\RM(\dso)$ region are combined with
each photon in the event to calculate the recoil mass difference
$\RMD(\dso\gamma)$.

To obtain the signal yield and disentangle the signal contribution
from the non-resonant $\dss \overline{D}{}^* K$ process, we perform fits to $\overline{D}{}^* K$
mass distributions in bins of $\RMD(\dso\gamma)$. The
$\RMD(\dso\gamma)$ bin width is chosen to be $4\ensuremath{\,{\mathrm{MeV}/c^2}}$. The signal
function is the sum of two  Gaussians with the mass and width fixed to the corresponding MC
values. The background is described with a threshold function
(Eq. \ref{ds1.bg}).
The signal yield for each bin  of the recoil mass difference for the $\dso\to\dstbkm$ and $\dso\to
\dstmks$ channels are plotted in Fig.~\ref{ds1.binfit.results}.

\begin{figure}[h!]
\begin{center}
\epsfig{file=ds1.results.eps,width=0.7\textwidth}
\end{center}
\caption{The signal yields in bins of $\RMD(\dso\gamma)$: a) for
the $\dso \to \dstbkm $ channel; b) for the $\dso \to \dstmks $
channel. The solid curves represent the fits described in the text.}
\label {ds1.binfit.results}
\end{figure}

The method used to extract the signal yield guarantees that  the
signal peak includes  only those events that contain both a \dso\
meson and a \dss\ decaying to $D_s\gamma$. Thus, the only possible
peaking background sources are  final states with extra pion(s)
({\emph e.g.} $\ee \to D^{*+}_s D_{s1}^- \pi^0$). The process with one
extra $\pi^0$ is suppressed due to isospin conservation. The possible
contribution of  $\ee \to \dso D_{sJ}^+(2460)$ followed by
$D_{sJ}^+(2460)\to \dss\pi^0$ decay, where the isospin violation
occurs in the $D_{sJ}$ decay, is suppressed by the recoil mass cut at the $4
\sigma$ level. The processes with two pions are also suppressed at
the $\sim3.5 \sigma$ level. Therefore, we conclude that the remaining
background is negligibly small; the extracted yields are
dominated by the signal.

The distributions in Fig.~\ref{ds1.binfit.results} are fitted to the
sum of a signal Gaussian (with the mass and width fixed to the MC
values) and a linear background function. The fits yield
\yieldone\ signal events for $\dso \to\dstbkm$ and \yieldtwo\ events
for $\dso\to\dstmks$. To determine the reconstruction efficiency the
same procedure is applied to the MC signal sample.

\subsection{\dss\ tag}
In the \dss\ tag,  the \dss\ of  the  process $\ee\to\dss\dso$  is fully reconstructed.
The $\ds\gamma$ mass distribution is shown in
Fig.~\ref{dsst1.mdsst} and fitted with the sum of a signal Gaussian
and a linear background function. The obtained central value of the
Gaussian $M=2113.5\pm 0.4\ensuremath{\,{\mathrm{MeV}/c^2}}$ is somewhat larger than the  PDG value, but
 in  good agreement with the MC expectation, which shows a similar shift.  The signal region is
chosen to be $|M(\ds\gamma)-2.1121|<0.010\ensuremath{\,{\mathrm{GeV}/c^2}}$ while the sidebands
are defined as $0.015\ensuremath{\,{\mathrm{GeV}/c^2}}<|M(\ds\gamma)-2.1121|<0.025\ensuremath{\,{\mathrm{GeV}/c^2}}$.

\begin{figure}[htb]
\begin{center}
\epsfig{file=dsst1.mass.eps,width=0.7\textwidth}
\end{center}
\caption{ The $\ds \gamma$ invariant mass distribution. The solid curve represents the fit described in the text.}
\label {dsst1.mdsst}
\end{figure}

The $\RM(\dss)$ spectrum for the signal region is shown  in
Fig.~\ref{dsst1.rm}. A wide enhancement around $2.1\ensuremath{\,{\mathrm{GeV}/c^2}}$ is formed by
overlapping $D_s^-$ and $D_s^{*-}$ signals. The signal process $\ee\to
\dss \dso$ contributes to the bump around $2.55\ensuremath{\,{\mathrm{GeV}/c^2}}$ together with the
process $\ee \to  \dss\dst$. We define the signal region by the
requirement $|\RM(\dss)-2.5354|<0.100\ensuremath{\,{\mathrm{GeV}/c^2}}$ ($\sim2\sigma$).

\begin{figure}[htb]
\begin{center}
\epsfig{file=dsst1.recoilmass.eps,width=0.7\textwidth}
\end{center}
\caption{The $\RM(\dss)$ distribution. The events in the \dss\ signal region are shown as points
with error bars.  The \dss\ sidebands are superimposed as hatched
histogram. The  $\RM(\dss)$ signal window is indicated by
vertical lines. }
 \label {dsst1.rm}
\end{figure}

To extract the $\ee \to  \dss\dso$ signal yield for \dss\ tag, we perform a fit to
the $\ds \gamma$ mass distributions in bins of  $\RMD(\dss K)$. The
signal function is a  Gaussian with its parameters fixed to the corresponding  MC
values. The background  is parametrized by a  linear function.
The yields returned by the fits for $\dso\to\dstbkm$ and $\dso\to
\dstmks$ channels are plotted in Fig.~\ref{dsst1.binfit.results}. We
fit the resulting distribution with the sum of a signal Gaussian (with
fixed mass and width) and a threshold function (Eq. \ref{ds1.bg}). The
fits yield \yieldthree\ signal events for the \dstbkm\ channel and
\yieldfour\ events for the \dstmks\ channel.
\begin{figure}[h!]
\begin{center}
\epsfig{file=dsst1.results.eps,width=0.7\textwidth}
\end{center}
\caption{The $\ee\to\dss\dso$ yields for \dss\ tag in bins of $\RMD(\dss K)$: a) for
the $\dso \to \dstbkm $ channel; b) for the $\dso \to \dstmks$
channel. The solid curves show the fits described in the text.}
 \label {dsst1.binfit.results}
\end{figure}

\section{Systematic errors}

The systematic uncertainties are summarized in Table \ref{systematics}.
\begin{table}[h!]
\caption{Contributions to the systematic error.}
\begin{tabular}{@{\hspace{0.3cm}}l@{\hspace{0.3cm}}||@{\hspace{0.3cm}}c@{\hspace{0.3cm}}}
\hline \hline
Source & Estimated error,  [\%] \\
\hline
Production angle     &  0.4  \\
$\overline{D}{}^*$ polarization   &  3.8  \\
ISR                  &  2.3  \\
Recoil mass cut      &  4.8  \\
Fit                  &  4.3  \\
Tracking/$\pi^0$ reconstr. efficiency             &  3.0  \\
${\cal B} (\overline{D}{}^{(*)})$ &  5.2  \\
MC statistics        &  2.3  \\
\hline
Total                &  10.1 \\
\hline
\hline
\end{tabular}

\label {systematics}
\end{table}

Different production and helicity angular distributions result in
different ratios of reconstruction efficiencies for the two tagging
procedures. A potential source of significant systematic error are the uncertainties,
associated with reconstruction efficiencies of the kaon from the \dso\ decay
and the photon from the \dss, which strongly depend on the mother particle
polarizations. As the kaon and photon are, however,  reconstructed in both tags, the
systematic uncertainties related to the reconstruction of these
particles cancel out.  The remaining  systematic uncertainty in  the
efficiency ratio from  the production angle $\theta_{\mathrm{prod}}$ (polar angle
of the momentum of \dso\ meson) is estimated from MC simulation by
comparing the regions $|\cos(\theta_{\mathrm{prod}})|>0.5$ and
$|\cos(\theta_{\mathrm{prod}})|<0.5$ and found to be small
($0.4\%$).

The  systematic error arising from unknown $\overline{D}{}^*$
polarizations contributes only to the uncertainty in the efficiency of
\dso\ tag and does not cancel out in the efficiency ratio. This
contribution is estimated from MC assuming  different $\overline{D}{}^*$
polarizations.

 The efficiency of the requirement on  recoil mass depends on the
fraction of events with an ISR photon as well as on the recoil mass
window cut. The error from the ISR photon is conservatively estimated
by changing the fraction of events containing ISR photons in the whole
MC sample by $30\%$. The resulting relative shift in ${\cal B}(\ds)$ is
found to be $2.3\%$. The error from the recoil mass cut is estimated
by varying the recoil mass window width in the range of
$80-120\ensuremath{\,{\mathrm{MeV}/c^2}}$. The difference in the calculated value of ${\cal B}(\ds)$ ($4.8\%$) is
conservatively considered as the systematic error.

The fit systematics is estimated assuming  different shapes for the
signal (single and double Gaussians) and for the background (different
order polynomials). The difference in the resulting branching fraction
is found to be $4.3$\% and is treated as the systematic uncertainty in  the
fit.

Other contributions come from the uncertainties in the track/$\pi^0$
reconstruction efficiencies, $\overline{D}{}^{(*)}$ branching fractions and MC statistics.

All of the systematic contributions   described above are added in quadrature to give $10.1$\%.

\section{Results}

Using the measured signal yields $N(\dss)$ and $N(\dso)$ with \dss\
and \dso\ tags, respectively, and taking into account the
efficiency ratio ($\frac{ \epsilon(\dso)}{\epsilon(\dss)}=0.50$ for
\dstbkm\ and $\frac{ \epsilon(\dso)}{\epsilon(\dss)}=1.65$ for
\dstmks\ channels, respectively), we obtain the \ds\ absolute
branching fraction:
\begin{eqnarray}
  \mathcal{B} (\ds \to K^+ K^- \pi^+ )= \frac{N(\dss)}
  {N(\dso)}\cdot \frac{ \epsilon(\dso)}{\epsilon
  (\dss)}\cdot \mathcal{B}(\overline{D}{}^{(*)}),
\end{eqnarray}
where $\mathcal{B}(\overline{D}{}^{(*)})$ is the product of
$\overline{D}{}^*$ branching fraction and those of sub-decays. Separate calculations  for the \dstbkm\ and \dstmks\ channels
yield $\mathcal{B}(\ds \to K^+ K^- \pi^+)$ of (\resultsone)\% and
(\resultstwo)\%, respectively.

\section{Summary}

In summary, we report a measurement of the $\ds \to K^+ K^- \pi^+ $
branching fraction using a new method of  double tag partial
reconstruction. The branching fraction is measured separately in two channels
$\ee \to \dss\dso(\to \dstbkm)$ and $\ee \to \dss\dso(\to \dstmks)$. The average value is $\mathcal{B}(\ds \to K^+ K^-
\pi^+)=(\results$)\%.

redpajama_set_name: RedPajamaArXiv