\section{...}, \subsection{...}, \subsubsection{...}, \paragraph{...}, \subparagraph{...}

\defcitealias{Herenz2017}{HW17}

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

Modern wide-field integral-field spectrographs \citep[IFS;][]{Bacon2017a} on
large ground based telescopes offer unprecedented capabilities for spectroscopic
surveys of faint emission line galaxies.  Currently these instruments are
represented by the Multi-Unit Spectroscopic Explorer on ESO's Very Large
Telescope ``Yepun'' at Cerro Paranal \citep[MUSE;][]{Bacon2014} and the Keck
Cosmic Web Imager on Keck II at Manua Kea \citep{Morrissey2018}.  The continuous
spatial and spectral coverage provided by wide-field IFS avert the classical
spectroscopic survey paradigm of target pre-selection for follow-up
spectroscopy, as these observations provide literally ``spectroscopy of
everything'' within the observed fields.  This concept is also utilised in the
Hobby Eberly Telescope Dark Energy Experiment (HETDEX), a dedicated IFS survey
facility for Lyman $\alpha$ (Ly$\alpha$\,$\lambda1216$) emitting galaxies
(LAEs), that aims at constraining the Dark Energy equation of state
\citep{Hill2008,Gebhardt2021,Hill2021}.  Given their pivotal importance in
addressing fundamental astrophysical and cosmological questions the next
generation of wide-field IFS are already envisioned \citep[e.g., as part of the
ESO Spectroscopic Facility;][]{Pasquini2018} or in early stages of planning
\citep[BlueMUSE for ESO/VLT;][]{Richard2019}.

Efficient detection and cataloguing of sources in IFS datasets is obviously a
crucial task for fully exploiting the survey capabilities of these instruments.
This analysis step follows in sequence to the data reduction of the
observational raw data.  The data reduction pipelines processes the instrumental
detector data of the on-sky exposures into science-ready data-structures (see,
e.g., \citealt{Weilbacher2020} for the MUSE data reduction pipeline and
\citealt{Turner2010} for general aspects of IFS data reduction).  The reduced
high-level data-products are three dimensional (3D) arrays, with two axes
mapping the spatial domain (on-sky position / Right Ascension and Declination)
and the third axis mapping the spectral domain (wavelength).  One element of
this 3D array is called a voxel (short form of volume pixel), and the collection
of all voxels at the same spatial coordinate are called spaxel (spectral pixel).
The jargon for such IFS data-products is ``datacubes'' in literature
\citep[e.g., review by][]{Allington-Smith2006}, but since generally the length
of the spatial and spectral axes differ we will here use the correct term
``data cuboid''.

In the last decade several works were concerned with the problem of detecting
and cataloguing faint emission line sources in data cuboids from deep MUSE
observations \citep[e.g.,][]{Bourguignon2012}.  This happened because already
first science case that motivated the construction of MUSE envisioned it as a
discovery machine for the faintest line emitting objects. \citep{Bacon2002}; a
dream that has now become reality
\citep[e.g.,][]{Maseda2018,Bacon2021,SanchezAlmeida2022}.  As of this writing
the following software implementations for detecting and cataloguing emission
lines in wide-field IFS data cuboids are freely\footnote{Another software is
  \texttt{CubeEx}, but this programme is only available from the author on
  request \citep[see the the ``code availability'' remark in][]{Ginolfi2022}.}
available: ``Source Emission Line FInder'' \citep[SELFI;][]{Meillier2016},
ORIGIN \citep{Mary2020}, MUSELET \citep[a wrapper of the popular imaging source
detection software SExtractor and part of the MUSE Python Data Analysis
Framework;][]{Piqueras2017}, and the ``Line Source Detection and Cataloguing
Tool'' \mbox{LSDCat} \citep[][hereafter
HW17]{Herenz2017}.

LSDCat distinguishes itself from ORIGIN or SELFI by a rather simplistic
detection and cataloguing approach.  It is based around a straightforward
two-step process of 3D linear (matched) filtering and thresholding; an optional
third step can be invoked to parameterise the found emission line sources
\citepalias{Herenz2017}.  The simplicity allows for fast processing of MUSE
data cuboids, which comprise typically data volumes of 4\,GByte per pointing
(flux data and corresponding variances).  Fast processing is a requirement when
working with a large number of such data cuboids, but it allows also for
efficient source insertion and recovery experiments to empirically determine the
selection function for some particular type of emission line galaxies.  Robust
selection function determinations are a must when estimating and modelling
fundamental galaxy population statistics.  Worked out examples in this respect
are the spatial clustering measurement \citep{HerreroAlonso2021} as well as the
luminosity function determination \citep{Herenz2019} of Ly$\alpha$ emitting
galaxies at redshifts $3 \lesssim z \lesssim 6$ in the MUSE-Wide survey.

In the present article we present an improvement of the filtering algorithm of
\mbox{LSDCat}\footnote{LSDCat can be obtained from the Astrophysics Source Code
  Library: \url{https://ascl.net/1612.002}.}.  The improved method accounts for
high-frequency variations of the noise in spectral direction within IFS
data cuboids.  A high incidence of such high-frequency variations at
considerable strength occurs at $\lambda \gtrsim 7000$\,\AA{} in ground based
observations due to the atmospheric spectral lines of OH-molecules \citep[see
Sect.~4.3 in][]{Noll2012}.  Our revision of the original algorithm is motivated
primarily by the desire to improve the detectability of the faintest Ly$\alpha$
emitting galaxies at the highest redshifts in the deepest MUSE-datasets,
i.e.,  the MUSE Hubble Ultra Deep Field
\citep[$t_\mathrm{exp} \gtrsim 30$\,h;][]{Bacon2017,Inami2017} and the MUSE
eXtreme Deep Field \citep[MXDF, $t_\mathrm{exp} =141$\,h;][]{Bacon2022}.  Here
it was found that robust detections with ORIGIN appear not as significantly
detected with LSDCat while others could even be missed at reasonable detection
thresholds (R.~Bacon, priv. comm.).  Another motivation is to provide a
deterministic detection algorithm, whose selection function is solely defined by
the noise properties of the data and the physical properties of the coveted
emission line sources (i.e., , their line flux, as well as their spatial-, and
spectral profiles).

While this article accompanies the release of the improved version of the
\mbox{LSDCat} software (\mbox{LSDCat2.0}), its aim is also to provide a robust
general formal framework for the detection of emission line sources in
integral-field spectroscopic datasets via matched filtering.  The remainder of
this article is structured as follows: In Sect.~\ref{sec:2} we review some basic
properties of source detection by matched filtering that are of relevance for
the task at hand (Sect.~\ref{sec:mf}).  We then briefly review the original
algorithm in LSDCat (Sect.~\ref{sec:classic}) before explaining the
modifications in the improved version (Sect.~\ref{sec:modified}).  Next we
contrast the improved algorithm to the previous version and we verify the
implementation against the detected emission lines within the first public data
release of the MUSE-Wide survey (Sect.~\ref{sec:verify}).
Section~\ref{sec:selfun} then describes the construction of the deterministic
selection function of the detection metod.  First, we derive the selection
function for the idealised case where the spectral profiles of the surveyed
emission line sources are exactly known and where the spatial profiles are
indistinguishable from the instrumental point spread function
(Sect.~\ref{sec:ideal-select-funct}).  This idealised selection function is then
verified against an empirical source insertion and recovery experiment in the
recently released MUSE data cuboid of the MUSE eXtreme Deep Field
(Sect.~\ref{sec:verif-analyt-expr}).  Lastly, we present a re-scaling method of
the idealised selection function to different source profiles
(Sect.~\ref{sec:rs}).  We close the paper with a summary and conclusion in
Sect.~\ref{sec:sc}.


\section{An Improved 3D Line Detection Algorithm}
\label{sec:2}

Source detection is, at its heart, a statistical test to rule out the hypothesis
$\mathcal{H}_0$ of no source being present somewhere in the observational
dataset under examination.  If $\mathcal{H}_0$ can be rejected with confidence,
we accept the reasonable alternative, $\mathcal{H}_1$, that the signature from
the sought astronomical source is present in our dataset.  \mbox{LSDCat}, like
virtually every astronomical source detection software, is based around linear
filtering of the data prior to deciding between $\mathcal{H}_1$ and
$\mathcal{H}_0$.

Arguably, the most suitable filter for the detection of isolated astronomical
emission line sources in a background noise dominated observational dataset is
the ``matched filter''; the name resonates that the response of the filter is
matched to the known or expected shape of the signal to be detected.  Matched
filters provide the optimal test statistic to decide between $\mathcal{H}_0$ and
$\mathcal{H}_1$, i.e.,  no other statistic provides a better decision criterion
\citep[see, e.g.,][for a textbook treament of matched filtering]{Das1991}.  A
recent review on advanced matched filtering concepts was presented by
\cite{Vio2021} and we follow in some parts the notation developed in there.  The
applications of the matched filter are widespread in science and engineering, e.g.,
detection of gravitational waves \citep{Jaranowski2012} or radar
\citep{Schwartz1975}, and in the following we describe how this concept is being
applied to the emission line source detection problem in IFS data.

\subsection{Matched Filtering for Emission Lines in IFS data cuboids}
\label{sec:mf}

Considering first the one-dimensional discrete case with Gaussian noise allows
us to formalise the line detection problem as follows.  We consider a line
profile that represented by a shape vector
\begin{equation}
  \label{eq:1}
  \bm{s} = (s_1, \dots, s_M)^T  \;\text{,}
\end{equation}
which is normalised, i.e.,  $ \sum_i s_i \equiv 1$.  This normalised shape,
multiplied by some amplitude $a > 0$, may or may not be present in a noisy data
vector
\begin{equation}
  \label{eq:2}
  \bm{f} = (f_1, f_2, \dots, f_i, \dots, f_{N-1}, f_N)^T
\end{equation}
at known position $i_0$.  The
noise of $\bm{f}$ is given by a Gaussian noise vector with zero mean,
\begin{equation}
  \label{eq:3}
  \bm{n}=(n_1, \dots, n_N)^T \;\text{,}
\end{equation}
that relates\footnote{For example, given a covariance matrix $\bm{\Sigma}$ a
  realisation the noise vector $\bm{n}$ follows from the decomposition
  $\bm{\Sigma} = \bm{C} \bm{C}^T$ and a vector of standard normal uncorrelated
  noise $\tilde{\bm{n}}$ via $\bm{n} = \bm{C} \tilde{\bm{n}}$.} to an
$N\times N$ covariance matrix $\bm{\Sigma}$.  In Eqs.~\eqref{eq:1} to
\eqref{eq:3} and below the superscript $T$ indicates the transpose operation.

Explicitly, we think of $\bm{f}$ in Eq.~\eqref{eq:2} as an observed 1D spectrum
with $i = 1, \dots, N$ indexing the grid of spectral bins, while
$\bm{s}$ in Eq.~\eqref{eq:1} describes the known (or expected) shape of the
emission line to be detected.  The position $i_0$ refers to some well defined
feature in $\bm{s}$, e.g., the position of the line peak or the first moment of
the line profile.  For the problem to be well defined the line width of $\bm{s}$
is assumed to be significantly smaller compared to the length of $\bm{f}$, i.e., 
$M \ll N$.  Lastly, $\bm{n}$ in Eq.~\eqref{eq:3} describes the noise of the
spectrum that results from the detector (read out and dark current) and, for
ground based observations, from the telluric background
\citep[e.g.,][]{Ferruit2010}.  The off-diagonal terms in $\bm{\Sigma}$ that
define the co-variance in $\bm{n}$ are caused by resampling of the detector
pixels, e.g., from multiple exposures and/or rectifying procedures, onto the
final data grid within the data reduction pipeline.  

Without loss of generality we now fix $M \bmod 2 = 1$ and assert that $i_0$
corresponds to the central pixel of the line signal in $\bm{f}$.  This allows us to
introduce a vector of length $N$, in which $\bm{s}$ is embedded and adequately
padded with zeros, i.e., 
\begin{equation}
  \label{eq:5}
  \bm{s}^{(i_0)} =
  \left (
    \bm{0}_{[i_0 - (M-1)/2 - 1]}^T \, ,
    \,
    \bm{s}^T \, ,
    \,
    \bm{0}_{[N - i_0 - (M-1)/2]}^T
  \right )^T  \;\text{,}  
\end{equation}
where $\bm{0}_{[x]}$ denotes a zero vector of length $x$.  So called ``edge
effects'' occur for $i_0 \leq (M-1)/2 - 1$ and $i_0 \geq N - (M-1)/2 $, and
for such lines the detection problem is not well defined.

In order to argue for the presence or absence of a source signal at $i_0$ we now
require a statistic to decide between $\mathcal{H}_0$ and $\mathcal{H}_1$, formally
\begin{align}
  \label{eq:4mod}
  \begin{split}
    \mathcal{H}_0 &\Leftrightarrow \bm{f} = \bm{n}    \\ \text{or} \quad
    \mathcal{H}_1 &\Leftrightarrow \bm{f} = a \cdot \bm{s}^{(i_0)} + \bm{n} \;\text{.}
  \end{split}
\end{align}
It is proven with mathematical rigour that the optimal
decision statistic for rejecting $\mathcal{H}_0$ based on the alternative
$\mathcal{H}_1$ in Eq.~\eqref{eq:4mod} is given by the scalar product
\begin{equation}
  \label{eq:6}
  \mathcal{T}^{(i_0)} = \left (\bm{T}^{(i_0)} \right )^T \bm{f} \;\text{,}
\end{equation}
with the vector
\begin{equation}
  \label{eq:7}
  \bm{T}^{(i_0)} = \alpha_{(i_0)} \, \bm{\Sigma}^{-1} \bm{s}^{(i_0)}
\end{equation}
being the \textit{matched filter} \citep[e.g.,][]{Vio2021}.  $\bm{\Sigma}^{-1}$
in Eq.~\eqref{eq:7} denotes the inverse of the co-variance matrix $\bm{\Sigma}$
and $\alpha_{(i_0)}$ is a normalisation constant that is ideally chosen as
\begin{equation}
  \label{eq:8}
  \alpha_{(i_0)} = 1 / \sqrt{\left (\bm{s}^{(i_0)} \right)^{T} \bm{\Sigma}^{-1}
    \bm{s}^{(i_0)}} \; \text{.}
\end{equation}
The matched filter from Eq.~\eqref{eq:7} normalised according to
Eq.~\eqref{eq:8} then defines $\mathcal{T}^{(i_0)}$ in Eq.~\eqref{eq:6} as
multiples of the standard deviation of the filtered noise.  We then may reject
$\mathcal{H}_0$ and decide on $\mathcal{H}_1$ when $\mathcal{T}^{(i_0)} > \gamma$ and
speak of a line detected at a ``$\gamma \times \sigma$'' significance.
$\mathcal{T}^{(i_0)}$ is therefore also called the \textit{detection significance}.

We note that standard conversion of this ``$\gamma\times\sigma$'' detection
significance into a probability of ruling out $\mathcal{H}_0$
\citep[e.g.,][]{Wall1979}, also known as ``false detection probability'', is
formally only valid if $i_0$ is known a-priori.  While such situations exist in
astronomy, i.e.,  in follow up spectroscopy to detect an emission line of an
already known source \citep[see, e.g.,][]{Loomis2018}, this is not the case in
blind searches for emission lines.  The resulting statistical subtleties for
calculating false detection rates in blind searches on matched filtered data
sets have been extensively studied by \cite{Vio2016}, \cite{Vio2017}, and
\cite{Vio2019} for sub-mm aperture synthesis imaging.  While the details are
involved, generally the real detection significance of an a priori unknown
source is lower than what the standard conversion would suggest.  Perhaps more
intuitively, since $\mathcal{T}^{(i_0)}$ represents the maximised
\textit{signal-to-noise ratio} (SN) of the emission line signal at $i_0$
according to the matched filtering theorem \citep[e.g.,][]{Das1991}, the
hypothesis test can also be regarded as a minimum SN criterion for the detection
of emission lines.  Therefore, $\mathcal{T}^{(i_0)}$ is sometimes also called
the SN of an emission line.

In a blind search for emission lines in $\bm{f}$ we use Eqs.~\eqref{eq:6} --
\eqref{eq:8} to compute the vector
\begin{equation}
  \label{eq:9}
  \bm{\mathcal{T}} = (\mathcal{T}^{(0)}, \dots, \mathcal{T}^{\left((M-1)/2 -1\right)} ,
  \dots, \mathcal{T}^{\left ( N - (M-1)/2 + 1 \right)}, \dots,  \mathcal{T}^{(N)})^T 
\end{equation}
and then select the peak (or peaks, if multiple lines are in the spectrum) that
pass the desired threshold.  Only the elements $\mathcal{T}^{(i)}$ with
$i \geq (M-1)/2 - 1$ and $i \leq N - (M-1)/2 + 1$ can provide a valid test
statistic that is not biased by edge effects.  Explicitly, the computation of
the elements $\mathcal{T}^{(i)}$ of $\bm{\mathcal{T}}$ in Eq.~\eqref{eq:9} is
achieved by replacing $i_0$ with $i$ in Eq.~\eqref{eq:6}:
\begin{equation}
  \label{eq:10}
  \mathcal{T}^{(i)} = \sum_{j = -(M-1)/2}^{(M-1)/2}  T^{(i)}_j f_{i-j} \;\text{.}
\end{equation}
Eq.~\eqref{eq:10} describes the discrete convolution of the of the data vector with the
matched filter.   
The vector $\bm{\mathcal{T}}$ containing the test for each
spectral bin in Eq.~\eqref{eq:9} can thus be computed via
\begin{equation}
  \label{eq:11}
  \bm{\mathcal{T}} = \bm{T} \cdot \bm{f} \; \text{,}
\end{equation}
where $\bm{T}$ is the matrix with elements $T^{(i)}_j$; $j$ indexes the columns
and $i$ indexes the rows of this matrix.  We understand $T^{(i)}_j$ as the $j$-th
element of the matched filter for an emission line source at position $i$ in
$\bm{f}$.  Thus, the matrix $\bm{T}$ collects all the filter profiles for each
$i$, and we can view it as a row stack of the vectors $\bm{T}^{(i)}$ from
Eq.~\eqref{eq:7} normalised by Eq.~\eqref{eq:8}.  Importantly, $\bm{T}$ is a
banded sparse matrix, since $T^{(i)}_j = 0$ for $|j| > (M-1)/2$.
Multiplications with sparse matrices can be computed efficiently numerically,
and high-level interfaces to such efficient implementations exist, e.g., in
SciPy's \texttt{sparse.csr\_matrix} routines \citep{Virtanen2020}.

From now on we assume uncorrelated noise, i.e.,  that the co-variance matrix is diagonal,
\begin{equation}
  \label{eq:12}
  \Sigma_{ij} = 0 \Leftrightarrow i \neq j \qquad \text{and} \qquad \Sigma_{ii}
  = \sigma_i^2 \;\text{.}
\end{equation}
In this case the elements of $\bm{n}$ (Eq.~\ref{eq:3}) are independent (but not
identically distributed) Gaussian random variables.  Moreover, the inverse of
the co-variance matrix is now
\begin{equation}
  \label{eq:13}
  \bm{\Sigma}^{-1} = \mathbb{1} \cdot (\sigma^{-2}_1, \dots, \sigma^{-2}_N)^T \;\text{,}
\end{equation}
with $\mathbb{1}$ being the identity matrix; $\bm{\Sigma}^{-1}$ is thus a
diagonal matrix with $\Sigma^{-1}_{ij} = 0$ and $\Sigma_{ii} = \sigma_i^{-2}$.
According to Eq.~\eqref{eq:7}, Eq.~\eqref{eq:8}, and Eq.~\eqref{eq:13} then the $T^{(i)}_j$ in Eq.~\eqref{eq:10}, i.e.,  the elements of the
matrix $\bm{T}$ in Eq.~\eqref{eq:11}, are
\begin{equation}
  \label{eq:14}
  T^{(i)}_j = \frac{1}{\sqrt{\sum_{j=-k}^{k}
      \frac{s_j^2}{\sigma_{i-j}^2}}} \times \frac{s_j}{\sigma_{i-j}^2} \;\text{.}
\end{equation}
For ease of notation we now replaced $(M-1)/2$ with $k$ in Eq.~\eqref{eq:14}.

Introducing the voxels $F_{x,y,z}$ of a 3D data cuboid $\bm{F}$ and the
associated voxels $\sigma^2_{x,y,y}$ of the variance data cuboid $\bm{\sigma}$,
both of dimensions $X\times Y\times Z$ ($X$ and $Y$ are the spatial dimensions,
while $Z$ is the spectral dimension), as well as the voxels $s_{x',y',z'}$ of
the normalised 3D line source template $\bm{S}$ (with
$\sum_{i,j,k} S_{i,j,k} \equiv 1$) of dimensions $X'\times Y'\times Z'$, with
$X'\ll X$, $Y'\ll Y$, and $Z'\ll Z$, the expressions in Eq.~\eqref{eq:10} and
Eq.~\eqref{eq:14} are then trivially generalised to 3D:
\begin{equation}
  \label{eq:15}
  \mathcal{T}^{(x,y,z)} =   \sum_{i = -d}^{d} \sum_{j =
    -e}^{e} \sum_{k = -f}^{f}  T^{(x,y,z)}_{i,j,k} F_{x-i,y-j,z-k}   \;\text{,}
\end{equation}
with
\begin{equation}
  \label{eq:16}
  T^{(x,y,z)}_{n,m,l} = \frac{1}{\sqrt{\sum_{i=-d}^{+d} \sum_{j=-e}^{e} \sum_{k=-f}^{f}
      \frac{S_{i,j,k}^2}{\sigma_{x-i,y-j,z-k}^2}}} \times \frac{S_{n,m,l}}{\sigma_{x-n,y-m,z-l}^2}  \;\text{.}
\end{equation}
Without loss of generality here $S_{0,0,0}$ was set as the central pixel of the
line source template, such that the summation indices $i$, $j$, and $k$ cover
the whole extend of $\bm{S}$, i.e.,   $2(d,e,f) + 1 = (X', Y', Z')$; in the
following the limits of the summations will be omitted and $\sum_{i,j,k}$ will
be used as a shorthand instead.

Equation~\eqref{eq:15} and Eq.~\eqref{eq:16} define the starting point of the in
\mbox{LSDCat} adopted matched filtering solution for detecting emission lines in
IFS data cuboids; $x$ and $y$ are the spatial coordinates (spaxel coordinates),
and $z$ indexes the cuboid layers in spectral direction (cuboid layer
coordinates).  Eq.~\eqref{eq:16} describes the voxels of the 3D matched filter
for the template $\bm{S}$ at position $x,y,z$ in $\bm{F}$.  Moreover,
$S_{0,0,0}$ equates with the peak of the template profile in \mbox{LSDCat}.  In
principle a reformulation of Eq.~\eqref{eq:15} into a matrix-vector dot-product
like Eq.~\eqref{eq:11} is also possible \citep[see Appendix A.3
of][]{Ramos2011}, but for uncorrelated noise this appears not especially
practical and, hence, is not further pursued here.

In analogy to the vector of Eq.~\eqref{eq:9}, we compute with Eq.~\eqref{eq:15}
the voxels of the cuboid $\bm{\mathcal{T}}$, where peaks above desired detection
threshold, $\gamma$, constitute ``$\gamma \times \sigma$'' detections.  In
LSDCat we dub $\bm{\mathcal{T}}$ as ``SN-cuboid'', in reference to the
SN-maximising characteristic of the matched filter, and we baptise $\gamma$ as
SN-threshold $\mathrm{SN}_\mathrm{thresh}$ in the following.

As in the original implementation of LSDCat \citepalias{Herenz2017}, the
template $\bm{S}$ is optimised for emission line sources, whose spatial and
spectral properties are independent:
\begin{equation}
  \label{eq:17}
  \bm{S} = \bm{S}^\text{spat} \otimes\bm{S}^\text{spec} \;\text{,}
\end{equation}
where $\bm{S}^\text{spat}$ denotes the spatial profile (dimensions
$X' \times Y'$), $\bm{S}_\text{spec}$ denotes the spectral profile (dimension
$Z'$), and $\otimes$ denotes the outer product. The voxels of $\bm{S}$ in
Eq.~\eqref{eq:17} are thus given by
\begin{equation}
  \label{eq:18}
  S_{x,y,z} = S^{\text{spat}}_{x,y} \, S^{\text{spec}}_{z}\;\text{.}
\end{equation}
As we will demonstrate below, this separation provides significant benefits for
the computation of $\bm{\mathcal{T}}$.  We provide a brief overview of the
parametric functions adopted for $S^{\text{spat}}_{x,y}$ and
$S^{\text{spec}}_{z}$ in the current version of LSDCat in \ref{sec:temp}.

The use of uncorrelated variances (Eq.~\ref{eq:12}), which allows to write the
matched filter with the simple expressions of Eq.~\eqref{eq:10} or
Eq.~\eqref{eq:16} for the 1D or 3D case, respectively, is born out of practical
necessity.  As explained above, resampling from detector space to the
data cuboid space certainly introduces co-variance terms (see Fig.~5 in
\citealt{Bacon2017} for a visualisation of the spatial correlation in MUSE), but
the inflation in data volume by storing this information is deemed
computationally not tractable in the current wide-field IFS data reduction
pipelines\footnote{Because of the rapid advancements in computer technology,
  this is not seen as a limitation in the future anymore and the pipeline of the
  planned BlueMUSE IFS \citep{Richard2019} will take into account for
  co-variances due to resampling \citep{Weilbacher2022}.}.  This enforces the
design of our current algorithm to work without co-variances.  A side effect of
ignoring the co-variances in the data reduction is that the absolute values of
$\sigma_{x,y,z}^2$ are underestimated.  Thus, a multiplicative re-scaling of the
variance cuboid is required in order to obtain correct variances peak values for
line detections.  Thus empirical methods for re-scaling the variance cuboid
exist (see Sect.~3.1.1 in \citealt{Herenz2017a}, Sect.~3.2.4 in
\citealt{Urrutia2018} and Sect.~3.1.5 in \citealt{Bacon2017}).  Moreover, all
formal considerations regarding the statistical significance of the detections
from the S/N peaks are only valid as long as $\bm{S}$ matches the real emission
line signal in the cuboid, but in reality a diversity of emission line shapes
will be encountered; we will deal explicitly with the effects of such mismatches
between filter and emission line signals in the context of MUSE survey data in
Sect.~\ref{sec:rs}.


\subsection{The original LSDCat Ansatz and its flaw}
\label{sec:classic}

The evaluation of Eq.~\eqref{eq:15} with $T^{(x,y,z)}_{i,j,k}$ from
Eq.~\eqref{eq:16} is computationally challenging due to the large dimensions of
the MUSE data (typically
$X \times Y \times Z \simeq 300 \times 300 \times 3800 = 3.42 \times 10^8$
voxels, but some studies require mosaiced ``super-cuboids'' with significantly
larger spatial dimensions, see, e.g., \citealt{SanchezAlmeida2022}).  In the
original implementation of LSDCat we thus computed instead
\begin{equation}
  \label{eq:19}
  \mathcal{T}_\text{classic}^{(x,y,z)} = \frac{\sum_{i,j,k}
    S_{i,j,k} F_{x-i,y-j,z-k}}{\sqrt{\sum_{i,j,k} S^2_{i,j,k}
      \sigma^2_{x-i,y-j,z-k}}} \; \text{.}
\end{equation}
With the help of Eq.~\eqref{eq:18} we could separate Eq.~\eqref{eq:19} into two
simple convolution operations, that can be computed quickly via Fast-Fourier
transformation for $\bm{S}^\text{spat}$ and sparse-matrix multiplication for
$\bm{S}^\text{spec}$.

We can justify the use of Eq.~\eqref{eq:19}, by noting that Eq.~\eqref{eq:15}
with the formal matched filter $T^{(x,y,z)}_{i,j,k}$ according to
Eq.~\eqref{eq:16} is well approximated by Eq.~\eqref{eq:19}, provided that
$\sigma^2_{x,y,z}$ does not vary strongly within the limits of summation.  We
demonstrate this for the 1D case, i.e.,  Eq.~\eqref{eq:10} with $T^{(i)}_j$
according to Eq.~\eqref{eq:14}.  Setting $\sigma_i = \sigma_\mathrm{const}$ in
Eq.~\eqref{eq:10} and Eq.~\eqref{eq:14} we find
\begin{equation}
  \label{eq:20}
  \mathcal{T}^{(i)} \approx \frac{\frac{1}{\sigma_\text{const.}^2} \sum_j s_j
    f_{i-j}}{\frac{1}{\sigma_\text{const.}} \sqrt{\sum_j s_j^2}} = \frac{\sum_j s_j
    f_{i-j}}{\sigma_\text{const} \sqrt{\sum_j s_j^2}} \;\text{.}
\end{equation}
And for the 1D-case of Eq.~\eqref{eq:19},
\begin{equation}
  \label{eq:21}
  \mathcal{T}^{(i)}_\text{classic} = \frac{\sum_j s_j f_{i-j}}{\sqrt{\sum_j s_j^2
    \sigma_{i-j}^2 }} \; \text{,}
\end{equation}
we have with $\sigma_i := \sigma_\mathrm{const}$
\begin{equation}
  \label{eq:22}
  \mathcal{T}^{(i)}_\text{classic} \approx \frac{\sum_j s_j
    f_{i-j}}{\sigma_\text{const} \sqrt{\sum_j s_j^2}} \;\text{,}
\end{equation}
which is identical to Eq.~\eqref{eq:20}.  However, when there are strong
variations of the variance values within the summation limits, then
Eq.~\eqref{eq:21}, or Eq.~\eqref{eq:19} for the 3D case, does not produce
correct results.   In this case we expect the detection significances
to be biased towards lower values.

We illustrate the difference between $\mathcal{T}^{(i)}_\mathrm{classic}$
(Eq.~\ref{eq:21}) and $\mathcal{T}^{(i)}$ (Eq.~\ref{eq:10} with Eq.~\ref{eq:14})
in Figure~\ref{fig:mfdemo} for an example of increased variances in the wings of
an idealised 1D emission line signal that is recovered with its matched filter.
Here the signal is assumed to be a perfect 1D Gaussian (dispersion = 4 bins;
peak amplitude = 2.5). For each bin a unit variance is assumed, except for 4
bins leftwards and rightwards of the peak position, where a variance of 6.25 has
been assumed; the total length of $\bm{f}$ is 81 bins and the peak of the line
is placed at $i_0 = 40$.  The ground truth according to this set-up is shown in
the upper left panel of Figure~\ref{fig:mfdemo}.  A particular realisation of
this set-up is shown in the upper right panel of this figure, where the line
signal clearly got buried in noise.  In the bottom two panels of
Figure~\ref{fig:mfdemo} we contrast the resulting vectors $\bm{\mathcal{T}}$
(Eq.~\ref{eq:9}) where the $\mathcal{T}^{(i)}$ have been calculated with the
classic formalism (Eq.~\ref{eq:21}; bottom left panel) and the correct matched
filter formalism (Eq.~\ref{eq:10} with Eq.~\ref{eq:14}; bottom right panel).
These results are from $10^3$ Monte-Carlo realisations of the set-up.  It can be
seen that the peak values $\mathcal{T}^{(i_0)}$, which is the relevant decisions
statistic for the presence or absence of a source, are clearly lower for the
simplified formulae ($\mathcal{T}^{(i_0)}_\mathrm{classic}$ -- mean: 3.3;
10-/90-percentile: 2.0/4.5) then for the proper matched-filter (mean: 4.8;
10-/90-percentile: 3.3/6.4).  Or, putting it differently, in 72\% of all
detection attempts with $S/N_\mathrm{thresh} = 4 $ the line would have been
rejected with the classic LSDCat Ansatz, whereas this happens only for 12\% of
all attempts with the same threshold in the correct formalism.

While the idealised set-up in Figure~\ref{fig:mfdemo} is contrived for
illustrative purposes, spectrally rapidly varying variances appear in the red
part of the optical spectral range (R and I-band) due to the air glow background
\citep[see Sect.~4 of][]{Noll2012}.  Here predominantly the molecular Meinel
OH-bands are contributing, which appear like combs in 1D spectral
representations \citep[see also Figures in][]{Hanuschik2003} and the spectral distance
between individual air glow lines is only a few \AA{}.  The result is that
emission line signals which are overlapping with the OH-lines will not be
optimally recovered by the classic algorithm.

\begin{figure}
  \centering
  \includegraphics[width=0.5\textwidth]{./figures/matched_filter_test_stat.pdf}
  \caption
   
    Illustration of the difference in SN between a simplified matched filter
    that assumes variances to be constant over the filter profile and the
    correct matched filter that correctly accounts for uncorrelated variances.
    the The \textit{top left} panel displays the ground truth: an emission line
    signal $\bm{s}$ (black line) and a corresponding variance vector (we show
    the square-root $\sigma$ as grey line).  The \textit{top right} panel shows
    a Monte-Carlo realisation of a vector $\bm{f}$ resulting from this set-up.
    The \textit{bottom left} panel shows the SN spectrum from the approximate
    matched filter ($\mathcal{T}^{(i)}_\mathrm{classic}$; Eq.~\ref{eq:21}),
    whereas the \textit{bottom right} panel shows the SN spectrum $\bm{T}$ from
    the correct matched ($\mathcal{T}^{(i)}$; Eq.~\ref{eq:10} with
    Eq.~\ref{eq:14}).  In both bottom panels we show the average from 1000
    realisations (\textit{black line}), the 90 percentile (\textit{grey shaded
      region}) and the SN spectrum of the example realisation from the top right
    panel (dotted line).  }
  \label{fig:mfdemo}
\end{figure}


\subsection{The improved Ansatz (LSDCat\,2.0)}
\label{sec:modified}

Our goal for the improved algorithm is now to account for spectrally varying
background in IFS data cuboids while maintaining the computational simplicity of
the original algorithm.  This can be achieved, when we assume that the variances
are not varying as a function of position, i.e.,  when we set
\begin{equation}
  \label{eq:23}
  \sigma_{x,y,z}^2 := \sigma_{z}^2 \;\text{.}
\end{equation}
Under this condition we can separate the template in Eq.~\eqref{eq:16} using
Eq.~\eqref{eq:18}, i.e., 
\begin{equation}
  \label{eq:24}
  T^{(x,y,z)}_{i,j,k} = T_k^{\{z\}} T_{ij}^{\{z\}} \;\text{,}
\end{equation}
with
\begin{equation}
  \label{eq:25}
  T_k^{\{z\}} = \frac{1}{\sqrt{\sum_n
      \frac{\left (S_n^{\{z\}} \right)^2}{\sigma^2_{z-n}}}} \times \frac{S^{\{z\}}_k}{\sigma_{z-k}^2}
\end{equation}
and
\begin{equation}
  \label{eq:26}
  T_{ij}^{\{z\}} = \frac{S_{ij}^{\{z\}}}{\sqrt{\sum_{lm}
      \left (S_{lm}^{\{z\}} \right )^2}}  \;\text{.}
\end{equation}
The superscript ${\{z\}}$ in Eqs.~\eqref{eq:24} -- \eqref{eq:26} is meant to
indicate the implicit dependence of our templates $\bm{S}^\text{spat}$ and
$\bm{S}^\text{spec}$ on wavelength and thus on the index $z$ in spectral
direction (see \ref{sec:temp} for the rationale).  We note that the spatial
filter in Eq.~\eqref{eq:26} is just the re-normalised spatial template and the
spectral filter in Eq.~\eqref{eq:25} is identical to the spectral filter in the
1D case introduced in Sect.~\ref{sec:mf} (Eq.~\ref{eq:14}).  Plugging
Eq.~\eqref{eq:24} into Eq.~\eqref{eq:15} we find that the voxels of the
SN-cuboid can be computed via
\begin{equation}
  \label{eq:27}
  \mathcal{T}^{(x,y,z)} = \sum_k T_k^{\{z\}} \left ( \sum_{ij} T_{ij}^{\{z\}}
    F_{x-i,y-j,z-k}  \right) \;\text{.}
\end{equation}

The numerical computation of Eq.~\eqref{eq:27} can be carried out efficiently.
First, we evaluate Eq.~\eqref{eq:26} for all $z$ for the parametric models for
$S_{ij}^{\{z\}}$ (Eq.~\ref{eq:43} or Eq.~\ref{eq:44} in \ref{sec:temp}) and we
then use the discrete fast Fourier transform to calculate the convolution sum in
brackets for each $z$.  This operation can be trivially parallelised by
partitioning $\bm{F}$ along the spectral axis into as sub-cuboids for each CPU
core.  For the outer sum we first pre-compute all $T_k^{\{z\}}$ according to
Eq.~\eqref{eq:25} for the given spectral template (Eq.~\ref{eq:46} in
\ref{sec:temp}).  Then we see the similarity between the outer sum and
Eq.~\eqref{eq:10}, i.e.,  we can compute for each spaxel a (sparse) matrix -
vector product according to Eq.~\eqref{eq:11}.  Again, this step is trivially
parallelised, now by partitioning the data cuboid over the spatial axes.

Above simplification only works under the provision that Eq.~\eqref{eq:23}
holds, i.e.,  that the variances do not vary as a function of position.  We found
that this is indeed the case for MUSE data cuboids obtained from observations
that follow the recommended dithering and rotation strategies.  In fact,
\cite{Herenz2017a}, when compiling the first catalogue of emission line galaxies
in the MUSE-Wide survey, already calculated an 1D empirical variance spectrum
$\sigma^2_z$, that was then blown up to a variance cuboid via
$\sigma^2_{x,y,z} = N^\mathrm{exp}_{x,y}/N_\mathrm{max} \cdot \sigma^2_z$, where
$N^\mathrm{exp}_{x,y}$ denotes the exposure count map that counts the average
number of exposures that contribute to each spaxel and $N_\mathrm{max}$ denotes
the maximum number of exposures, i.e., 
$N_\mathrm{max} = \max_{x,y} N^\mathrm{exp}_{x,y}$.  The procedure to calculate
such an empirical variance spectrum was then refined by \cite{Urrutia2018}, but
again $N^\mathrm{exp}_{x,y}$ was used for rescaling.  Significant differences of
$N^\mathrm{exp}_{x,y}$ with respect to $N_\mathrm{max}$ occur near the edges of
the field of view, but here edge effects render the decision statistic from the
matched filter unreliable anyway.

\subsection{Verification of the improved method}
\label{sec:verify}

We verified the expected improvement of detection significances with the method
presented in Sect.~\ref{sec:modified} by a re-analysis of the publicly available
data of MUSE-Wide survey \citep{Urrutia2018}.  To this aim we re-processed all
released data cuboids in the same way as in the \citeauthor{Urrutia2018} study,
i.e., , we used the template parameters provided in Table~2 of \cite{Urrutia2018},
but now with the \mbox{LSDCat2.0} routines.  For each detected emission line we
then compared the peak S/N values between the original Ansatz
(Sect.~\ref{sec:classic}) and the improved Ansatz (Sect.~\ref{sec:modified}) for
all 3057 emission lines in the \citeauthor{Urrutia2018} catalogue.  The result
of this exercise is presented in Figure~\ref{fig:comp}, where we show for each
catalogued detection the relative difference between the new and the old method,
$\mathrm{SN}_{2.0}/\mathrm{SN}_\mathrm{classic} - 1$, as a function of detection
wavelength.  Here $\mathrm{SN}_{2.0}$ refers to the peak values from
Eq.~\eqref{eq:27}, while $\mathrm{SN}_\mathrm{classic}$ refers to the peak
values according to Eq.~\eqref{eq:19}.  We also inset an arbitrarily scaled
effective variance spectrum into Figure~\ref{fig:comp} to illustrate the
position and amplitude of the air-glow lines in the MUSE spectral range.

As can be seen from Figure~\ref{fig:comp}, the improved Ansatz indeed boosts the
SN of many detected emission lines especially in the forest of air glow lines.
The overall improvement also becomes apparent in the histogram of relative SN
differences shown in Figure~\ref{fig:hist}.  There we also separated the counts
in each bin by their spectral position in the data cuboid, and again it can be
seen how, especially in the red part of the spectrum, a significant number of
sources experiences a boost in their detection significances.  In numbers: 1892
detections have $\mathrm{SN}_{2.0} > \mathrm{SN}_\mathrm{classic}$, whereas for
1165 lines $\mathrm{SN}_{2.0} < \mathrm{SN}_\mathrm{classic}$.  However, for the
vast majority of detections where the improved algorithm provides a lower SN,
the difference is marginal, whereas the boost for the former can be substantial.
For example, while only 99 detections have
$\mathrm{SN}_{2.0} < 0.95 \times \mathrm{SN}_\mathrm{classic}$, 722 show
$\mathrm{SN}_{2.0} > 1.05 \times \mathrm{SN}_\mathrm{classic}$.

Nevertheless, we also find in this comparison that ten emission lines show a
peak value that is marginally lower than the original detection threshold.  All
these catalogued lines are overlapping with air-glow lines, one even with the
extremely strong [\ion{O}{I}] $\lambda5577$ line.  For those lines the
$\mathrm{SN}_\mathrm{classic}$ value got artificially boosted due to a strong
positive noise contribution that gets adequately suppressed in the new
formalism.  Eight of those ten lines are secondary lines of emission line
galaxies that have been detected in other, stronger, lines, i.e.,  these galaxies
are not lost from the catalogue.  The remaining two single-line detections (ID
129026144 and 146018244) were supposed Ly$\alpha$ emitters, but these low SN
detections were assigned also with the lowest confidence (i.e.,  $\approx$50\%
probability of being false detections according to \citealt{Urrutia2018}).
Thus, our new algorithm does not remove significant detections from the previous
catalogue, but significantly improves the detection significances of emission
lines that are overlapping with air-glow lines.  Hence, this little experiment
verifies indeed that the new algorithm performs as expected.

We limited ourselves here to only check emission lines that were already
catalogued.  Nevertheless, because of the overall improvement of low-SN
detections within the forest of the OH-bands we foresee that the
\mbox{LSDCat2.0} algorithm also uncovers new lines in those spectral regions.
However, identifying and classifying new detections still requires manual
inspection and verification by a team of experts (see \citealt{Inami2017} and
\citealt{Urrutia2018} for descriptions of this process).  Such a complete
re-analysis of the MUSE-Wide dataset, clearly beyond the scope of the present
work, may be desired for subsequent data MUSE-Wide releases.

\begin{figure*}[t!]
  \centering
  \includegraphics[width=0.85\textwidth, trim=10 10 20
  10,clip=True]{./figures/plot_compare_v1_v2.pdf}
  \caption{ Comparison between of the SNs obtained with the improved filter,
    $\mathrm{S/N}_{2.0}$ (Eq.~\ref{eq:27}), to those obtained with the classic
    filter, $\mathrm{S/N}_\mathrm{classic}$ (Eq.~\ref{eq:19}), for all
    catalogued emission lines of the MUSE-Wide data release by
    \cite{Urrutia2018}.  The relative difference of $\mathrm{S/N}_{2.0}$ with
    respect to $\mathrm{S/N}_\mathrm{classic}$ is shown as a function of peak
    detection wavelength $\lambda_\mathrm{S/N_\mathrm{peak}}$, and the absolute
    $\mathrm{S/N}_{2.0}$ is colour coded as indicated by the colourbar
    (detections with $\mathrm{S/N}_{2.0} > 40$ are not assigned to a different
    colour than $\mathrm{S/N}_{2.0} = 40$ detections).  The 10 emission lines
    for which $\mathrm{S/N}_{2.0}$ is below the detection threshold of the
    original catalogue are marked by red crosses.  An arbitrarily scaled
    variance spectrum is plotted as a grey line to indicate the spectral
    variation of the telluric background. }
  \label{fig:comp}
\end{figure*}

\begin{figure}[t]
  \centering
  \includegraphics[width=0.45\textwidth,trim=20 10 20
  0,clip=True]{./figures/plot_compare_v1_v2_hist.pdf}
  \caption{Stacked histogram (bin width = 0.05) of the relative differences of
    $\mathrm{S/N}_{2.0}$ with respect to $\mathrm{S/N}_\mathrm{classic}$
    (cf. Fig.~\ref{fig:comp}); emission lines at
    $\lambda_\mathrm{SN_\mathrm{peak}} \leq 7000$\,\AA{} are counted in blue,
    those at $\lambda_\mathrm{SN_\mathrm{peak}} 7000$\,\AA{} are counted in
    orange.  }
  \label{fig:hist}
\end{figure}



\section{The selection function of the 3D matched filter}
\label{sec:selfun}

A selection function, $f_C$ (sometimes also referred to as completeness
function, e.g., in \citealt{Caditz2016}), encodes the probability that a source
with given attributes will enter our catalogue.  Selection functions are needed
whenever we want to construct a model of reality given the catalogue data
\citep[see][]{Rix2021}.  Relevant applications in the context of surveys for
emission line objects with IFS are, e.g., galaxy luminosity functions or
modelling the galaxy power spectrum via correlation analysis
\citep[e.g.,][]{Blanc2011,Herenz2019,HerreroAlonso2021,Zhang2021}.  In terms of emission
line source detection in IFS data cuboids the most important parameters are the
wavelength, $\lambda$, of the detection and the flux of the emission line,
$F_\mathrm{line}$.

Processing the flux data cuboids with a 3D matched filter (Eq.~\ref{eq:15} or,
more specifically, Eq.~\ref{eq:27}), can be understood as a mapping of
$F_\mathrm{line}$ to SN space.  As long as the observations are background
limited this mapping is proportional to the emission line flux:
$S/N \propto F_\mathrm{line}$ (we recall from Sect.~\ref{sec:mf} that
$\mathcal{T}^{x,y,z}$ can be understood as the maximised SN for a emission line
at position $x,y,z$ in the data cuboid).  The proportionality factor depends on
the noise properties encountered at the spatial ($x,y$) and spectral ($\lambda$
or $z$) position of the emission line source in the data cuboid and on the level
of congruence between template $\bm{S}$ and actual emission line source.

For the construction of the filter we explicitly assumed in
Sect.~\ref{sec:modified} that the noise properties are not an implicit function
of position ($x,y$) and that the effect of varying depth over the field of view
can be accounted for by rescaling.  Moreover, because of the randomness of the
noise, the value of SN is stochastic and thus we can abbreviate the mapping as
\begin{equation}
  \label{eq:28}
  E[S/N_\mathrm{\lambda}] = C(\lambda) \times F_\mathrm{line}\;\text{,}
\end{equation}
where $E[S/N_\mathrm{\lambda}]$ denotes the expectation value of
$S/N_\mathrm{\lambda}$ and $C(\lambda)$ is in units of inverse flux.  We also
assumed that the noise is normally distributed and uncorrelated, thus the
distribution of $S/N_\mathrm{\lambda}$ for a given $F_\mathrm{line}$,
$g(S/N_\lambda | F_\mathrm{line})$, is also normal and with variance
$V[S/N_\mathrm{\lambda}] = 1$, i.e., 
\begin{equation}
  \label{eq:4}
  g(S/N_\lambda | F_\mathrm{line}) = \frac{1}{\sqrt{2}} \exp \left ( - \frac{(S/N_\lambda -
      C(\lambda) F_\mathrm{line})^2}{2} \right ) \;\mathrm{.}
\end{equation}
Now the selection function,
$f_C(F_\mathrm{line}, \lambda | \mathrm{S/N}_\mathrm{thresh}) \in [0,1]$,
encodes the cumulative probability that $F_\mathrm{line}$ will result in
$S/N_\mathrm{\lambda}(F_\mathrm{line}) > \mathrm{S/N}_\mathrm{thresh}$, i.e., 
\begin{equation}
  \label{eq:49}
  f_C(F_\mathrm{line}, \lambda | \mathrm{S/N}_\mathrm{thresh}) =
  \int_{\mathrm{S/N}_\mathrm{thresh}}^\infty g(S/N_\lambda | F_\mathrm{line})
  \,\mathrm{d}(S/N_\mathrm{\lambda}) 
\end{equation}
for a given a detection threshold $\mathrm{S/N}_\mathrm{thresh}$.  Using the
error-function
$\erf (x) = \frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}\,\mathrm{d}t = 1 -
\frac{2}{\sqrt{\pi}} \int_x^\infty e^{-t^2} \, \mathrm{d} t$, we can express
Eq.~\eqref{eq:49} as
\begin{multline}
  \label{eq:29}
  f_C(F_\mathrm{line}, \lambda | \mathrm{S/N}_\mathrm{thresh}) 
  = \\\frac{1}{2}
  \left [  1 + \erf \left ( \frac{C(\lambda) \times F_\mathrm{line} -
        \mathrm{S/N}_\mathrm{thresh}}{\sqrt{2}}  \right )  \right ] \;\text{.}
\end{multline}
As can be seen, the selection function is completely determined if the
proportionality factor $C(\lambda)$ is known.

Due to the symmetry of the normal distribution, Eq.~\eqref{eq:29} is also
symmetric around $F_\mathrm{line}$, and this flux thus equates with the 50\,\%
completeness limit which we now denote $F_{50}$, i.e., 
$f_C(F_{50}, \lambda | \mathrm{S/N}_\mathrm{thresh}) \equiv 0.5$ and
\begin{equation}
  \label{eq:30}
  F_{50} = \frac{\mathrm{S/N}_\mathrm{thresh}}{C(\lambda)}
  \;\text{.}
\end{equation}


\subsection{The idealised selection function}
\label{sec:ideal-select-funct}

We define the idealised selection function as the selection function for an
emission line source in the data cuboid that is perfectly matched by the source
template $\bm{S}$ (Eq.~\ref{eq:17} and Eq.~\ref{eq:18}).  For the calculation of
$C(\lambda)$ we calculate the response of the filter for such a source at some
position $x',y',z'$ in the data cuboid $\bm{F}$:
\begin{equation}
  \label{eq:31}
  F_{x,y,z} = A \cdot S_{x-x',y-y'}^\mathrm{spat} S^\mathrm{spec}_{z-z'} + n_{x,y,z}\;\text{,}
\end{equation}
where $A$ is the amplitude at which the normalised source profile and
$n_{x,y,z}$ describes a source-free noise cuboid.  Then the matched filter of
Eq.~\eqref{eq:27} with the template according to Eq.~\eqref{eq:25} and
Eq.~\eqref{eq:26} has the expectation value
\begin{equation}
  \label{eq:32}
  E_\mathcal{I} \left [ \mathcal{T}^{(x',y',z')} \right ] = A \cdot \sqrt{ \sum_{ij} \left ( S_{ij}^{\{z' \}}
    \right )^2 } \cdot \sqrt{\sum_{k}  \frac{\left ( S_k^{\{z' \}}  \right
      )^2}{\sigma_{z'-k}^2} } \; \text{,}
\end{equation}
where we used $E[n_{x,y,z}] = 0$.  In anticipation of the discussion in
Sect.~\ref{sec:rs} below we wrote $E_\mathcal{I}$ in Eq.~\eqref{eq:32} to
signify that this is the expectation value for the idealised scenario.

A voxel in the flux data cuboid $F_{x,y,z}$ is assumed to encode flux density.
For example, MUSE DRS provides fluxes in erg\,s$^{-1}$cm$^{-2}$\AA{}.   The
amplitude A in Eq.~\eqref{eq:31} then relates to the emission line flux
$F_\mathrm{line}$ in erg\,s$^{-1}$cm$^{-2}$ via
\begin{equation}
  \label{eq:33}
A = F_\mathrm{line}/\Delta \lambda \;\text{,}  
\end{equation}
where $\Delta \lambda$ is the native wavelength sampling of the data cuboid
(default $\Delta \lambda = 1.25$\,\AA{} for MUSE).  The wavelength $\lambda$
relates linearly to the spectral coordinate $z$ (Eq.~\ref{eq:47}) and we
abbreviate it by $\lambda[z]$.  We now equate Eq.~\eqref{eq:28} with Eq.~\eqref{eq:32}
to find
\begin{equation}
  \label{eq:34}
  C_\mathcal{I}(\lambda [z]) = \frac{1}{\Delta\lambda}  \cdot 
  \sqrt{ \sum_{ij} \left ( S_{ij}^{\{z \}}
    \right )^2 } \cdot \sqrt{\sum_{k}  \frac{\left ( S_k^{\{z \}}  \right
      )^2}{\sigma_{z-k}^2} } \;\text{,}
\end{equation}
where again the subscript $\mathcal{I}$ signifies that this is the
proportionality constant for the idealised selection function.

Equation~\eqref{eq:29} with $C(\lambda[z])$ from Eq.~\eqref{eq:34} provides an
expression for the idealised selection function for the line search with some
template $\bm{S}$ (Eq.~\ref{eq:17}) given an effective variance spectrum
$\sigma^2_z$.  


\subsection{Verification of the analytical expression for the idealised
  selection function}
\label{sec:verif-analyt-expr}

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.45\textwidth,trim=10 0 3 3
  ,clip=True]{./figures/plot_compare_exp_theo.pdf} \vspace{-0.5em}
  \caption{Comparison between the idealised selection function according to
    Eq.~\eqref{eq:29} with $C(\lambda)$ from Eq.~\eqref{eq:34} ($f_C$,
    \textit{top panel}) to the idealised selection function from a source
    insertion and recovery experiment ($f_C^\mathrm{exp}$, \textit{centre
      panel}) for a $z > 6$ LAE search with LSDCat2.0 in the MXDF \citep{Bacon2022}
    using the template given in
    Table~\ref{tab:sn2} and a detection threshold of
    $S/N_\mathrm{thresh} = 6.41$. The \textit{bottom panel} shows the difference
    $f_C - f_C^\mathrm{exp}$. }
  \label{fig:2dsel}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.5\textwidth]{./figures/plot_lsd20pred_vs_f50_pq.pdf}
  \caption{Analytical 50\% completeness estimate in comparison to the 50\%
    completeness achieved in a source insertion and recovery experiment,
    $F_{50}^\mathrm{exp}$, for the LSDCat2.0 line search at
    $\lambda > 8500$\,\AA{} in the MXDF with $S/N_\mathrm{thresh} = 6.41$ under
    the assumption of idealised sources described by the 3D template in
    Table~\ref{tab:sn2}.  The \textit{top panel} displays both the analytical
    $\log F_{50}$ (Eq.~\ref{eq:30} with $C(\lambda)$ from Eq.~\ref{eq:34};
    \textit{dashed black line}) and the $\log F_{50}^\mathrm{exp}$ from the
    numerical experiment (\textit{red line}), whereas the \textit{bottom panel}
    displays the difference $\Delta = \log F_{50} - \log F_{50}^\mathrm{exp}$. }
  \label{fig:f50test}
\end{figure}

\begin{figure*}
  \centering
  \includegraphics[width=\textwidth,trim=30 11 45 0,clip=True]{./figures/plot_lsd20_pres_vs_ins_crvs_pq.pdf}
  \caption{Analytical completeness curves (\textit{dashed black lines};
    Eq.~\ref{eq:29} with $C(\lambda)$ from Eq.~\ref{eq:34}) in comparison to
    the completeness curves from a source insertion and recovery experiment (red
    lines) for the \mbox{LSDCat2.0} line search in the MXDF (see text and
    captions to Figures~\ref{fig:2dsel} and \ref{fig:f50test} for details) at
    four representative wavelength layers.  In each panel we print the
    difference $ \Delta = \log F_{50}^\mathrm{ana} - \log F_{50}^\mathrm{exp}$
    (cf. Fig.~\ref{fig:f50test}) and we inset a zoom-in of the extracted average
    noise spectrum around the layers of interest.}
  \label{fig:f50crvs}
\end{figure*}

We performed a source insertion and recovery experiment to check whether
Eq.~\eqref{eq:29} with $C(\lambda)$ from Eq.~\eqref{eq:34} can indeed predict
the idealised selection function of the algorithm.  For this experiment we used
the recently released data cuboid of the MUSE eXtreme Deep Field (MXDF) by Bacon
et al. (2022).  This dataset is not only the deepest blank field ever obtained
with MUSE ($t_\mathrm{exp} = 141$\,h), but it also used a non-standard observing
strategy to reduce potential systematics due to the rectangular alignment of the
image slicer and the instrument field splitter.  The result is that the final
combined field of view within the data cuboid is circular.  As groundwork for
future population studies of unprecedented faint Ly$\alpha$ emitters in this
dataset we constructed an emission line emitter catalogue with LSDCat For
reference we list the provide the template parameters used in this search in
\ref{sec:templ-param-used} (Table~\ref{tab:sn2}).  As detection threshold we
chose $S/N_\mathrm{thresh} = 6.41$; the rationale for this value will be
detailed in Herenz et al. (in prep.).  For the demonstration here we only
consider the extreme red of the MUSE spectral range, where the density and
amplitude of the telluric background emission lines is highest.

For the recovery experiment we implanted emission line sources into the MXDF
that are exactly described by this template.  We used a regular grid of 70
insertion positions and we considered only the central region of $> 140$\,h
depth \citep[see Fig.~3 in][]{Bacon2022}.  The distance between neighbouring sources
was 6\arcsec{}.  At each wavelength layer $\lambda[z]$ we then inserted the
artificial line sources with fluxes from
$\log F_\mathrm{line}^\mathrm{ins}$[erg\,s$^{-1}$cm$^{-2}$]$ = - 18.6$ to
$\log F_\mathrm{line}^\mathrm{ins}$[erg\,s$^{-1}$cm$^{-2}$]$= -17.1$ at an
increment of $3 \times 10^{-20}$erg\,s$^{-1}$cm$^{-2}$; in each insertion
experiment all sources have the same flux.  After each insertion we ran the full
source detection chain that was used for the construction of the catalogue,
i.e.,  we subtracted a running median to remove continuum emission and we ran the
3D matched-filter with the template given in Table~\ref{tab:sn2}.  We then
counted the number of recovered lines above $S/N_\mathrm{thresh}$ which, divided
by the number of inserted sources, provides us with the experimentally
determined completeness
$f_C^\mathrm{exp.}(F_\mathrm{line}^\mathrm{ins}, \lambda[z])$ for each
($F_\mathrm{line}^\mathrm{ins}, \lambda[z]$)-pair.  We excluded detections of
real sources in our tally by excluding insertion positions that overlapped with
detections of real emission line sources in the original search, i.e.,  in some
layers the denominator of inserted sources was less than 70 for the calculation
of $f_C^\mathrm{exp.}(F_\mathrm{line}^\mathrm{ins}, \lambda[z])$.

In Figure~\ref{fig:2dsel} we show a comparison between the analytic prediction
for the idealised case, i.e.,
$f_C(\lambda,F_\mathrm{line} | S/N_\mathrm{thresh} = 6.41)$ according to
Eq.~\eqref{eq:29} with $C(\lambda)$ from Eq.~\eqref{eq:34} for the template
parameters given in Table~\ref{tab:sn2}, and the experimentally determined
selection function, $f_C^\mathrm{exp}$, from the above described insertion and
recovery procedure. To compute the analytic prediction we evaluated
Eq.~\eqref{eq:34} with an averaged variance spectrum taken from 100 spaxels of
the variance cuboid in the centre of the field of view.  The visual comparison
between the analytical prediction (top panel of Figure~\ref{fig:2dsel}) with the
experimental result (middle panel of Figure~\ref{fig:2dsel}) does not reveal
large discrepancies.  Only the sampling in $\log F$ of the numerical experiment
becomes apparent as pixelated structure in the middle panel.  Obviously, the
analytic prediction can be evaluated on a much finer flux grid than what would
be feasible for an insertion and recovery experiment.  The impression of
congruity between prediction and experiment is verified in the bottom panel of
Fig.~\ref{fig:2dsel}, where we show the difference between the $f_C$ and
$f_C^\mathrm{exp}$.  Only very few $(F_\mathrm{line},\lambda)$-bins show a
difference greater than 20\%, and generally the experimental and analytical
curves are close to each other.

The congruence between prediction and experiment can also be appreciated in
Figure~\ref{fig:f50test}, where we compare the analytic 50\% completeness
(Eq.~\ref{eq:30}) to the 50\% completeness from the experiment.  The latter was
calculated via linear interpolation using the two flux insertion values for
which $f_C^\mathrm{exp}$ had the smallest absolute difference to $0.5$.  Again,
both curves are nearly identical, with the average
$\Delta = \log F_{50} - \log F_{50}^\mathrm{exp}$ being $0.01$.  Inspecting
Figures~\ref{fig:2dsel} and \ref{fig:f50test} in detail reveals that the
experiment results in slightly deeper flux limits, but the discrepancies are
marginal and dwarfed in comparison to the uncertainties on the selection
function due to mismatches between filter template and real sources (see
Sect.~\ref{sec:rs} below).  Importantly, the fluctuations on $\Delta$ are random
and do not correlate with features in the sky background.

Lastly, we compare in Figure~\ref{fig:f50crvs} the analytic
$f_C(F_\mathrm{line}, \lambda)$ curves with their experimental counterparts for
four different wavelength layers.  These layers are chosen to be representative
for the typical situations in the background noise spectrum.  The first panel
(from the left) demonstrates a wavelength layer where the background noise is
lower, but sharp spikes rise in redward and blueward layers; this situation is
akin to the example discussed in Sect.~\ref{sec:classic}
(Figure~\ref{fig:mfdemo}).  The second and fourth panel show layers that
corresponds to the wings of a [OH]-line, and the third panel shows a layer where
the surrounding noise spectrum is not varying strongly.  Again, all experimental
curves show an good level of agreement with the analytic predictions.



\subsection{Towards a realistic selection function}
\label{sec:rs}


\begin{figure}[t]
  \centering \includegraphics[width=0.45\textwidth,trim=0 10 0
  5,clip=True]{./figures/plot_xi_ps_sn2.pdf}
  \caption{Example for the ratio $C_{\widetilde{S}}/C_\mathcal{I}$ according
    to Eq.~\eqref{eq:39} for the search with a 3D Gaussian template
    (Eq.~\ref{eq:36} with $S_{xy}$ from Eq.~\ref{eq:43} and $S_z$ from
    Eq.~\ref{eq:46}) search for Lyman $\alpha$ emitting galaxies in the MXDF
    (Sect.~\ref{sec:ideal-select-funct}, Table~\ref{tab:sn2}, and
    $C_\mathcal{I}(\lambda)$ from Eq.~\ref{eq:34}) and point sources described
    by the Moffat model (Eq.~\ref{eq:44}) reconstructed with \texttt{muse-psfr}
    (Table~\ref{tab:mof}).  The \textit{top panel} shows
    $C_{\widetilde{S}}/C_\mathcal{I}$ for point sources described by this Moffat
    template (\textit{solid line}) and by a Gaussian template of same FWHM as
    the Moffat template (\textit{dashed line}).  The \textit{bottom panel} shows
    the FWHM difference between the \texttt{muse-psfr} reconstructed point
    spread function and the search template.  The dotted line in both panels
    indicates $\lambda = 8700$\,\AA{} which is the wavelength used for
    example shown in Figure~\ref{fig:xhserc} (Example 2). }
  \label{fig:ximoff}
\end{figure}

\begin{figure}[t]
  \centering
  \includegraphics[width=0.45\textwidth,trim=0 10 0 5,clip=True]{./figures/plot_sersic_test_Xh.pdf}
  \caption{Example for the ratio $C_{\widetilde{S}}/C_\mathcal{I}$ for the
    search with a 3D Gaussian template (see caption to Fig.~\ref{fig:ximoff} or
    Sect.~\ref{sec:verif-analyt-expr} for details) and sources described by a
    PSF core + S\`{e}rsic halo profile (Eq.~\ref{eq:42}), where we used the PSF
    Moffat model according to the parameters in Table~\ref{tab:mof}.  The plot
    shows the evaluation of Eq.~\eqref{eq:39} at $\lambda = 8700$\,\AA{} with
    the profile from Eq.~\eqref{eq:42} for the three average halo models of
    \citeauthor{Wisotzki2018} (\citeyear{Wisotzki2018}; their Extended Data Table 1) as indicated by the
    legend. }
  \label{fig:xhserc}
\end{figure}

The idealised situation analysed in the previous section, where astronomical
sources match exactly the source template used for the filtering, is never
encountered in the reality of blind surveys.  Each real world emission line
source will have different spectral and spatial profiles.  A motivated strategy
for the line search is to chose a template that maximises the SN for a large
fraction of the desired sources.  This may be achieved iteratively by varying
the template parameters and checking the SN response of sources found in a
previous iteration or even a previous survey in the same field (a possible
diagnostic diagram in this respect is shown in Fig.~7 of
\citetalias{Herenz2017}).  Potentially, one might even consider to construct a
catalogue from multiple search templates that are optimised for different source
characteristics.  In this case, additional bookkeeping is required to track which
template resulted in a source entering the catalogue, and multiple entries for
an individual sources recovered with different templates need to be weeded out.
We will not consider this case here.  The aim of this section is to provide
ideas towards the answer the question: ``\textit{What is the realistic selection
  function for a given source population in the catalogue that results from the
  search with a template $\bm{S}$?}''.  We understand as ``source population''
here a class of emission line objects, e.g., high-redshift LAEs, whose spatial
and spectral properties are understood in a statistical sense.

In order to obtain a realistic selection function we first have to calculate the
response of the 3D filtering operation in Eq.~\eqref{eq:27} to a source that is
not matched by the source template in Eq.~\eqref{eq:17}.  This can be done for
any emission line source profile $\bm{\widetilde{S}}$ by directly
evaluating the expectation value of Eq.~\eqref{eq:27} for the particular profile
and given line flux $F_\mathrm{line}$ and then Eq.~\eqref{eq:28} provides
$C_{\widetilde{S}}$ needed in Eq.~\eqref{eq:29}.  Here the subscript
$\widetilde{S}$ is meant to signify that this is the proportionality factor for
an non-template matching source profile.

To guide such calculations we now derive some useful analytic expressions for
$C_{\widetilde{S}}(\lambda)$ under the assumption that the source's spectral-
and spatial profile can be separated like our filtering template
(Eq.~\ref{eq:17} and Eq.~\ref{eq:18}).  In the spirit of Eq.~\eqref{eq:31}, we
place this source at position $x',y',z'$ into the data cuboid,
\begin{equation}
  \label{eq:35}
  F_{x,y,z} = A \cdot \widetilde{S}_{x-x',y-y'}^\mathrm{spat}
  \widetilde{S}_{z-z'}^\mathrm{spec} + n_{x,y,z} \; \text{,}
\end{equation}
where
\begin{equation}
  \label{eq:36}
  \widetilde{S}_{x,y,z} = \widetilde{S}_{xy}^\mathrm{spat}
\widetilde{S}_{z}^\mathrm{spec}
\end{equation}
is the normalised 3D profile
($\sum_{ijk} \widetilde{S}_{ijk} \equiv \sum_{ij}
\widetilde{S}_{ij}^\mathrm{spat} \equiv \sum_k \widetilde{S}_{k}^\mathrm{spec}
\equiv 1$) of the non-template matching source
($\widetilde{S}_{ij}^\mathrm{spat} \neq {S}_{ij}^\mathrm{spat}$ and
$\widetilde{S}_{k}^\mathrm{spec} \neq {S}_{k}^\mathrm{spec}$ for
${S}_{ij}^\mathrm{spat}$ and ${S}_{k}^\mathrm{spec}$ from Eq.~\ref{eq:17}) and
with the amplitude $A$ as defined in Eq.~\eqref{eq:33}.  The expectation value of
the filter in Eq.~\eqref{eq:27} at position $x',y',z'$ is then
\begin{equation}
  \label{eq:37}
  E[\mathcal{T}^{(x',y',z')}] =
  \frac{A^2}{E_\mathcal{I}[\mathcal{T}^{(x',y',z')}]} \cdot \sum_k
  \frac{S_k^{\{z'\}} \widetilde{S}_{k}^\mathrm{spec}}{\sigma^2_{z'-k}} \cdot
  \sum_{ij} S_{ij}^{\{z'\}} \widetilde{S}_{ij}^\mathrm{spat} \;\text{,}
\end{equation}
with $E_\mathcal{I} [\mathcal{T}^{(x',y',z')}]$ according to Eq.~\eqref{eq:32}.  Equating Eq.~\eqref{eq:37}
with Eq.~\eqref{eq:28} and using Eq.~\eqref{eq:34} we find
\begin{equation}
  \label{eq:38}
  C_{\widetilde{S}}(\lambda[z]) = \frac{1}{\Delta \lambda^2
    C_\mathcal{I}(\lambda[z])}
  \sum_k
  \frac{S_k^{\{z\}} \widetilde{S}_{k}^\mathrm{spec}}{\sigma^2_{z-k}} \cdot
  \sum_{ij} S_{ij}^{\{z\}} \widetilde{S}_{ij}^\mathrm{spat} \;\text{.}
\end{equation}
This equation can be evaluated for any assumed source profile described by
Eq.~\eqref{eq:36} to obtain the selection function with Eq.~\eqref{eq:29}.

Considering the case where only the spatial part of the source profile differs
from that of the template and that the spectral part is a perfect match
($\widetilde{S}_{ij}^\mathrm{spat} \neq {S}_{ij}^\mathrm{spat}$, but
$\widetilde{S}_{k}^\mathrm{spec} \equiv {S}_{k}^\mathrm{spec}$) we find
\begin{equation}
  \label{eq:39}
  C_{\widetilde{S}}(\lambda[z]) = C_\mathcal{I}(\lambda[z]) \cdot
  \frac{\sum_{ij} S_{ij}^{\{z\}} \widetilde{S}_{ij}^\mathrm{spat}}{\sum_{ij}
    \left ( S_{ij}^{\{z\}} \right )^2} \; \text{,}
\end{equation}
with $C_\mathcal{I}$ from Eq.~\eqref{eq:34}.

The particular simple expression for $C_{\widetilde{S}}(\lambda)$ in
Eq.~\eqref{eq:39} arises because we assumed the noise variance to be spatially
invariant (Eq.~\ref{eq:23}).  The problem of filter mismatch under this
assumption has been analysed by \citeauthor{Zackay2017} (\citeyear{Zackay2017};
their Appendix D) for 2D Gaussian profiles (Eq.~\ref{eq:43}).  For this case,
and under conditions that allow to replace the summations in Eq.~\eqref{eq:39}
by integrals, we have
\begin{equation}
  \label{eq:40}
  C_{\widetilde{S}}/C_\mathcal{I} = 2\kappa / (\kappa^2 + 1) \;\text{,}
\end{equation}
with $\kappa$ denoting the ratio between source dispersion,
$\sigma_{\widetilde{S}}$, and the template dispersion, $\sigma_G$, i.e., 
$\kappa = \sigma_{\widetilde{S}}/\sigma_G$ (see Fig.~1 in \citealt{Zackay2017}).
For more general mismatches between spatial templates and source profiles
Eq.~\eqref{eq:39} has to be evaluated numerically.  We now illustrate such
numerical evaluations of Eq.~\eqref{eq:39} in two examples.  Both examples are
calculated for the line search in the MXDF described in
Sect.~\ref{sec:verif-analyt-expr}.

\textbf{Example 1:} \textit{Moffat point sources.}  The curve shown in
Figure~\ref{fig:ximoff} demonstrates the ratio $C_{\widetilde{S}}/C_\mathcal{I}$
for perfect point sources whose spectral profile is given by a 1D Gaussian of
250\,km\,s$^{-1}$ width (Eq.~\ref{eq:46} with Eq.~\ref{lsd_eq:11}.  The point
sources were modelled as Moffat functions (Eq.~\ref{eq:44}), which provide an
adequate description of adaptive-optics assisted MUSE observations
\citep{Fusco2020}.  In particular we use the Moffat parameters provided by the
\texttt{muse-psfr} tool by \cite{Fusco2020}.  This software uses the telemetry
data of the ground layer adaptive optics system to derive the Moffat parameters
that are provided in Table~\ref{tab:mof} in \ref{sec:templ-param-used} for
reference.  It can be seen in our example the ratio varies slowly as a function
of wavelength and that the curve follows the behaviour of the FWHM difference
between assumed perfect point sources and used recovery template.  Recalling the
definition of $C(\lambda)$ in Eq.~\eqref{eq:28} we can read from this plot that
a point-source reaches the same SN expectation value at about 116\% the flux of
the ideal Gaussian source.  Of course, this is not the optimal SN for this
source, as this would be reached only if the Moffat profile had been used as a
search template.  However, this would lower the detectability of more extended
sources.  We also show in Figure~\ref{fig:ximoff} the case, where we use the
FWHM parameterisation of \texttt{muse-psfr} for Gaussian point sources (dashed
line), i.e., the mismatch described by Eq.~\eqref{eq:40}.  These Gaussian
sources reach the same expectation value in S/N at about 114\% the flux of the
perfect template matching sources.

\textbf{Example 2:} \textit{Moffat core + S\'{e}rsic halo profiles.}  The
spatial distribution of Ly$\alpha$ emission from high-redshift Ly$\alpha$
emitting galaxies can almost always be characterised by a spatially unresolved
``point-like'' core component and an extended halo component.
\citep[e.g.,][]{Steidel2011,Wisotzki2015,Leclercq2017,LujanNiemeyer2022}.
\cite{Wisotzki2018} parameterised the average profile of these sources with a
\cite{Sersic1968} profile\footnote{An English translation of the relevant
  chapter in the seminal \cite{Sersic1968} atlas that was published in Spanish
  can be found at \url{https://doi.org/10.5281/zenodo.2562394}.} for the halo
component:
\begin{equation}
  \label{eq:41}
  I_H(r) = I_H(0) \times \exp \left ( - \left [ \frac{r}{R_e} \right
    ]^{1/n} \right ) \; \text{.}
\end{equation}
The parameters $R_e$ and $n$ in Eq.~\eqref{eq:41} define the characteristic
radius and the kurtosis of the profile \citep{Graham2005}.  The observable
profile of the ``point-like + halo'' model can then be written as
\begin{equation}
  \label{eq:42}
  \widetilde{S}(r,\lambda) = I_0 \,\times\, \mathrm{PSF}(r,\lambda) \star \left \{ (1 - X_\mathrm{H} ) \delta(r)
    + X_\mathrm{H} I_H(r) \right \} \; \text{,}
\end{equation}
where $\delta(r)$ is the 2D delta-function, $\star$ denotes the 2D convolution,
$\mathrm{PSF}(r)$ denotes the model of the point spread function\footnote{The
  effects of PSF convolution on the profile in Eq.~\eqref{eq:41} have been
  analysed by \cite{Trujillo2001a} and \cite{Trujillo2001} for a Gaussian and a
  Moffat PSF model, respectively.}, $I_0$ is a constant chosen such that the
integral
$\int_0^{2\pi} \int_0^{\infty} S(r,\lambda) \, r^2 \,\mathrm{d}r\,\mathrm{d}\phi
\equiv 1$, and $X_\mathrm{H}$ ($\in [0,1]$) denotes the halo-flux fraction.  An
example for the required rescaling of this type of profile we show in
Figure~\ref{fig:xhserc} the ratio $C_{\widetilde{S}}/C_\mathcal{I}$ of
Eq.~\eqref{eq:39} calculated for the MXDF emission line search with a 3D
Gaussian template (Sect.~\ref{sec:verif-analyt-expr} and Table~\ref{tab:sn2}).
In particular, we here evaluated Eq.~\eqref{eq:39} at $\lambda = 8700$\,\AA{}
for the profile given in Eq.~\eqref{eq:42} with the Moffat PSF (Eq.~\ref{eq:44})
using the parameterisation from \texttt{muse-psfr} in Table~\ref{tab:mof} and
for three different S\`{e}rsic profiles (Eq.~\ref{eq:41}) with
$R_e = (0.86\arcsec{}, 0.9\arcsec{}, 1.67\arcsec{})$ and $n=(2.8, 3.3, 6.5)$.
These are the average best-fit parameters of the analytic profiles according to
Eqs.~\eqref{eq:41} and \eqref{eq:42} for Lyman $\alpha$ emitting galaxies
analysed by \citeauthor{Wisotzki2018} (\citeyear{Wisotzki2018}; their Extended
Data Table 2).  We ensured that the numerical evaluations of Eq.~\eqref{eq:41}
and Eq.~\eqref{eq:42} were not affected by sampling effects\footnote{The
  resulting profiles from those equations do not represent physical reality if
  the PSF and the S\'{e}rsic profile would only be evaluated at the centre of
  each spatial pixel, especially when $R_e$ is of similar dimensions as the
  spaxels and when $n$ becomes large \citep[see also][]{Peng2002}.  We note,
  that this is the case with the 2D S\'{e}rsic model provided in the Astropy
  package \citep{Astropy-Collaboration2018} in the most recent version 5.1. The
  model \texttt{astropy.modeling.functional\_models.Sersic2D} must therefore not
  be used to model observed profiles.  Here we performed the calculations on a
  grid where 625 spatial pixels correspond to one MUSE spaxel, as we found that
  the results did not change significantly when using finer grids.  We also
  assumed a zero pixel-phase, i.e.,  the peak is located at the centre of the
  central pixel.  The so obtained 2D profiles of Eq.~\eqref{eq:42} were then
  downsampled to the native MUSE resolution, and those resampled 2D profiles
  were then used for the evaluation of Eq.~\eqref{eq:39}.}.
Figure~\ref{fig:xhserc} shows, how the ratio increases for increasing halo flux
fractions.  For the typical halo flux fraction of 80\%-90\% found in
\cite{Wisotzki2018} the MXDF emission line search template introduced in
Sect.~\ref{sec:verif-analyt-expr} recovers sources with the same SN at around
70\% flux of a perfect template matching sources.  We note that this is not the
optimum, as in principle the profile Eq.~\eqref{eq:42} could be used as a
spatial template.  According to the matched filter theorem \citep[e.g.][see
also Sect.~\ref{sec:mf}]{Das1991} this would increase the depth even more for
such extended sources.  However, the use of such a template would negatively
impact the detectability of more compact emission line sources.

Above considerations provide us with $C_{\widetilde{S}}(\lambda)$ for individual
sources.  However, the numerical evaluation of the relevant expressions is
nearly instant, even when considering the required sub-sampling and re-binning
schemes needed when dealing with S\'{e}rsic type profiles.  Thus it appears
feasible to consider large ensembles of source templates $\widetilde{S}_i$ from
some population and then to calculate $C_{\widetilde{S}_i}$ for each of those
sources.  Next, Eq.~\eqref{eq:29} needs to be evaluated for each
$C_{\widetilde{S}_i}$ to obtain $f_{C_{\widetilde{S}_i}}$, and the final
selection function is then given by $ f = \sum_i w_i f_{C_{\widetilde{S}_i}}$,
where the weights ($\sum_i w_i \equiv 1$) have to be chosen according to the
expected occurrence rate of $\widetilde{S}_i$ within the population.
Alternatively, the source parameters or profiles are directly drawn from the
underlying population, and in this case
$f = N^{-1} \sum_i f_{C_{\widetilde{S}_i}}$, where $N$ is the number of draws.
The resulting selection function then provides a robust estimate of the
realistic selection function for the catalogue entries of such a population in a
IFS data cuboid.


\section{Summary and Concluding Remarks}  
\label{sec:sc}

This article presented a discrete matched filtering approach in three dimensions
that correctly accounts for rapidly varying variances along one dimension while
in the two other dimensions the variances are fixed.  The motivation for this
approach is the search for faint astronomical emission line signals in wide
field integral field spectroscopic datasets.  An implementation of the method is
provided in an updated version in the open source Python software \mbox{LSDCat}
that can be obtained from the Astrophysics Source Code Library:
\url{https://ascl.net/1612.002}.

As with the original \mbox{LSDCat}, the matched filtering routines are
implemented in \texttt{lsd\_cc\_spatial.py} and \texttt{lsd\_cc\_spectral.py}.
For \texttt{lsd\_cc\_spatial.py}, which computes the inner sum of
Eq.~\eqref{eq:27}, only a flux data cuboid and the template parameters
(\ref{sec:temp}) are needed as input.  The outer sum of Eq.~\eqref{eq:27} is
implemented in \texttt{lsd\_cc\_spectral.py}, and here the resulting temporary
cuboid from the inner sum and the effective variance spectrum are required as an
input.  If desired, the user can still use the classic algorithm
(Eq.~\ref{eq:19}) in order to reproduce results obtained with the original
implementation of LSDCat.

We demonstrated in this article, by making use of the publicly available
MUSE-Wide Data Release 1, that the updated algorithm provides indeed better SN
for emission line signals in the spectral vicinity of telluric air-glow lines.
We then discussed a useful property of the emission line search with a
matched-filter, namely that in such a search the selection function is
deterministic for a given variance spectrum provided that the noise is behaving
according to the expectations.  Under this provision the selection function can
be expressed by a simple analytic formula (Eq.~\ref{eq:29}), with a factor of
proportionality, $C(\lambda)$, that only depends on the congruence between
assumed source profile and actual source profiles.  We provide an expression for
$C(\lambda)$ for the idealised case where template and source match exactly
(Eq.~\ref{eq:34}).

We then verified the analytic selection function against a source insertion
experiment.  This experiment was carried out in the data cuboid from the deepest
MUSE survey ever obtained, the $t_\mathrm{exp} = 141$\,h MUSE eXtreme Deep
Field.  In this dataset the analytic expression for an idealised selection
function, where the sources match the template exactly, shows excellent
congruence with a selection function derived from a source insertion and
recovery experiment.  We remark that the MXDF dataset is somewhat special in
the sense that the observing strategy was chosen to homogenise the background
noise and to remove any potential residual systematics.  Nevertheless, the
assumption of spatially invariant noise is also met to a sufficient degree in
shallower MUSE data, and it was shown for the standard dithering strategy of the
deep MUSE observations of the MUSE Hubble Deep Field South survey
\citep{Bacon2015} that the noise scales only with some moderate deviation from
the expected $1/\sqrt{N_\mathrm{exp}}$ behaviour.  Nevertheless, a full verification
of the here presented approach with shallower data cuboids is desirable in the
future.

We also presented ideas regarding how the deterministic selection function for
the idealised case can be used to obtain realistic selection functions for line
source profiles that are not fully congruent with the template.  In particular,
we provided a formula that allows to calculate the effect on the factor
$C(\lambda)$ for source profiles that differ from the template profiles 
(Eq.~\ref{eq:38}).  We then analysed two example situations of spatial profile
mismatches.  We argued that this method can be extended to obtain realistic
selection functions for catalogue entries of a population of emission line
sources, where the spatial and spectral properties have been understood in a
statistical sense.

Using an algorithm with a deterministic selection function removes the need for
computationally cumbersome source insertion and recovery experiments.  Moreover,
such an algorithm may also be used to efficiently design IFS survey
observations.  State-of-the art telescope facilities require a reliable estimate
of the required exposure time already at the application stage.  To this aim
interactive exposure time calculators are provided.  The calculations of those
tools use a model of the detectors and the atmosphere to provide a reliable
estimate of the background noise.  Such estimates of the background noise could
then be used with the here presented approach to predict, e.g., catalogue
incidences of particular emission line galaxies given their line luminosity
function.

Some important aspects of the line source detection problem in integral field
spectroscopic data were not discussed here, but should be mentioned briefly.
First, we did not make direct statements regarding the reliability (sometimes
also dubbed purity) on the samples from our catalogue.  The reliability $R$ is
defined as the complement of the false detection probability $p_\mathrm{F}$:
$R = 1 - p_\mathrm{F}$ \citep[e.g.][]{Hong2014}.  As noted in
Sect.~\ref{sec:mf}, the standard conversion between $p_\mathrm{F}$ and detection
threshold is formally not correct.  A series of articles addressed this issue in
2D \citep{Vio2016,Vio2017,Vio2019}.  While the details are mathematically
involved, the upshot is that the detection threshold has to be chosen more
conservatively than what the standard expectation based on a Gaussian
distribution would provide.  An often used empirical approach to quantify the
reliability to use negated datasets, which works well under the assumption of
the noise being symmetric \citep[see also][]{Serra2012}.  Second, the
construction of the selection function is based on the ground truth for
$F_\mathrm{line}$, which however is unknown for observed sources, where we
measure $F_\mathrm{line}^\mathrm{obs} \pm \Delta F_\mathrm{line}^\mathrm{obs}$
to characterise the distribution of possible line fluxes.  This needs to be
taken into account, when we use the selection function for modelling purposes,
and the well known Eddington-Malmquist bias in luminosity function
determinations results from this effect \citep[e.g., Chapter 5.5
in][]{Ivezic2014}.  One way to avoid such biases are additional cuts of the
sample by only using sources where the errors on the flux measurements do not
correspond to significant changes in the selection function, i.e., where
$f_C(F_\mathrm{line}^\mathrm{obs} \pm \Delta F_\mathrm{line}^\mathrm{obs},
\lambda) \simeq f_C(F_\mathrm{line}^\mathrm{obs}, \lambda)$.  However, this may
drastically reduce the sample, and an alternative is to account for the
measurement errors in the modelling process \citep[][]{Rix2021}.  Last, a
current bottleneck in the construction of line emitter catalogues based on
wide-field IFS data is the manual classification step.  Related to this are
spurious line signals due to imperfect continuum subtraction of sources with
strong continuum emission.  These residuals also have to be weeded out manually
from catalogues.  We advocate future research on those latter issues by
investigating machine learning techniques for classification and better
continuum subtraction methods that do not leave strong residuals.  While matched
filtering is a useful method to uncover the faintest emission line galaxies in
wide-field IFS datasets, ultimately it is only the first step in the detection
and analysis chain, and subsequent steps are required before astrophysical facts
can be inferred from the data.  Arguably, however, a good understanding of this
first step in the chain is required for successful statistical analyses of IFS
survey data.


\section*{Acknowledgments}

I wish to express my thankfulness to Lutz Wisotzki for sparking off the initial
idea for the work presented here.  I thank Leindert Boogaard from the
Max-Planck-Institut f\"{u}r Astronomie and Friedrich Anders from the Departament
de F\'{i}sica Qu\`{a}ntica i Astrof\'{i}sica at the University of Barcelona for
careful proofreading the manuscript prior to submission.  My special thanks go
to the staff working at ESO/Chile and the Paranal Observatory, especially
Fernando Selman and Fuyan Bian, for a interesting and educational time.  I wish
also to express my gratitude to Pascale Hibon for organising the ``Joint
Observatories Kavli Science Forum in Chile'', where I could discuss the ideas
expressed in this article publicly for the first time
\citep{2022joks.confE..12H}.  Lastly, I thank Roland Bacon for envisioning and
leading the MUSE project.

\normalsize



