# Data Combination: Interferometry and Single-Dish Imaging in Radio Astronomy

Adele Plunkett National Radio Astronomy Observatory, 520 Edgemont Road,
Charlottesville, VA 22903, USA  Alvaro Hacar Department of Astrophysics,
University of Vienna, Türkenschanzstrasse 17, 1180, Vienna (Austria)  Leiden
Observatory, Leiden University, PO Box 9513, 2300-RA Leiden (The Netherlands)
Lydia Moser-Fischer  Argelander-Institut fuer Astronomie, Auf dem Huegel 71,
D-53121 Bonn, Germany  Dirk Petry European Southern Observatory, Karl-
Schwarzschild-Strasse 2, 85748 Garching, Germany  Peter Teuben Astronomy
Department, University of Maryland, 4296 Stadium Dr., College Park, MD 20742,
USA  Nickolas Pingel Department of Astronomy, University of Wisconsin–Madison,
475 N. Charter Street, Madison, WI 53706, USA  Research School of Astronomy
and Astrophysics, The Australian National University, Canberra, ACT 2611,
Australia  Devaky Kunneriath National Radio Astronomy Observatory, 520
Edgemont Road, Charlottesville, VA 22903, USA  Toshinobu Takagi  Japan Space
Forum, 3-2-1 Kandasurugadai, Chiyoda-ku, Tokyo 101-0062, Japan  Yusuke
Miyamoto  Department of Electrical and Electronic Engineering, Fukui
University of Technology, 3-6-1, Gakuen, Fukui, 910-8505, Japan  National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Emily Moravec Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA
Sümeyye Suri Department of Astrophysics, University of Vienna,
Türkenschanzstrasse 17, 1180, Vienna (Austria)  Kelley M. Hess Instituto de
Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía s/n, 18008 Granada,
Spain  ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990
AA, Dwingeloo, The Netherlands  Melissa Hoffman Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD, 21218, USA  Brian Mason
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA
22903, USA

(Received 2022 December; Revised 2023 January; Accepted 2023 February)

###### Abstract

Modern interferometers routinely provide radio-astronomical images down to
sub-arcsecond resolution. However, interferometers filter out spatial scales
larger than those sampled by the shortest baselines, which affects the
measurement of both spatial and spectral features. Complementary single dish
data is vital for recovering the true flux distribution of spatially resolved
astronomical sources with such extended emission. In this work, we provide an
overview of the prominent available methods to combine single dish and
interferometric observations. We test each of these methods in the framework
of the CASA data analysis software package on both synthetic continuum and
observed spectral datasets. We develop a set of new assessment tools which are
generally applicable to all radio-astronomical cases of data combination.
Applying these new assessment diagnostics we evaluate the methods’ performance
and demonstrate the significant improvement of the combined results in
comparison to purely interferometric reductions. We provide combination and
assessment scripts as add-on material. Our results highlight the advantage to
use data combination to ensure high-quality science images of spatially
resolved objects.

Millimeter astronomy(1061) — Submillimeter astronomy(1647) —
Interferometry(808) — Fast Fourier transform(1958) — Spectroscopy(1558) —
Astronomical technique(1684)

††journal: PASP††facilities: ALMA††software: CASA (The CASA Team et al.,
2022), astropy (Astropy Collaboration et al., 2013, 2018), numpy (Harris et
al., 2020), scipy (Virtanen et al., 2020), matplotlib (Hunter, 2007),
analysisUtils (https://casaguides.nrao.edu/index.php/Analysis_Utilities)

##  1 Introduction

A well-known challenge of interferometric imaging in the radio and sub-
millimeter regimes is that it relies on aperture synthesis in which the
Fourier or _(u,v)_ plane is irregularly or incompletely sampled. While
interferometry is a powerful technique for resolving structures on smaller
scales than can be attained with a single dish telescope, on the other hand, a
single dish telescope is necessary for recovering extended structure. This is
commonly known as the “short-spacing” problem that we show in Figure 1,
visualized as a “hole” in the inner region of the _(u,v)_ coverage (e.g.
Wilner & Welch, 1994; Kurono et al., 2009); visualizations in the literature
are given by Mason (2020) (see their §2.2 and Figures 2-4) and Leroy et al.
(2021) (Figures 24-25), among others.  The lack of short spacings in the
interferometer results in the inability to detect the flux of large spatial
scales and thus leads to an incomplete representation of the true sky
brightness distribution in the image. Ideally, observations would sample the
complete _(u,v)_ space corresponding to the spatial scales under study, and
recover consistent amplitudes at overlapping _(u,v)_ distances if multiple
interferometry configurations and single dish observations are required (see
e.g., Figure 3 from Koda et al., 2019). As described in Rau & Cornwell (2011),
the impact of the zero-spacing problem is not restricted to continuum images.
Spectra can also be strongly affected, as shown in molecular line observations
(e.g. Pety et al., 2013; Hacar et al., 2018, see also Sect.6.2).

Figure 1:  Central _(u,v)_ coverage of the M100 dataset presented in this work
showing all baselines of $\leq$ 50m observed by the 12m (blue), 7m (red), and
Total Power (gray) arrays. Different ellipses highlight the overlapping uv-
coverage between the 12m and 7m (hatched orange area) and between the 7m and
TP (hatched black area) arrays, respectively. Significant overlap between
different arrays is preferred for data combination.

The solution is to complement the interferometric observations by a separate
observation with an instrument sensitive to the large scales, typically a
single-dish telescope, using data combination. For example, in the case of the
Atacama Large Millimeter/sub-millimeter Array (ALMA), this is done by
combining data from the main ALMA-12m array and/or the ALMA-7m Compact Array
(ACA) with data from the dedicated ALMA single-dish (SD, or Total Power
denoted as TP) telescopes. A broad literature highlights the critical role of
data combination in both Galactic and extragalactic studies including a large
variety of interferometer plus SD combinations not only within ALMA (for a
very thorough example, see Leroy et al., 2021, among others) but also using
CARMA and NRO 45m/IRAM-30m/FCRAO (Koda et al., 2011; Kong et al., 2018;
Plunkett et al., 2013, 2015), SMA and APEX (Kauffmann et al., 2017), (E)JVLA
and GBT (Pineda et al., 2011; Williams et al., 2018; Wang et al., 2020),
PdBI/NOEMA and IRAM-30m (Pety et al., 2013; Beuther et al., 2022), ALMA-12m
and MOPRA/IRAM-30m (Peretto et al., 2013; Hacar et al., 2018), or ASKAP and
Parkes (Pingel et al., 2022) telescopes, to cite some examples.

Interferometric image synthesis and data combination is not a new topic.
Stanimirovic (2002) provides a comprehensive introduction to the problem and
the history of solution attempts. An update is given by Mason (2020)
concentrating on the “feathering” technique of data combination. Most
generally, methods exist for combination at different steps of the imaging
process, either in the spatial frequency (Fourier) domain or in the image
domain. However, a systematic and quantitative comparison of the results of
these methods is still missing.

In this work we summarize current data combination techniques and provide new
tools to evaluate their performance. After setting the foundation by briefly
describing important aspects of imaging purely interferometric data in §2, we
present the prominent existing methods for combining interferometric and
single-dish observations in §3. In §4, we specify a simulated and an observed
test dataset for our evaluation of methods. In §5, we describe in detail our
new assessment metrics for interferometric and combined images which are
implemented in our new performance evaluation tools. We then run the described
combination methods on the test datasets and evaluate the results. Finally, in
section §6 we derive general recommendations for the use of the combination
methods. The script suite for running the combination methods and our
evaluation tools, and the data to reproduce the figures in this paper are
available from https://github.com/teuben/DataComb.

##  2 Imaging Interferometric Data with CLEAN

There is a Fourier transform relationship between the signals from the
antennas making up an interferometer and the sky brightness distribution. The
correlated signals, known as visibilities, are visualised in the $(u,v)$
plane, and an inverse Fourier transform is used to provide an image of the
spatial plane (see Mason, 2020). The point spread function (PSF) of the
interferometer is the Fourier transform of the visibilities, and as these do
not fill the $(u,v)$ plane, the PSF can be highly irregular. An initial step
of imaging produces the so-called ‘dirty’ image, which is the convolution of
this PSF with the true sky brightness distribution. PSF sidelobes manifest as
image artifacts. Thus, subsequent deconvolution is necessary to optimise image
quality. Presently, the most commonly used deconvolution algorithms are based
on CLEAN (Högbom, 1974; Clark, 1980).

CLEAN deconvolves the data from the effects of the PSF of the interferomenter
by iteratively modeling the emission as a set of point sources (delta
functions). The final image and residuals are convolved with an idealised
restoring beam based on the FWHM of the PSF. CLEAN achieves a good image
quality for astronomical objects comparable in size to the angular resolution.
For objects significantly larger, typically multi-scale (MS) CLEAN (Cornwell,
2008) is used, which employs inverted truncated paraboloids instead of delta
functions to model the emission. CLEAN algorithms also incorporate options for
weighting visibilities and for interpolating into the missing spacings in the
$(u,v)$ plane, which provide a trade-off between resolution and sensitivity,
but can themselves introduce artifacts if the sky brightness has structure on
the unsampled scales such that the interpolation is inaccurate. Further
extentions of CLEAN (MT-MFS CLEAN, Rau & Cornwell, 2011) account for the
frequency dependence of the _(u,v)_ coverage.

For actual data analysis work in this study, we use exclusively the CASA
(Common Astronomy Software Applications) package (The CASA Team et al., 2022)
version 6.1. CASA is the primary data analysis package for ALMA and the Karl
G. Jansky Very Large Array (VLA). Other available radio-astronomical data
analysis packages include AIPS (van Moorsel et al., 1996), GILDAS (Pety, 2005;
Gildas Team, 2013), MIRIAD (Sault et al., 1995) and more recently developed
programs such as wsclean (Offringa et al., 2014), resolve (Arras et al.,
2021), and purify (Carrillo et al., 2014).

In CASA, CLEAN and its extensions are all accessible via the tclean task.
Additionally, all the data combination methods presented in §3 can be
implemented by calling CASA tasks. Some are already available as a single CASA
task (such as sdintimaging and feather) or add-on packages for CASA (tp2vis).
Each method in §3 uses tclean (or the underlying libraries of tclean) at some
step in the imaging and data combination process.

The user can control the behaviour of tclean via a large set of input
parameters which are explained in the CASA
documentation111https://casadocs.readthedocs.io/. Some of the most critical
parameters are those related to (a) the number of deconvolution iterations,
(b) the CLEAN threshold (i.e. the minimum flux justifying the placement of a
CLEAN component), (c) the spatial masking (where the CLEAN components can be
placed), and (d) the handling of the relative weighting of the interferometric
visibilities.

To facilitate the comparison among the combination methods, we fix all tclean
parameters to common values. We set the number of iterations to 100,000. We
use Briggs weighting (Briggs, 1995) with robust value of 0.5 (i.e. balance
between natural and uniform weighting) for all runs of tclean as well as
gridder = ‘mosaic’; in this study we deal exclusively with mosaic maps, but
gridder = ‘mosaic’ applies for the case of combining data from multiple arrays
of different antenna sizes, even if they are single pointing observations.

We explored several masking options, as careful masking is important to
mitigate constructive interference of the sidelobes of the synthesized beam
(Condon et al., 1998). The mask is generated by first utilizing the SD image
to mask (include) all regions above 0.3 times the max value in the SD image,
and in the next step include any additional regions selected by the “auto-
multithresh” algorithm in tclean (Kepley et al., 2020), using the standard
“auto-multithresh” parameters. This masking technique is determined to capture
the regions of larger-scale emission (as in the SD image), as well as any
regions of more compact emission that might appear primarily in the
interferometry data. Finally, the spectral definition mode can be chosen to be
either ‘mfs’ (multi-frequency synthesis, for continuum imaging with only one
output image channel) or ‘cube’ (for spectral line imaging with one or more
channels), depending on the dataset and desired image. In this manuscript, we
present continuum and spectral data, as an image and cube, respectively.

##  3 Combination Methods

In this section, we describe the theory and execution of the currently
available data combination methods which we have selected for evaluation in
this study (see also Table 1). The methods incorporate one or more of the
following steps in order to combine the interferometry dataset(s) with SD
observations: combine the SD data with a clean interferometer image; combine
the SD data with the dirty interferometer image before joint deconvolution;
convert the single-dish data to pseudo-visibilities, combine in the ($u,v$)
plane and then FT the combined data and CLEAN. The methods we present for
evaluation are: Feather (§3.1), SDINT (§3.2), a model-assisted CLEAN plus
Feather (hereafter MACF) (§3.3), TP2VIS (§3.4), and Faridani et al. (2018)’s
Short Spacing Correction (hereafter FSSC) (§3.5). For each of these methods we
include a list of the control parameters.

We note that the prerequisites common to all combination methods presented
here are the following:

  * •

all data have been fully calibrated (this would also include the possibility
of self-calibration if applicable);

  * •

astrometry should align;

  * •

images have overlapping, well-sampled spatial frequencies, i.e. adequate
coverage of the _(u,v)_ -plane among combined datasets;

  * •

images have well-defined beams: primary beam (PB) of low-resolution image,
$\text{PB}_{low-res}$ and primary beam of the high-resolution image,
$\text{PB}_{hi-res}$; and

  * •

bandwidth is narrow enough that the frequency dependent scaling of the
baseline lengths and the primary beam can be ignored, although SDINT can in
fact take this into account.

We assume that flux scales are aligned, but further fine-tuning can be done
for specific use cases. Generally, every input image is assumed to be
corrected for all instrumental effects including PB-response. We refer to the
documentation of each method (i.e. CASA Guides, Github, etc.) for guidance on
handling units and the specific treatments of the PB-response during the
respective steps of the data combination (sometimes referred to as DC)
process. Finally, we use test cases without bright emission that extends or
emerges beyond the area of the well-characterised PB-response. The same
methods should, in principle, be adaptable to other scenarios.

Table 1: Summary of data combination methods

Methodology | Domain | Method | Task name | Input | Output
---|---|---|---|---|---
[1] | [2] |  |  | interferometry | single dish |
| F/I | SDINT | sdintimaging | vis | image | image
Before |  |  |  |  |  |
| F | TP2VIS |  tp2vis[3] | vis | SD image | pseudo-visibilities
|  |  | tclean | vis | pseudo-visibilities | image
During | F+I | MACF | tclean | vis | image as model | image
|  |  | feather | image | image | image
| F | Feather | feather | image | image | image
After |  |  |  |  |  |
| I | FSSC | (script) | image | image | image

Note. — [1] indicates combination before, during, or after image
deconvolution; [2] Fourier (“F”) or image (“I”) domain in which the method
operates; [3] only available in CASA after import of the tp2vis package

###  3.1 Feather

####  3.1.1 Feather algorithm

The Feather algorithm is probably the most well-known method for combining SD
and interferometeric data. The algorithm is implemented as a CASA task of the
same name. This method inverts the input images via Fourier transform (FT) and
combines them in Fourier space (i.e. the _(u,v)_ -plane), weighting them
according to the spatial frequencies of the response of each telescope. The
combined data is then Fourier transformed back to the image plane.

Although feather does regrid, in some cases, regridding may be better to do
outside the task, for which case imregrid and specsmooth can be used. For a
more comprehensive discussion of feather, see Cotton (2017); the CASA
documentation fully explains additional parameters.222An alternative
implementation of the Feather algorithm is possible with the uvcombine package
in Python, which implements a similar approach to CASA’s feather task, but
with additional options. See https://github.com/radio-astro-tools/uvcombine.
immerge in Miriad offers the option feather, allowing the combination of two
images with Gaussian PSFs (see https://github.com/baobabyoo/almica for
implementation). Additionally, J-comb (Jiao et al., 2022) is a new linear
combination technique that has been benchmarked against both CASA’s feather
and Miriad’s immerge tasks, but is not further explored in this work. See Jiao
et al. (2022) for additional information about their investigation into the
feather and immerge implementations.

In more detail, the step-by-step order of operations undertaken with the
feather task can be enumerated by the following:

  1. 1.

Regrids the SD/low-resolution image spatially and spectrally to the high-
resolution image.

  2. 2.

FT both images, such that the next steps operate on the visibilities.

  3. 3.

Scales the FT low-resolution image by the ratio of clean beams (CBs):

$CB_{hi-res}/CB_{low-res},$ (1)

in order to account for different beam-sizes.

  4. 4.

Adds the FT of the high-resolution image $\times(1-\omega t)$ to the scaled
FT-ed low-resolution image, where $\omega t$ is the FT of $CB_{low-res}$. In
other words, the images have been weighted and combined in the $(u,v)$-plane.

  5. 5.

FT back to the imaging plane.

####  3.1.2 Feather Parameters

From the CASA documentation, the feather command inputs include those listed
in Table 2. Two particularly important parameters to consider for feather are
sdfactor and effdishdiam. sdfactor is used to adjust the flux scale of the SD
image. It should be constrained by comparing flux where there is overlap
between the spatial-frequencies in the low and high resolution images, and
empirically, the value is most commonly set close to 1.0-1.2. The parameter
sdfactor has parallels to sdgain in sdintimaging in SDINT (see § 3.2), but
they are not identical mathematically. effdishdiam is the effective diameter
of the SD telescope, corresponding to the “low-res” image in the feather
command; the effective dish diameter depends on the aperture efficiency, and
is generally determined/reported for a given telescope. When Feathering, the
weighting function for the SD data is usually the Fourier transform of the SD
primary beam (PBSD), but the weighting function can be altered by indicating a
reduced SD diameter, corresponding to a different PBSD. We did not further
investigate the impact of lowpassfiltersd.

Table 2: Parameters for method Feather CASA Task | Parameter | Default | Description
---|---|---|---
feather | imagename | = ‘’ | name of output feathered image
| highres | = ‘’ | name of high resolution (interferometer) image
| lowres | = ‘’ | name of low resolution (single dish) image
| sdfactor | = 1.0 | scale factor to apply to SD image
| effdishdiam | = -1.0 | effective diameter (in meters) of the SD telesope.
| lowpassfiltersd | = False | filter out the high spatial frequencies of the SD image

###  3.2 Joint Deconvolution: SDINT

####  3.2.1 The SDINT algorithm

The _SDINT_ algorithm (Rau et al., 2019) permits joint deconvolution of
wideband SD and interferometer data, and has been implemented as the task
sdintimaging in CASA333The task sdintimaging is incorporated in CASA 5.7/6.1
and subsequent releases.. We refer the reader to Rau et al. (2019) for full
details and explanatory diagrams, and we briefly describe the algorithm here.
The main inputs to the algorithm are the SD image and PSF (as a cube, if
applicable), as well as the interferometric dataset. Note that if no SD PSF is
provided, sdintimaging can solve for the beam information based on the image
header.

The sdintimaging algorithm follows these steps:

  1. 1.

The process starts as for standard “CLEANing” with a “major cycle” in which
the interferometric data are gridded and Fourier-transformed into an image
cube. At the same time, the corresponding PSF is formed.

  2. 2.

Then the SD cube is combined with the interferometric cube in a Feathering
step (see §3.1).

  3. 3.

A similar combination is applied to the interferometric PSF and the
representative SD kernel to form a joint PSF. The joint image and PSF cubes
then form inputs to a standard CLEAN “minor cycle” deconvolution.

  4. 4.

As the minor cycle limits are reached, model images from the deconvolution are
translated back (depending on the chosen deconvolution algorithm) to model
image cubes which are then (a) subtracted from the SD image cube to form a new
residual SD image and (b) Fourier-transformed and degridded to form a new set
of residual interferometric visibilities.

  5. 5.

These visibilities together with the residual SD image then form the input to
the next major cycle and Feathering. For mosaicking, appropriate application
of the primary beam correction takes place prior to deconvolution.

Because it is based on the standard CLEAN family of algorithms, the task
sdintimaging shares many parameters with the task tclean in CASA, as well as
controls for a Feathering step to combine interferometer and SD datasets
within the imaging iterations. The most important parameter to control the
Feathering step is sdgain. This parameter decides on the relative weight of
interferometric and SD data similar to the relative weight for the
visibilities from different baselines in pure interferometric imaging. The
default setting for sdgain is 1.0, i.e. equal weight to both contributions. A
value smaller than 1.0 will give the SD contribution less weight than the
interferometric one, while a value greater than 1.0 gives it more weight. The
flux scale of the joint image, however, will be kept constant. The best way to
understand the role of sdgain is to think of the combination as forming a
weighted mean of two measurements. In standard error propagation, the relative
weight of different measurements of the same quantity should be

$w_{i}=\frac{\sigma_{i}^{-2}}{\sum_{j}\sigma_{j}^{-2}}$ (2)

where $w_{i}$ is the weight and $\sigma_{i}$ is the error on the $i$th
measurement. The natural sdgain value is thus a function of the noise or
sensitivity of the SD and the interferometric contributions and should be
close to 1.0 if the sensitivities are similar, larger than 1.0 if the SD
observation is more sensitive and smaller than 1.0 if it is less sensitive
than the interferometric one.

####  3.2.2 SDINT Parameters

Table 3: Parameters of method SDINT CASA Task | Parameter | Default | Description
---|---|---|---
sdintimaging | usedata | ‘sdint’ | output image type (int, sd, sdint)
| vis | ‘’ | input interferometric visibility data
| sdimage | ‘’ | input single dish image
| sdpsf | ‘’ | input single dish PSF image
| sdgain | 1.0 | A factor or gain to adjust single dish flux scale
| dishdia | 100.0 | single dish diameter (in meters)
| selectdata | ‘’ | enable data selection parameters, as in tclean

Here we specify parameters specifically applicable when usedata=‘sdint’, which
is the case we explore for combination in this work. dishdia should be
considered as the effective dish diameter of the SD, and if sdpsf is given as
a dish diameter, this parameter will be ignored. For a summary of pertinent
parameters see Table 3.

###  3.3 Model-assisted CLEAN plus Feather method (MACF)

There is another method which is in essence a variant of feather, and we refer
to it as the Model-assisted CLEAN plus Feather method or MACF (for a list of
the input parameters see Table 4). It utilizes feather to ultimately combine
the interferometric image and the SD image (as in §3.1), but takes a slightly
different approach to generate the interferometric image used as input. The
parameter startmodel is used when invoking tclean on the interferometric data,
providing the SD image (in units of Janskys per pixel) as the starting model
image444A similar procedure to our MACF is documented by J. Kauffmann and is
available at http://tinyurl.com/zero-spacing .

The advantage of cleaning with the single dish image as a starting model is
that the extended emission is conveyed in the model image, and then
interferometer-based clean components are incorporated, before convolving with
a clean beam for the final image restoration (adding residuals). A caveat is
that the zero-spacing flux, although initialized, is not constrained by
tclean, so this should be considered an intermediate image only.

The next step is to combine this interferometry-with-startmodel image with the
SD image using feather. This step in the combination procedure is identical to
that described in §3.1. Examples of the application of the MACF method to
observations can be found in Dirienzo et al. (2015), Kauffmann et al. (2017)
and Hacar et al. (2018).

Table 4: Parameters for method MACF (model-assisted CLEAN plus Feather) CASA Task | Parameter | Default | Description
---|---|---|---
tclean | startmodel | = ‘’ | Name of low resolution (single dish) image
feather | imagename | = ‘’ | Name of output feathered image
| highres | = ‘’ | Name of high resolution (interferometer) image,
|  |  | generated after using low resolution image as ‘startmodel’
| lowres | = ‘’ | Name of low resolution (single dish) image

###  3.4 TP2VIS

The package Total Power Map to Visibilities (TP2VIS)555The TP2VIS package is
available via Github at: https://github.com/tp2vis/distribute was developed
for joint deconvolution of ALMA-12m, 7m, and TP data, and is thoroughly
presented by Koda et al. (2019). We reiterate that Total Power, or TP, is
another name used by radio astronomers for SD data. TP2VIS operates in CASA to
convert a SD map into visibilities; these so-called “pseudo-visibilities” can
then be used as input for joint deconvolution along with interferometric
visibilities. The technique has been presented and utilized by Rodríguez-
Fernández et al. (2008, & references therein); Kurono et al. (2009, &
references therein); Pety & Rodríguez-Fernández (2010, & references therein);
Koda et al. (2011, & references therein), with the TP2VIS package being the
first and only implementation into CASA.

TP2VIS is roughly structured in four steps, with much more detail on each
available in Koda et al. (2019):

  1. 1.

The SD map is converted into the sky brightness distribution, to be observed
as if by a virtual interferometer.

  2. 2.

The tp2vis function converts the brightness distribution into visibilities.

  3. 3.

Next, there is an option to set the weights of the TP visibilities according
to a few weighting schemes, such as that they represent the rms noise of the
original TP map or are set to a constant value.

  4. 4.

Finally, combination is done by deconvolving the pseudo-visibilities and the
interferomentric visibilities jointly using the tclean task.

The TP2VIS package also offers the function tp2vistweak, which attempts to fix
the problem of discrepant beamsizes of the dirty and clean beams in the image
space after deconvolution, as detailed in Jorsater & van Moorsel (1995).

In our tests with CASA 6,666We note that TP2VIS was developed with an earlier
version of CASA, as reported by Koda et al. (2019). TP2VIS turned out to be
less straight-forward to use effectively and consistently across different
datasets, compared with the other combination methods described here. In order
not to further delay the publication, we decided not to include it in the
evaluation section of this paper. We have, however, included it in the
accompanying software (§3.6).

###  3.5 Faridani’s Short Spacing Correction (FSSC) method

Implementations which utilize Fast Fourier Transforms (FFTs), such as the
Feather method, may introduce artifacts in the final combined image due to
structures that extend up to or beyond the edge of the image, creating a sharp
discontinuity. This is the case for the feather task. Faridani et al. (2018)
introduced a combination method that works purely in the image domain to avoid
the introduction of such artifacts; they provide a stand-alone Python code
that we have also incorporated into our suite of combination methods in order
for the implementation to be done consistently. The method is succinctly
summarized in Equation 6 of Faridani et al. (2018):

$\displaystyle I_{\rm missing}$ $\displaystyle=I^{\rm reg}_{\rm sd}-I^{\rm
conv}_{\rm int}$ (3) $\displaystyle I_{\rm comb}$ $\displaystyle=I_{\rm
int}+\alpha\cdot I_{\rm missing}.$

The missing flux ($I_{\rm missing}$) of the interferometry observations is
determined by first taking the difference between the SD image at its native
resolution $I^{\rm reg}_{\rm sd}$, and the interferometer-only image that has
been convolved and regridded to the angular resolution and pixel scale of the
SD image, $I^{\rm conv}_{\rm int}$. The combined image is then the sum of the
original interferometer-only image $I_{\rm int}$ and this missing flux image,
also scaled by $\alpha$ which is the ratio of the clean interferometer beam to
the SD beam to account for the difference in beam sizes (our Equation 1). In
addition to mitigating potential FFT artifacts, this approach is also less
computationally expensive since no FFTs are required. Similar to feather, an
additional scaling factor can be applied to the SD image to compensate for
differences in the flux scale beyond the differences in beam sizes.

We refer to this process as the Faridani’s Short Spacing Correction method
(hereafter FSSC) and note that the output images are very similar to those
from the Feather technique. While the concepts are equivalent, the FSSC method
operates purely in the image plane, with no Fourier transforms, contrary to
Feather. This technique has successfully been used to combine images from the
Large Millimeter Telescope, Planck, and the Caltech Submillimeter Observatory
to create a 1.1 mm image of the Central Molecular Zone of our Galaxy sensitive
to all angular scales down to 10.5″ (Tang et al., 2021a, b).

###  3.6 Software developed in this study

In this work, we provide a single data combination script (called “DC_script”,
that runs with Python and incorporates Astropy, Astropy Collaboration et al.
(2018)) that performs (a) the data combination methods described in §3 on
simulated and observational datasets (§4.1 and §4.2) and (b) assessment
methods described in §5. The script, data, and documentation are available at
https://github.com/teuben/DataComb. When running the script for the first
time, follow the steps in the README file. Also available in the code
repository are an overview of all steps, a ”quick start” guide, and template
parameter files. “DC_script” is optimized to run non-interactively, but the
user can choose to run interactively if needed in order to change certain
parameters. “DC_script” also includes the full set of quality assessment
metrics described in this work to evaluate the performance of the different
data combination methods (§ 5).

##  4 Test datasets

###  4.1 Simulated continuum observations

Only with simulated observations can the accuracy of the image reconstruction
be properly measured, since in that case the “true” intensity distribution on
the sky is known. For efficiency we try to capture a large range of possible
source structures in a single test image to study. An elegant way of achieving
this is to emulate a multi-scale molecular cloud with a power spectrum that
follows a power law. We follow the procedure described in Koda et al. (2019).
In addition we insert two bright point sources, one in a bright and one in a
faint region. We call the resulting model the “skymodel”. The model has
$4,096^{2}$ pixels with a size of 0.05″ each, resulting in an image size of
204.8″ each side, shown in Fig. 2.

Figure 2: The synthetic ”skymodel” image used as input for the simulated
observations.

With the “skymodel” as input, we ran the simobserve task in CASA version 5.6.1
to create simulated observational data in MeasurementSet format (MSs).
Following the ALMA Technical Handbook (Remijan et al., 2019), we adopted two
configurations for 12m array (C-4 and C-1 in cycle 7) along with ACA 7m and
SD. The expected beam size for this configuration is $\sim$1″ for an observing
frequency of 115 GHz, and we covered 2 GHz bandwidth. The center coordinate of
the field was set as (RA,Dec) $=$ (12h00m00.0s, $-35$d00m0.0s). For 12m (7m)
array, there are 52 (17) fields to cover the simulated image and each field is
set to be visited 30 (90) times per observing session on different days with
the integration time of 10 seconds per field. The SD observation was set up
with 169 fields, each visited 9 times per day. Consequently, the total
integration time to cover the skymodel image with 12m/7m/SD array is
1040/1020/3549 minutes, respectively. Both configurations of the ALMA-12m
array have the same integration time. This particular setup is to simulate a
realistic ALMA observation and optimise _(u,v)_ coverage. To ensure ideal
_(u,v)_ coverage even for the 7m array at each field, observations for 12m and
7m arrays are set to be carried out over four separate days with a different
hour angle at the start of observations, which is $-2$, $-1$, 0, and 1$h$. For
the SD array, observations are set to be carried out over 14 days, due to the
large number of fields. The achieved beam sizes are 57″ for the SD data and
1.2″ for the interferometric data. Relevant parameters for the simulated
observations are summarized in Table 5.

Figure 3: Maps of the simulated observations based on the skymodel image in
Fig. 2 obtained with different combination methods. From left to right and
from top to bottom: int-CLEAN (interferometric-only), SD-only, Feather, FSSC,
SDINT, and MACF. All interferometric and combined images are convolved into a
final circular beam of 1.7″, and are displayed within the same colour scale to
facilitate their comparison. The beamsize is indicated in the lower left
corner of each image.

Figure 3 shows the images obtained from the simulated skymodel observations
with interferometry-only (hereafter int-CLEAN), the SD data only, and with the
different data combination methods (Feather, FSSC, SDINT, and MACF) described
in §3 using both the simulated interferometric and the simulated SD data.

###  4.2 Observational, spectrally resolved data

Working with real observational data, one can only make assumptions about the
true sky brightness distribution. However, one condition which the combined
interferometric and SD image must at least fulfill is that it recovers most of
the flux visible in the SD image because in that image, no filtering of
spatial frequencies has taken place.

Furthermore, as mentioned in the introduction, there is a danger of the
deformation of spectral features when no data combination is applied. So, the
flux-recovery condition needs to be tested for each spectral channel
separately and the results have to be compared between channels. Considering
this, we therefore choose the spectrally resolved observation of the galaxy
M100 as an example of real observational data.

The galaxy M100 (NGC 4321) was observed by ALMA as Science Verification
(SV)777https://almascience.org/alma-data/science-verification data, and has
been featured in the CASA Guides series as an example of data combination
using the Feather method 888https://casaguides.nrao.edu/index.php/M100_Band3.
The ‘grand-design’ barred spiral morphology of this galaxy, seen relatively
face-on, reveals structures with different physical size scales. The data cube
that we present here reveals in 50 channels the CO ($1-0$) (115.271 GHz, ALMA
Band 3) molecular emission. On each side of the spectral line, additional 10
channels are included in the cube to show the spectral baseline. Observations
were made of a 47 pointing mosaic with the 12m Array, a 23 pointing mosaic
with the 7m Array, and an on-the-fly map with the SD array, effectively
recovering a range of scales down to the 12m Array resolution of $<2\arcsec$.

Figure 4: Integrated emission (moment 0, over all 70 channels) maps of our
M100 data comparing the interferometry-only (left) and the SD image (right).
The circles in the lower left represent the respective angular resolution
(beam size).

This dataset allows us to showcase the necessary preparatory steps of header
manipulation, axis reordering, and regridding. The observational details are
given in Table 5. Figure 4 shows the interferometric and the SD M100 images
(in each case the integrated emission over all 70 channels).

Table 5: Details of the datasets Parameter | Skymodel | M100
---|---|---
Phase center RA | 12h00m00.0s | 12h22m54.900s
Phase center DEC | -35d00m0.0s | +15d49m15.000s
Rest frequency | 115GHz | 115.27120GHz
$V_{LSR}$ |  $\cdots$ | 1575 km/s
$\Delta V$ (channel width) |  $\cdots$ | 5 km/s
Velocity range imaged |  $\cdots$ | [1400-1745] km/s
Map size | 204$\farcs$8 $\times$ 204$\farcs$8 | 400″$\times$ 400″
12m-array pointings | 52 | 47
7m-array pointings | 17 | 23
12m-array beamsize | 1$\farcs$2 $\times$ 1$\farcs$0 (PA=-84$\fdg$8) | 4$\farcs$0 $\times$ 2$\farcs$6 (PA=-88$\fdg$2)
Range of 12m baselines | [15 - 784] m | [15 - 200] m
Range of 7m baselines | [9 - 45] m | [9 - 42] m

##  5 Assessment metrics

Assessing the quality of the output image of a data combination is of utmost
importance in the application and interpretation of the different methods
(e.g., see Pety et al., 2013). The intrinsic definition of quality refers to
the accuracy in which the combined data reproduce the true sky emission. In
empirical research, the definition of an absolute quality assessment is
hampered by the lack of a-priori expectations for this true sky distribution.
Instead, the quality assessment of a given target image $I_{\nu}(x,y)$ can be
evaluated against the emission constraints present in a reference image
$R_{\nu}(x,y)$. Here $I_{\nu}(x,y)$ is not necessarily identical to the data
combination output image but can be a modified version of it, e.g. a version
of this image convolved to the SD resolution.

This approach is motivated by our observational products. In real observations
the emission of any interferometric or data combination map can be compared
against its SD counterpart in order to evaluate the total flux and emission
structure recovered at the SD resolution. Following this same approach,
synthetic interferometric data and combinations can also be evaluated against
simulations at the maximum resolution provided by the interferometer once both
are regridded into a common spatial and spectral coordinate system. Our
analysis examines both the synthetic data from §4.1 and the real observations
from §4.2.

We define as reference image $R_{\nu}(x,y)$ either the SD image or the
synthetic model used for comparison, while the target image $I_{\nu}(x,y)$ is
assumed to be the combination result image convolved to the resolution of the
reference image. For spectrally resolved data (cubes) our quality assessments
will also detect variations of the quality with frequency, i.e. between
channels.

We note that the effective comparison of the two images, $R_{\nu}(x,y)$ and
$I_{\nu}(x,y)$, must be performed pixel-by-pixel at the same angular
resolution (typically corresponding to one defined by $R_{\nu}(x,y)$, i.e. SD
or synthetic model convolved to low-resolution). Thus, if $I_{\nu}(x,y)|_{0}$
is the image at native interferometric resolution or the data-combined image,
then $I_{\nu}(x,y)$ refers to the image convolved and (spatially and
spectrally) regridded to the reference image parameters.

The spatial frequency filtering which occurs in interferometric observations
usually affects the structure, spectral distribution, and absolute flux of the
recovered emission. Assessing the quality of a given dataset therefore
requires considering fundamental parameters such as the total recovered flux
and its dependence on scale and intensity, as well as the variation of these
indicators across the target field. These analyses should be carried out in
statistical terms in both maps and cubes in order to make the assessment
objective. In the following sections we introduce a new set of image quality
assessment metrics for the analysis of interferometric observations and their
combination between multiple arrays and SD data, and illustrate their use on
the comparison of int-CLEAN and Feather method results. Feather is used as a
benchmark of any combination technique presented here, since it is the most
common method used in CASA. Comparison of combination methods will be done in
a following section.

###  5.1 Assessment mask: definition and adaptive thresholding

A detailed analysis of the image quality requires a careful definition of the
so-called assessment region, that is the sub-sample of (x,y) pixels (or
(x,y,$\nu$) voxels for cubes) in which our assessment metrics will be
evaluated. Not necessarily all pixels in a map are relevant for the assessment
of recovered flux. Regions devoid or showing low levels of emission should be
used specifically for determining the RMS noise of an image, as they are
dominated by pure instrumental noise, and otherwise excluded (masked) from the
image quality assessment. Ideally, any quality assessment should also consider
the structure of the astronomical object as well as include the instrumental
response. The definition of assessment areas (e.g. rectangular region) or
global intensity thresholds can simplify this analysis. To improve these basic
– and sometimes arbitrarily defined – criteria, here we introduce a more
general assessment mask $AM_{\nu}(x,y)$ depending on the properties of the
corresponding data products.

The definition of $AM_{\nu}(x,y)$ considers three conversion factors depending
on the beam sizes $\theta_{0}$ and $\theta_{R}$ of the (native)
$I_{\nu}(x,y)|_{0}$ and $R_{\nu}(x,y)$ images, respectively, as well as the
primary beam map $PB_{\nu}(x,y)$. The mask is determined according to the
following steps:

  * •

We define the “mask threshold” at a level of $3\times RMS_{mask}$, where
$RMS_{mask}$ is the off-source noise level of the target (in this case,
interferometric) image at its native resolution $\theta_{0}$.

  * •

This single RMS value is corrected for the effective antenna response in each
pixel defined by the primary-beam correction (PB) per pixel, that is,
$1/PB_{\nu}(x,y)$. The response per-pixel is therefore
$RMS_{mask}/PB_{\nu}(x,y)$.

  * •

Next, the mask threshold needs to be evaluated at the reference resolution
$\theta_{R}$ for our assessments and therefore be corrected by a factor
$\propto 1/\sqrt{N}$ where $N$ is the number of $\theta_{0}$ beams within
$\theta_{R}$, that is,
$\sqrt{\theta_{0}^{2}/\theta_{R}^{2}}=\theta_{0}/\theta_{R}$.

  * •

Finally, the resulting mask threshold value, originally in units of flux per
$\theta_{0}$ beam, needs to be converted into $\theta_{R}$ units ($\propto
1/\theta_{R}^{2}$) requiring an additional $\theta_{R}^{2}/\theta_{0}^{2}$
factor.

Altogether, we thus define our assessment mask $AM_{\nu}(x,y)$ as:

$AM_{\nu}(x,y)=\left(3\times
RMS_{mask}\right)\cdot\frac{1}{PB_{\nu}(x,y)}\cdot\frac{\theta_{R}}{\theta_{0}},$
(4)

where the $\nu$ sub-script denotes that the mask can be applied to a single
image, or plane of a cube. We apply the same $AM_{\nu}(x,y)$ mask to both
$R_{\nu}(x,y)$ and $I_{\nu}(x,y)$ images, respectively, blanking all
corresponding pixels where emission $R_{\nu}(x,y)<AM_{\nu}(x,y)$ at the
reference resolution; in other words, the so-called “mask” retains the values
where the assessment should apply. Compared to a single intensity threshold,
$AM_{\nu}(x,y)$ effectively applies an adaptive thresholding taking into
account the different image resolutions and the spatial instrumental response,
as well as channel-by-channel variations in emission structure (in the case of
cubes), where the same mask is applied to all respective images.

###  5.2 Accuracy parameter & Fidelity: assessing flux recovery

The assessment of recovered flux has been quantified in several distinct but
related ways in the literature since radio interferometers became available. A
commonly used parameter is image ”fidelity” (see, e.g. Cornwell et al., 1993,
but a detailed definition follows below). This quantity is unsigned (i.e.
doesn’t indicate if too little or too much flux is recovered) and has no
optimal value (its best value is “as large as possible”). In this work, we
introduce a new assessment parameter which we call the Accuracy parameter
(A-par). A-par is defined as:

$\begin{split}A_{\nu}(x,y)&=\frac{I_{\nu}(x,y)-R_{\nu}(x,y)}{|R_{\nu}(x,y)|},\;\;\;-\infty<A_{\nu}(x,y)<\infty\end{split}$
(5)

A-par represents the relative difference between the flux of the input image
$I_{\nu}(x,y)$ and reference image $R_{\nu}(x,y)$. Values of A-par range
between negative and positive infinity. The ideal value is zero. A-par is the
signed relative error of the recovered emission after combination with respect
to the reference image, where its sign indicates whether the emission is
underestimated ($A_{\nu}<$ 0) or overestimated ($A_{\nu}>$ 0). A perfect
matching between the input $I_{\nu}(x,y)$ and reference $R_{\nu}(x,y)$ images
thus corresponds to minimizing A-par in all pixels.

A-par can be understood as the approximate inverse of the Fidelity parameter
($F_{\nu}$), such that $A_{\nu}(x,y)\approx F_{\nu}(x,y)^{-1}$. Note that
Fidelity is defined slightly differently by different authors. A compilation
of the different conventions for ALMA, VLA, and ngVLA, can be found in Mason
(2021). We adopt the following definition:

$F_{\nu}(x,y)=\left|\frac{R_{\nu}(x,y)}{I_{\nu}(x,y)-R_{\nu}(x,y)}\right|,\;\;\;0<F_{\nu}(x,y)<\infty$
(6)

When used for image quality assessment, higher Fidelity values indicate
stronger similarities between the input $I_{\nu}(x,y)$ and reference
$R_{\nu}(x,y)$ images (if $I_{\nu}(x,y)\sim R_{\nu}(x,y)$ then $F_{\nu}\gg
0$). Fidelity becomes a powerful assessment method for the relative flux
comparison of distinct input images (e.g. $F_{\nu,1}>F_{\nu,2}$). However, as
mentioned above, the target value of this Fidelity parameter is by
construction undefined (i.e. if $I_{\nu}(x,y)\rightarrow R_{\nu}(x,y)$ then
$F_{\nu}\rightarrow\infty$), which limits its use during the assessment of
individual images in absolute terms. In comparison, A-par provides additional
information about the flux recovered (i.e. considering both A-par value and
sign), and therefore we consider it more informative. Still, our analysis
incorporates both A-par and Fidelity assessments for a full description.

Figure 5: Accuracy parameter (A-par) assessment for the int-CLEAN reduction of
our synthetic continuum skymodel data (§4.1): (upper left panel) skymodel used
as reference image $R_{\nu}(x,y)$, (upper right panel) feather image used as
target $I_{\nu}(x,y)$, (lower left panel) A-par map, and (lower right panel)
A-par histogram. Note that all results are obtained at a common resolution of
2.0″ and display all pixels within the assessment mask $AM_{\nu}(x,y)$ (in
color). Total flux values for skymodel (reference; top right) and int-CLEAN
(target; top left) images, as well as the amount of flux recovered (bottom
left) are indicated in the legends. The histogram also includes the mean and
standard deviation values for the observed A-par distribution. The vertical
dashed lines show the mean A-par (red) and its desired value zero (black).
Figure 6: Accuracy parameter (A-par) assessment for the Feather of our
synthetic skymodel data (see also Fig. 5).

We calculate both A-par and Fidelity parameters pixel-by-pixel (or voxel in
cubes) in our images, that is, $A_{\nu}(x,y)$ and $F_{\nu}(x,y)$, for all
pixels passing the assessment mask $AM_{\nu}(x,y)$ (see §5.1). We carry out
these calculations by convolving and regridding both $I_{\nu}(x,y)$ and
$R_{\nu}(x,y)$ images into a common beam $\theta$ and grid (x,y),
respectively, typically corresponding to those of the reference image
$R_{\nu}(x,y)$. Additionally, we prefer to make assessments per-pixel (Eqs. 5
and 6), in order to portray the particular regions of the image that are more
or less accurate, rather than effectively averaging over an entire image.

Figure 7:  Fidelity assessment of using Feather on our Skymodel data (§4.1).
Panels are analogous to those for A-par in Fig. 6 except that here a mean
value is not plotted since the shape of the distribution is not well
characterized by it. An ideal value does not exist by definition.

We show the A-par maps for the assessment of the data combination of our
synthetic Skymodel data (§4.1) using int-CLEAN and Feather (§3.1) in Figures 5
and 6, respectively. In both figures we display the reference (i.e. the true
Skymodel; upper left panel) and target (Feather image; upper right panel)
images, both convolved into a circular beam of 2″. An assessment mask as
described by Eq. 4 has been applied. The masked pixels are shown with white in
the upper two panels of the figure. The lower panels show the A-par map (lower
left panel) and histogram obtained from the pixel values of the A-par map
(lower right panel). Figure 7 shows similar plots but for the Fidelity
parameter. This representation allows us to intuitively identify variations of
flux recovered in different regions of our maps (A-par or Fidelity maps) as
well as to statistically describe these properties over the target field
(histograms).

Figure 8: Comparison of the A-par distributions for the interferometer-only
CLEAN (red) and Feather-combined (orange) results on our skymodel data. (Left
panel) Histogram of all the A-par values within our assessment area, where a
perfectly recovered pixel would contribute a value of 0. (Right panel) Flux
recovered per flux bin shown in individual bins (upper sub-panel) and
cumulative (lower sub-panel). In both plots, the corresponding flux of the
reference image, i.e. the expected recoverable flux, is shown in grey and
corresponds to the right axis, also in gray. In the upper right panel, the
flux recovered per bin is shown in percentage with respect to the total flux
in the image (see right vertical axis). Horizontal dashed lines, with guiding
values of 100$\pm$10% in the upper sub-panel and 50, 90, and 100% in the lower
sub-panel, respectively, are included to improve the readability of the plot.

The analysis of the histograms and distribution of the A-par values in our
images allows us to quantify the quality of a data combination technique (see
Figs. 5 & 6, lower right panel). The mean value and dispersion of such
distributions provides information about the typical fraction of recovered
flux per pixel. High quality data combination techniques would produce A-par
histograms with mean values close to zero ($A=0$) and a small width
($\sigma(A)\rightarrow 0$) showing a narrow width distribution that is
approximately Gaussian.

Similar comparisons using Fidelity would produce distributions with
increasingly high peaks and more pronounced skewness towards large Fidelity
values. In both cases the inspection of the A-par and Fidelity maps (in Figs.
5, 6, and 7) can be used to identify spatial variations and systematic effects
on the quality of the data.

Multiple histograms can be overplotted in order to compare different images
obtained from the same observation. We show an example of such a comparison in
Figure 8 (left panel) which illustrates the very significant improvement on
the recovery of the true sky emission produced by applying a data combination
technique (here Feather, blue histogram; A-par$\approx-0.08\pm 0.24$) in
comparison to using the interferometric-alone int-CLEAN reduction (red
histogram; A-par$\approx-0.60\pm 0.43$). These improvements translate into the
total amount of flux recovered in each case. While int-CLEAN only recovers
$\sim$ 32% of the total reference flux, Feather manages to improve these
values up to 89%.

###  5.3 A-par vs flux: flux recovery as a function of intensity

Figure 9: A-par values as function of intensity for both int-CLEAN, i.e. pure
interferometric data, (left) and Feather (right) methods for our skymodel
synthetic data. Lower panels: pixel-by-pixel comparison between the reference
flux $R_{\nu}(x,y)$ and recovered flux in the target image $I_{\nu}(x,y)$ for
pixels included in the assessment mask $AM_{\nu}(x,y)$ (red points). The
dashed lines indicate the expected correlations for different A-par values.
Upper panels: Individual A-par $A_{\nu}(x,y)$ values as function of the
reference flux $R_{\nu}(x,y)$. We represent the mean and standard deviation in
10 intensity bins in dark blue. Also, we display the interquartile [Q10,Q90]
range in cyan. In both lower and upper panels, we indicate the expected
distribution limits for the case of perfectly recovered emission (i.e. $A=0$)
showing only white noise deviations with $\sigma=0.1$ (in image units).

The distributions of A-par and Fidelity parameters provide information about
the global properties of the flux recovery in our data. However, the fraction
of recovered flux in the target image depends on the sky brightness
temperature. This is particularly true if the target field exhibits large
regions with low-intensity but extended flux in comparison to bright and
compact emission spots (see variations of the corresponding A-Par and Fidelity
maps in lower left panels of Fig. 6 and Fig. 7, respectively). In these cases
it is instructive to investigate not only the total A-par distribution (§5.2)
but also the variation of this parameter as a function of the reference flux.

We investigate the dependence of the flux recovery as a function of flux in
Fig. 8 (right panel). We represent both the fraction (percentage) of the flux
recovered in individual (upper sub-panel) and cumulative (lower sub-panel)
diagrams. Critical for their interpretation, we overplot the contribution of
the reference flux (again per flux bin and cumulative, respectively) in both
sub-panels. Maximizing the flux recovered is key for any data combination
method, and most significant in the case of flux bins with large contributions
to the total image flux. In our Skymodel example, bins with intensities
between $\sim$ 0.4 and 4.0 (in image units) each contribute to $>$10% of the
total flux of the reference image (see grey distribution in upper sub-panel,
with reference to the right y-axis also in grey) becoming the most relevant
flux bins.

In Figure 9 we show a direct pixel-by-pixel comparison between reference flux
$R_{\nu}(x,y)$ in the Skymodel and the recovered flux $I_{\nu}(x,y)$ by our
Feather combination within $AM_{\nu}(x,y)$ in the lower panel of this figure.
Similarly, we show the variation of the A-par values as function of
$R_{\nu}(x,y)$ in the upper panel, including the mean and sigma values (dark
blue) and inter-quartile [Q10,Q90] limits (cyan) in 10 flux bins. In this
representation, a perfect data combination should produce a linear correlation
between the $R_{\nu}(x,y)$ and $I_{\nu}(x,y)$ fluxes (i.e. a straight line in
the lower panel) with a narrow distribution independent on the reference flux
(i.e. A-apar consistent with noise) in the upper panel.

The analysis of these diagrams permits a direct evaluation of the dependence
of the flux recovery with intensity in real data combination techniques such
as the one shown in Figure 9. Our quality assessments demonstrate how pure
interferometric data fails to recover the expected reference flux distribution
at all intensities (A-par $<0$ in all flux bins; see Fig. 9, left panel).
These issues also affect regions with highest intensity and compact emission
with average flux losses of ca. 25%. More importantly, the performance of pure
interferometric data processed with CLEAN shows a strong dependence of A-par
on flux, producing worse results at lower signal values. In comparison,
Feather shows a clear overall improvement on the flux distribution recovery at
all intensities with A-par values $\approx 0\pm 0.5$ in all our signal bins,
as expected from the analysis of the total flux recovery ($A=-0.07$ or -7%
mean offset, see Fig. 8). Despite these improvements, our analysis reveals an
increase of the dispersion of A-par values at lower signal values indicating
that, although better than a pure CLEANed interferometric data, the
performance of Feather may be hampered in regions of low intensity emission
(see a more detailed discussion in § 6).

The enhanced image quality indicated by A-par both globally (Fig. 8) and per
signal bin (Fig. 9) corresponds to an improved flux recovery by a combination
technique, in this case Feather (e.g. see Fig. 6), compared to
interferometric-only maps. An ideal data combination should in principle
maximize the flux recovery at all scales. In practice, however, this may not
be possible or may not be critical since, depending on the emission
distribution, not all signal bins may contribute equally to the total flux
budget. This is illustrated in Figure 8 (right panel) where we show the
percentage of flux recovered per signal bin independently (upper sub-panel)
and cumulatively (lower sub-panel) respect to the reference image. The dashed
grey line in the upper sub-panel indicates the contribution of each individual
signal bin to the total flux budget. We point out that for the skymodel data,
bins with intensities of $>0.2$ (in image units) contribute to most of the
emission within the assessment area ($>80$ % of the total). Recovering the
flux in those bins therefore guarantees the recovery of a significant fraction
of the total image sky brightness temperature. The analysis of these plots
also highlights the poor performance of int-CLEAN producing negative flux
values (due to PSF sidelobes plus filtering effects) in regions with extended
emission showing emission values $<1$ (in image units) and leading to flux
recovery values of only 25% of the total flux.

The high sensitivity of the above comparisons demonstrates the robustness of
our assessment metrics. Additional tests indicate that deviations at different
flux intensities may also reveal potential issues in the data combination
process (e.g. recovery of extended emission) and/or limitations of such
procedures (e.g. lack of baselines overlap). Moreover, this analysis permits
to investigate the performance of different data combination techniques and
their response to images with different dynamic range in emission (see
discussion of the data combination results in §6).

###  5.4 Power Spectra (PS): spatial scale sensitivity

Figure 10:  Spatial Power Spectra (SPS) obtained from the pure interferometric
data (int-CLEAN, red) and the Feather-combined interferomentric and SD data
(orange) compared with the original power spectrum of our skymodel simulation
(dark blue). Note the logarithmic scales and the several orders of magnitude
of flux present in the input model image. Vertical lines show the final image
resolution (dashed line, right) and largest angular scale effectively sampled
by the SPS (1/3 of the map size; dot-dashed line, left), beyond which this
simulated instrument is not sensitive.

The primary goal of each data combination method is the successful recovery of
emission at the angular scales filtered out by the interferometer. Thus, it is
crucial to describe how the recovered emission is distributed as a function of
angular scale. One such tool is the spatial power spectrum (SPS) defined as

$P\left(k\right)=\mathscr{F}\left[I_{\nu}\left(x,y\right)\right]\mathscr{F}^{*}\left[I_{\nu}\left(x,y\right)\right],$
(7)

where $k$ is the wavenumber ($k=\frac{1}{\theta}$ where $\theta$ is the angle
on the sky in rad) and $\mathscr{F}$ is the Fourier transform of the input
image, $I_{\nu}\left(x,y\right)$, under study. This is a two-point correlation
function that measures how power (i.e., structure) is distributed across
spatial scales. In practice the distribution of power is measured by computing
the 2D Fourier transforms of the integrated intensity images and measuring the
median in progressively larger annuli. We chose the median to mitigate the
bias introduced by ringing artifacts of the sinc function along the axis of
the 2D FT when emission extends toward the edge of the image. While phase
information is lost due to the multiplication by the complex conjugate,
comparing the SPS profiles derived from the images produced by the various
combination methods to the input reference (or total power) demonstrates how
successfully the emission is recovered across the full range of angular
scales. The ideal image would recover the flux of the input model image on all
scales.

As an example, we present in Fig. 10 the SPS profiles measured for the data
combination of our skymodel using Feather in comparison to the interferometer-
only data processed with CLEAN. Also shown is the profile of the original
input image of the simulation (Reference). These three profiles agree
remarkably well at larger angular scales (left side of the plot). Feather and
reference continue to agree all the way down to the smaller angular scales
accessible to the instrument. But the int-CLEANed interferometer-only data
start to diverge from the reference as the scale size decreases. This result
can be understood by considering the way the power is measured:  The
insufficient recovery of the larger scales effectively boosts the intensity
contrast of the already existing small-scale structure in the interferometer-
only image leading to higher power (i.e. larger differences) at small angular
scales compared to the reference image. Conversely (not shown in this
example), when the SD data is given too much weight in the combination, the
power decreases at small scales because the low resolution SD data effectively
washes out the observed small-scale structure. The variation in power at small
scales highlights subtle deficiencies that may not be immediately and clearly
visible in the images. Scales below the beam size (right of the dashed line in
the figure) are dominated by correlated noise in the pixels or convolution
residuals and should generally be ignored.

Figure 11:  A-par Spatial Power Spectra (A-par SPS) for the same int-CLEAN and
Feather data reductions shown in Fig. 10. For comparison we add the expected
(flat) A-par power spectrum of our reference image with Gaussian noise added.

The above SPS analysis can also be applied to the case of our A-par maps (e.g.
Fig. 6, lower left panel). The expected A-par power spectrum provides
information on the variance of the relative errors, $\sigma$(A-par), as a
function of scales (e.g. high A-par power denotes large systematic errors at a
given scale). In addition to the minimization of the (global) A-par mean and
dispersion values (e.g. in Fig. 8), the goal of any data combination method is
therefore to minimize the A-par SPS, that is, all systematic errors as a
function of scale. We show the results of the A-par SPS for our int-CLEAN and
Feather reductions in Fig. 11. Our results illustrate the reduction of the
A-par variance at all scales once a Feather combination is applied to the
original int-CLEAN reduction. The interpretation of these results in
comparison with other combination methods and the addition of Gaussian noise
to our reference skymodel is discussed in § 6.4.

###  5.5 Using single-dish observations as reference image

Unlike the case of our synthetic Skymodel as presented in our previous
assessments, the quality assessment of real scientific datasets is usually
hampered by the lack of information for reference at high spatial resolution.
Instead, when considering total emission on all scales, interferometric
observations can only be cross-checked with the flux measurements obtained at
SD resolution. Such comparison requires the original interferometric image at
its native (high-)resolution $I_{\nu}(x,y)|_{0}$ to be convolved into a final
(low-)resolution image $I_{\nu}(x,y)$ similar to the SD data $R_{\nu}(x,y)$.

We quantify the effects of this low-resolution comparison in our assessment
metrics in Figure 12 showing the results of int-CLEAN (red) and Feather (blue)
for our skymodel simulation, this time convolved into the SD resolution (57″),
for the case of the A-par statistics (left panel) and the recovered flux-per-
signal-bin (right panel). Figure 12 can be directly compared with Fig. 8
carried out at much higher resolution (1.7″). Several results can be taken
from this comparison. The number of pixels for the quality assessment is
largely reduced due to the use of a much larger beamsize as seen in the
overall pixel statistics. Changes in the beamsize also alter the local A-par
estimates due to the averaging of regions with/without emission within the
much smoother SD beam. Nonetheless, our assessment metric preserves the
relative differences between methods. Both the A-par statistics (left panel)
and the recovered flux (right panel) indicate the improved performance of data
combination techniques such as Feather in comparison with interferometric-only
CLEAN deconvolution.

Several caveats should be considered during the quality assessment of images
at the SD resolution. First, and by construction, the SD observations only
provide information of the total flux per SD beam but the emission
distribution at smaller scales remains unconstrained. Second, large
convolutions (such as the one shown in Fig. 12) can lead to averaging effects
of multiple emission features at interferometric resolutions (e.g. compact
sources and negative sidelobes) unresolved at the SD beamsize. Third, the
limited amount of resolution elements in the SD images restricts the use of
some size-dependent assessments such as SPS. Treated with caution at low-
resolutions, however, our assessment metrics provide a powerful toolkit to
assess the quality of both continuum and spectral interferometric
observations.

Figure 12: Same as Figure 8 but using the SD image as reference rather than
the simulation input image as it would be done when working with real
observational data.

###  5.6 A-par spectrograms: analysis of spectral cubes

Figure 13: A-par “spectrogram” for our M100 data obtained from the
interferometric data alone with int-CLEAN (top) and from the Feather combined
interferometric data and SD image cube (bottom) for all voxels within our
assessment mask. The detected emission in each spectral channel of this CO
line observation represents emission from gas at a specific velocity range
with its own spatial distribution.

In order to extend the previous assessments, which were of a continuum (single
channel) image, to a data cube, we must visualize the variation in the spatial
distribution as a function of frequency (or velocity) channel. This is
important because the astronomical objects are likely to have distinct
morphology for material at different velocities. Hence, any combination
method, or perhaps no combination, might be sufficient for the channels with
more compact emission, whereas the channels with more extended emission depend
much more on the combination for an acceptable image quality.

We extend the A-par analysis from the previous continuum case (e.g. Figure 10)
to that of a cube in the form of an “A-par spectrogram”, shown in Figure 13.
This spectrogram shows the distribution of A-par values for all the channels
in the target cube that satisfy our assessment mask $AM_{\nu}(x,y)$ criteria
(left panel). For a given velocity channel (x axis), the spectrogram shows the
number of voxels within this frequency plane (color scale) with a given A-par
value (y axis). Its analysis provides a visual comparison of the emission
recovered in multiple channels and allows the evaluation of local (per
channel) and global (mean distributions) deviations as well as issues (e.g.
outliers and/or problematic channels) in the data. Similar to the analysis of
the continuum images, the goal of any data combination method is to minimize
the value of this “A-par spectrogram” for all channels within a cube. We
complement these results with additional histograms (right panel) for A-par
showing all individual channel values (grey histograms) as well as the mean
distribution (red histogram).

These spectrograms turn out to be a powerful and intuitive tool when comparing
multiple data combination techniques in large spectral cubes. In Fig. 13 (top
panel) we illustrate again the limitation of the interferometer-only emission
recovery (A-par$<$ 0) at different velocities/channels. Moreover, our
spectrograms illustrate how the flux losses also depend on velocity, since the
brightness distribution also varies with velocity (see 6.2). Depending on the
selected channel, A-par varies between $\sim$ -0.1 (10% flux losses) and $>$
-0.5 ($>$ 50% flux losses). In contrast, data combination with Feather (lower
panel) produces significantly better results at all velocities/channels with
typical values of $|$A-par$|\sim$ 0.2 (or $\pm$ 20% error) for most channels
(between $[10,60]$). The improvement is also seen in the mean A-par vaues in
the corresponding histograms. The broad and shifted distribution in the int-
CLEAN reduction with mean 50% flux losses (A-par$=-0.51\pm 0.07$) is
significantly improved by Feather showing a narrower and more centered
distribution showing differences of less than 10% respect to the SD flux per
channel (A-par$=-0.07\pm 0.04$).

A variation of A-par with channel directly translates into a deformation of
spectral line features as will be discussed in §6.2. That is the case for our
int-CLEAN reduction as seen in Fig. 13 (top panel). In other words, data
combination is also relevant for spectroscopy.

##  6 Combination Method Comparison and Discussion

###  6.1 The indispensable short-spacing information

Figure 14: Comparison of the true sky brightness distribution with the
recovered emission of interferometry-only and two different data combination
methods in the central region of our Skymodel: (top left) Skymodel image (i.e.
reference), (top right) int-CLEAN (i.e. only interferometric data), (bottom
left) Feather, and (bottom right) SDINT, all convolved to a common resolution
of 1.7″ (see beamsize in the lower left corner) and shown with the same
intensity scale.

The results presented in § 5 demonstrate the large impact of the short-spacing
information on the analysis of interferometric observations. All assessment
metrics show significant improvements of the flux recovered at different
spatial scales and input signals using a standard data combination technique
such as Feather when compared to pure interferometric-only CLEAN (int-CLEAN)
reductions (see Figs. 6-13). Here, we explore the relative performance of
additional data combination techniques.

Figure 15: Compact emission recovered by different data combination methods at
the central region of skymodel image all convolved into a final resolution of
2 ″(see bottom left corner of each image) and plotted within the same
intensity range and with the same contours every 0.5 Jy beam-1. From left to
right and from top to bottom: skymodel, Int-CLEAN, Feather, SDINT, and MACF
results. The lower right panel shows the intensity profile indicated with a
dotted line in our maps obtained by these methods.

We highlight the enhanced image quality of different data combination
techniques in Figure 14 by comparing our input Skymodel emission (reference;
top left panel) with the results obtained by int-CLEAN (top right panel),
Feather (bottom left panel), and SDINT (bottom right panel) in the central
region of the Skymodel image (see also Fig. 3; FSSC and MACF maps are not
shown as they produce visually similar results to Feather and SDINT,
respectively). This region is selected because it incorporates both compact
and extended emission. A direct inspection of these images reveals how pure
int-CLEAN deconvolutions (top right panel) miss large fractions of the true
sky emission at all scales and flux values due to the incomplete sampling of
the _(u,v)_ -space in interferometric-alone observations when extended sky
emission is present. The resulting interferometric filtering becomes non-
intuitive, showing highly non-linear effects in the case of complex sources
and leading to a flux recovery of only $\sim$ 10-20% of the true sky emission
(see Fig. 8, red lines). These effects can be sometimes identified in our
images by the appearance of regions with negative emission (negative
sidelobes; dark areas in the int-CLEAN image)999We note here that the absence
of negative emission troughs in the image does not guarantee the absence of
filtering, as the effects of filtering can be compensated by the presence of
additional positive emission depending on the source structure..

Interferometric-alone observations (int-CLEAN) are expected to efficiently
recover the total flux of unresolved, isolated sources with sizes comparable
to the interferometer beamsize. However, this situation is not the same in the
presence of smooth, extended emission. In the latter cases, it is sometimes
assumed that the lack of short-spacing information only affects the recovered
emission beyond the Maximum Recoverable Scale (MRS)101010See Remijan et al.
(2019) for a definition. The MRS depends on the details of the interferometric
array setup. For the ALMA 12m nominal configurations with 43 antennas, the MRS
is approximately a factor 10 larger than the achieved angular resolution.
sampled by the interferometer. But our results clearly demonstrate how the
interferometric filtering has a significant impact at all scales including
those compact regions with sizes comparable to the beamsize (i.e. angular
scale $\ll$ MRS).

The addition of the short-spacing (SD) information dramatically improves the
image quality. Standard techniques such as Feather (Fig. 14, lower left panel)
recover up to $\sim$ 90% of the input Skymodel emission (see Fig. 8, blue
lines) and remove most of the negative emission patches in the field of
interest.

The addition of the short-spacing information at low-resolutions also has a
direct impact on the peak flux in many of the compact sources in this region,
potentially affecting column density and mass estimates in subsequent
scientific analysis. We demonstrate these issues in Fig. 15 by comparing the
emission recovered by int-CLEAN, Feather, SDINT, and MACF at the central
region of our skymodel showing multiple compact emission features with sizes
comparable to the beamsize. For illustrative purposes, we compare the emission
profile of the brightest region within this field (see lower right panel).
Although compact, the peak flux of this source in the interferometric-alone
int-CLEAN observations is reduced by a factor of $\sim$ 40% compared to the
true emission in our skymodel. Flux losses are also seen in other nearby
sources together with negative sidelobes in regions with fainter emission (see
also upper right panel). Still, the effects on our maps depend on the location
of the source and the structure of their surroundings, making impossible a
direct correction without the zero-spacing information. The use of data
combination techniques significantly improve the flux recovery in these
compact regions with, in our particular case, SDINT and MACF producing the
best results. These reported differences could have a significant impact on
derived properties such as mass and column densities of compact sources in
complex environments. These results underline the absolute need for short-
spacing information and data combination in the analysis of particular high-
dynamic range observations provided by state-of-the-art interferometers such
as ALMA.

###  6.2 Effects of the short-spacing information in line profiles

Figure 16: Examples of the data combination on lines profiles obtained from
our analysis of two representative regions of our M100 dataset: (Left panel)
Line profile in one region shows effects on the line peaks and relative
intensities; (Right panel) Line profile in another, lower-flux region shows
effects on the line shape.

So far the discussion focused on the analysis of the integrated emission maps.
As already introduced in § 5.6, the extent of the effects of lack of short-
spacing information will depend on the target structure in each frequency
channel. We illustrate these effects in Fig. 16 in two representative regions
in our M100 dataset. As shown in Fig. 16 (left panel), the interferometric
filtering produces changes of up to 50% on the flux per channel in this
spectrum. The true spatial structure of the emission and thus the fraction of
missing observed flux is however not necessarily homogeneous but can vary from
channel to channel. According to Fig. 13 (upper panel), these losses can be as
high as 90% in some channels. This can alter line ratios by several ten
percent, an effect already reported in previous studies (e.g., Pety et al.,
2013).

Figure 16 (right panel) shows how frequency-dependent, spatial filtering
effects can also modify the line Full Width at Half Maximum (FWHM) and line
centroid of the observed gas components. As illustrated in the int-CLEAN
spectrum, negative sidelobes can selectively alter different channels (see
negative feature in the spectrum) both reducing the line FWHM and shifting the
line peak in comparison with the Feather and SDINT reductions. Data
combination is therefore essential to recover not only the total flux in
spectral cubes but also to preserve the kinematic information in them. For
wider frequency coverage, one should also ensure that the frequency dependence
of _(u,v)_ spacing and PB are taken into account.

###  6.3 Comparison between data combination methods

Figure 17: Comparison of the data combination methods for the case of our
Skymodel: int-CLEAN (red), Feather (orange), SDINT (yellow), and MACF (light
blue). (Left panel) A-par statistics. (Right panel) Flux recovered by flux bin
compared to our skymodel as reference (grey triangles). See also Fig. 8, as
gray horizontal lines are the same.

Inspecting Fig. 14 in detail suggests that more elaborate methods such as
SDINT (lower right panel) may be able to reproduce the original sky brightness
distribution (upper left panel) even better than Feather. We quantitatively
investigate the performance of the different data combination techniques
introduced in § 3, all of which are available in the script we provide. In
Fig. 17 (left panel) we show the overall A-par values as well as the total
recovered flux in the images (see also figure caption) obtained by the int-
CLEAN, Feather, SDINT, and MACF methods, as applied to the Skymodel. The FSSC
method is not shown but produced results nearly identical to Feather in this
test. Feather recovers $\sim$ 89% of the total flux with a mean accuracy of
A-par $=-0.08$ and an accuracy dispersion of $\sigma($A-par$)=0.24$. The total
flux recovery is improved by the SDINT and MACF methods with total flux values
above $\sim$ 92% with respect to the true sky emission. While both are better
than Feather, we notice quantifiable differences between the SDINT and MACF
methods. While SDINT shows a smaller dispersion of the A-par values
($\sigma($A-par$)=0.13$, seen as a narrower distribution in this histogram),
MACF shows a better mean A-par value (A-par = -0.08, leading to a more
centered distribution). These differences are likely connected to the slightly
different responses of these methods to the flux intensity (see Figure 17,
right panel). Compared to Feather (and FSSC), both SDINT and MACF maximize the
flux recovery in those high intensity flux bins contributing to most of the
total sky emission.

###  6.4 Flux and A-par power spectra

Figure 18: Power spectrum of image intensity as introduced in § 5.4 and Fig.
10, here also showing the results for SDINT (yellow) and MACF (light blue)
with our skymodel as reference (dark blue). As in Fig. 10, vertical lines show
the range above and below which this simulated instrument is not sensitive.
Figure 19: Power spectrum of the A-par map as introduced in § 5.4 and Fig. 11,
here also showing the results for SDINT (yellow) and MACF (light blue).

Continuing the discussion in § 5.4 we show in Figs. 18 and 19 the intensity
power spectrum and the A-par power spectrum for the results on our simulated
Skymodel observation obtained by interferometry only and the different
combination methods.

The intensity power spectrum (Fig. 18) of all combination methods traces quite
closely that of the reference at first glance but with small variations.
Between the largest scales sampled by our SPS (left vertical dot-dashed line)
and log(angular scale)$\sim$0.7, all results are on or above the reference
with SDINT showing the closest power. Between log(angular scale)$\sim$0.7 and
$\sim$ 0.4, the results are constantly below the reference. Here, MACF is
closer to the reference. Shortly before reaching the smallest accessible
angular scales (right vertical dashed line) the results rise again above the
reference. And in this range SDINT performs best. Although tracing scales
rather than only flux recovery, the results of these SPS are consistent with
the analysis of Fig. 17 (right panel).

The A-par map power spectra shown in Fig. 19 speak a clearer language because
in this type of plot, the deviations from the reference occur in the angular
scale range where the instrumental sensitivity is actually missing. All
results with or without data combination, and no matter which method is used,
are far above the reference and only come close to it at the smallest angular
scales (right side of the plot). This is because the reference compared to
itself of course has an A-par of zero and the non-zero A-par values in the
plot are only caused by the added (Gaussian) noise which has no spatial
structure. The deviations of the reconstructed images from the reference,
however, do have a spatial structure, and A-par is extremely sensitive to this
(see also Fig. 6, lower left panel). In other words, the slopes seen in these
A-par SPS compared to a pure white noise (dark blue line in the plot)
demonstrate how the deviations of any data combination method compared to the
original reference image are not Gaussian but depend on scale and intensity
(see distribution of the A-par variations in Fig. 6, lower left panel). With
this sensitivity, we can see that SDINT is performing increasingly better
(lower, flatter) than the other combination methods as the angular scale
increases (towards the left). This corresponds to the narrowness of the A-par
distribution seen in the previous section.

Our experiments demonstrate how interferometers can significantly alter the
original power spectrum of a science target such that in the Skymodel. The
incomplete sampling of the (u,v)-plane information introduces non-linear
effects at different scales (see A-par SPS in Fig. 18) modifying both the
absolute value and slope of the recovered SPS that persist even after data
combination (see Fig. 19). Similar results are found in other structure
analysis techniques such as Probability Distribution Functions (PDF; e.g. see
Ossenkopf-Okada et al., 2016). The physical interpretation of the recovered
SPS in interferometric observations (e.g. absolute values and changes in the
power spectrum slope) should be treated with extreme caution.

###  6.5 Single-dish flux recovery

Figure 20: A-par spectrograms for SDINT (upper panel) and MACF (lower panel)
methods for the case of M100. See also Fig. 13.

In § 5.5 we explained that for actual observational data, the SD image needs
to serve as the reference image. Fig. 13 already showed the A-par spectrograms
for the int-CLEAN image cube and the Feather combined image cube of our M100
data using the M100 SD image cube as reference. In Fig. 20 we show the
corresponding plots for the remaining combination methods SDINT and MACF.
Similar to the case of our continuum maps (§ 6.1), all data combination
methods produce significant improvements in the image quality of the spectral
cube.

In the case of our M100 dataset all methods produce satisfactory results with
deviations from the expected flux per channel reduced to $\lesssim$10% in all
cases. The best results are produced by Feather and MACF methods both showing
narrow A-par distributions ($\sigma($ A-par$)=0.04$) with an average flux
recovery of $\gtrsim$ 95% ($|$A-par$|<$ 0.05) while SDINT shows slightly
larger differences with respect to the SD flux ($|$A-par$|=-0.11$ and
$\sigma($ A-par$)=0.06$). Comparisons within these 10% variations are however
difficult to interpret given the low number statistics available ($\sim$ 20
voxels per channel within the assessment mask) once these analyses are carried
out at the SD resolution.

Compared with the results obtained on our Skymodel, where SDINT and MACF
produced the best reductions (§ 6.3), our analysis suggests that the
performance of each combination method may slightly vary depending on spatial
structure of the astronomical target, which may vary between velocity
channels; there may be additional dependence on the depth and type of
observation. While our analysis does not allow to establish the clear
superiority of any of the data combination methods in all situations, our
assessment metrics appear to be a useful diagnostic toolkit to systematically
investigate and compare their performance, statistically quantify their
ability to recover the true sky emission at different scales, and identify
potential combination issues in both continuum and spectral datasets. For
simplicity, our data combinations assumed standard values for reduction
parameters such as briggs (CLEAN), sdfactor (Feather and MACF) or sdgain
(SDINT), among others (see § 2 and 3). The assessments presented here can also
be extended to quantify the effect and sensitivity on the flux and scales
recovered by other parameter choices and techniques. Users of these methods
are therefore recommended to evaluate their results using these diagnostics.

##  7 Summary

In this paper we underline with synthetic and real-world examples the
importance of the combination of single-dish (SD) data with interferometric
data in radio astronomy when the objects observed emit on scales that are
larger than a few times the angular resolution of the interferometric
observation. We then briefly describe the prominent data combination methods
available today (Feather, TP2VIS, MACF, SDINT, FSSC) and test them on two
datasets, one simulated and one observational, with a new set of tools
producing several novel diagnostic plots which permit one to quantify and
compare the performance of the combination methods. These tools are mostly
based on the accuracy parameter, “A-par,” which is essentially the relative
difference between the output image of the method and the original input image
of the simulation (in the case of simulated test data) or the original SD
image (in the case of a real observation). We also discuss the commonly used
Fidelity parameter, and we find A-par to have higher diagnostic value. An
additional valuable diagnostic (closely related to A-par) is the fraction of
recovered flux as a function of reference flux either per reference flux bin
or cumulative.

For a given observation, the A-par can be inspected in four different ways: as
a spatial map, as a histogram of the A-par value for all pixels, as a power
spectrum of the A-par map (as a function of angular scale), and for spectrally
resolved data (image cubes) as a spectrogram which shows the per-channel A-par
histograms for each channel next to each other. For individual images (a
single spectral channel) the A-par maps and histograms can clearly reveal the
need for data combination. For spectrally resolved data, we demonstrate that
not applying data combination can result in significantly deformed spectral
line shapes and intensity ratios and this is clearly visible in the A-par
spectrograms.

As a result of our analysis, we recommend obtaining additional SD observations
for all interferometric observations of significantly resolved objects amidst
extended emission, and to produce diagnostic A-par plots in order to verify
(a) whether combination with SD data is required and (b) whether the
combination was performed optimally.

Our results on the performance of the tested data combination methods show
quantitatively that the application of any method is better than using the
interferometric data on its own. Interferometric-alone reductions can miss up
to 90% (A-par=-0.9) of the flux depending on the scale and source morphology.
These flux losses are present at all scales, affecting both the diffuse
emission and compact sources down to scales similar to the beamsize. In
contrast, the use of data combination techniques such as Feather, SDINT, or
MACF allows the recovery of most of the true sky emission within errors of
$\lesssim$ 10% ($|$A-par$|<$ 0.1). Our tests show that the enhanced
performance of some of these methods under certain circumstances can reduce
the mean value ($|$A-par$|<$ 0.05) and variance ($\sigma($A-par$)<$ 0.05) of
these differences to $<$ 5% with respect to the actual target flux.

All combination methods have a parameter which permits tuning the contribution
of the SD data, the choice of which depends on the relative sensitivity
achieved in the interferometric and SD observations. While this parameter
should be chosen on first principles, one can also tune the relative weighting
of the input data depending on their sensitivities and the scientific goal of
the experiment. Moreover, the goodness of these different combinations can be
further quantified with our assessment diagnostics.

Recovering the true sky emission at all scales has fundamental implications
for the physical interpretation of scientific images. For instance, the
continuum flux is directly connected to measurements of the (gas and dust)
column densities, probability distribution functions, spectral indices, and
masses. Similarly, the line emission profiles provide information on the gas
excitation conditions, molecular abundances, or gas chemistry, among others.
The addition of the zero-spacing information is revealed as fundamental to
produce accurate, high-dynamic range images in radio-astronomical,
interferometric observations of spatially resolved targets with complex
emission substructure. Our results demonstrate the importance of implementing
advanced data combination techniques in current (e.g., NOEMA, SMA, and ALMA)
and future (e.g., ngVLA and SKA) interferometric studies. The reader is
encouraged to try our demonstration software package available from
https://github.com/teuben/DataComb.

The authors gratefully acknowledge the thorough review and constructive
comments from the very knowledgeable anonymous referee. The authors
acknowledge the Lorentz Center, which enabled the launch of this initiative
through funding, organizing and hosting the workshop “Improving Image Fidelity
on Astronomical Data: Radio Interferometer and Single-Dish Data Combination”
held on 12 - 16 August 2019. The European Southern Observatory provided
additional support. We particularly thank Hauyu Liu, Yanett Contreras, and
Alvaro Sanchez-Monge for co-organizing this 2019 Lorentz Center workshop.
Including, and in addition to, participants of the Lorentz Center workshop, we
acknowledge in particular: Amanda Kepley, Urvashi Rau, Sandra Burkutean, Adam
Ginsburg, and Jens Kauffmann, who contributed to development and discussion of
techniques used in this manuscript. This paper makes use of the following ALMA
data: ADS/JAO.ALMA#2011.0.00004.SV. ALMA is a partnership of ESO (representing
its member states), NSF (USA) and NINS (Japan), together with NRC (Canada),
MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. The National Radio Astronomy Observatory is a facility of the National
Science Foundation operated under cooperative agreement by Associated
Universities, Inc. AH and SS acknowledge the support and funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement Nos. 851435). AH
acknowledges assistance from Allegro, the European ALMA Regional Center node
in the Netherlands. LMF acknowledges funding by the German Federal Ministry of
Education and Research under the funding code 05A20PD1. YM acknowledges
support from grant JSPS KAKENHI (Grant Number JP20K04034). KMH acknowledges
financial support from the grant SEV-2017-0709 funded by MCIN/AEI/
10.13039/501100011033, from the coordination of the participation in SKA-
SPAIN, funded by the Ministry of Science and Innovation (MCIN), and financial
support from the grants RTI2018-096228-B-C31 and PID2021-123930OB-C21 funded
by MCIN/AEI/ 10.13039/501100011033, by “ERDF A way of making Europe” and by
the “European Union”. This research made use of Astropy
(http://www.astropy.org), a community-developed core Python package for
Astronomy (Astropy Collaboration et al., 2013, 2018).

## References

  * Arras et al. (2021) Arras, P., Bester, H. L., Perley, R. A., et al. 2021, A&A, 646, A84, doi: 10.1051/0004-6361/202039258
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Beuther et al. (2022) Beuther, H., Wyrowski, F., Menten, K. M., et al. 2022, A&A, 665, A63, doi: 10.1051/0004-6361/202244040
  * Briggs (1995) Briggs, D. S. 1995, PhD thesis, New Mexico Institute of Mining and Technology
  * Carrillo et al. (2014) Carrillo, R. E., McEwen, J. D., & Wiaux, Y. 2014, MNRAS, 439, 3591, doi: 10.1093/mnras/stu202
  * Clark (1980) Clark, B. G. 1980, A&A, 89, 377
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693, doi: 10.1086/300337
  * Cornwell (2008) Cornwell, T. J. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793, doi: 10.1109/JSTSP.2008.2006388
  * Cornwell et al. (1993) Cornwell, T. J., Holdaway, M. A., & Uson, J. M. 1993, A&A, 271, 697
  * Cotton (2017) Cotton, W. D. 2017, PASP, 129, 094501, doi: 10.1088/1538-3873/aa793f
  * Dirienzo et al. (2015) Dirienzo, W. J., Brogan, C., Indebetouw, R., et al. 2015, AJ, 150, 159, doi: 10.1088/0004-6256/150/5/159
  * Faridani et al. (2018) Faridani, S., Bigiel, F., Flöer, L., Kerp, J., & Stanimirović, S. 2018, Astronomische Nachrichten, 339, 87, doi: 10.1002/asna.201713381
  * Gildas Team (2013) Gildas Team. 2013, GILDAS: Grenoble Image and Line Data Analysis Software, Astrophysics Source Code Library, record ascl:1305.010.  http://ascl.net/1305.010
  * Hacar et al. (2018) Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77, doi: 10.1051/0004-6361/201731894
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jiao et al. (2022) Jiao, S., Lin, Y., Shui, X., et al. 2022, Science China Physics, Mechanics, and Astronomy, 65, 299511, doi: 10.1007/s11433-021-1902-3
  * Jorsater & van Moorsel (1995) Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037, doi: 10.1086/117668
  * Kauffmann et al. (2017) Kauffmann, J., Pillai, T., Zhang, Q., et al. 2017, A&A, 603, A89, doi: 10.1051/0004-6361/201628088
  * Kepley et al. (2020) Kepley, A. A., Tsutsumi, T., Brogan, C. L., et al. 2020, PASP, 132, 024505, doi: 10.1088/1538-3873/ab5e14
  * Koda et al. (2019) Koda, J., Teuben, P., Sawada, T., Plunkett, A., & Fomalont, E. 2019, PASP, 131, 054505, doi: 10.1088/1538-3873/ab047e
  * Koda et al. (2011) Koda, J., Sawada, T., Wright, M. C. H., et al. 2011, ApJS, 193, 19, doi: 10.1088/0067-0049/193/1/19
  * Kong et al. (2018) Kong, S., Arce, H. G., Feddersen, J. R., et al. 2018, ApJS, 236, 25, doi: 10.3847/1538-4365/aabafc
  * Kurono et al. (2009) Kurono, Y., Morita, K.-I., & Kamazaki, T. 2009, PASJ, 61, 873, doi: 10.1093/pasj/61.4.873
  * Leroy et al. (2021) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43, doi: 10.3847/1538-4365/ac17f3
  * Mason (2021) Mason, B. 2021, NGVLA Memo 67
  * Mason (2020) Mason, B. S. 2020, arXiv e-prints, arXiv:2006.06549.  https://arxiv.org/abs/2006.06549
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606, doi: 10.1093/mnras/stu1368
  * Ossenkopf-Okada et al. (2016) Ossenkopf-Okada, V., Csengeri, T., Schneider, N., Federrath, C., & Klessen, R. S. 2016, A&A, 590, A104, doi: 10.1051/0004-6361/201628095
  * Peretto et al. (2013) Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112, doi: 10.1051/0004-6361/201321318
  * Pety (2005) Pety, J. 2005, in SF2A-2005: Semaine de l’Astrophysique Francaise, ed. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 721
  * Pety & Rodríguez-Fernández (2010) Pety, J., & Rodríguez-Fernández, N. 2010, A&A, 517, A12, doi: 10.1051/0004-6361/200912873
  * Pety et al. (2013) Pety, J., Schinnerer, E., Leroy, A. K., et al. 2013, ApJ, 779, 43, doi: 10.1088/0004-637X/779/1/43
  * Pineda et al. (2011) Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2011, ApJ, 739, L2, doi: 10.1088/2041-8205/739/1/L2
  * Pingel et al. (2022) Pingel, N. M., Dempsey, J., McClure-Griffiths, N. M., et al. 2022, PASA, 39, e005, doi: 10.1017/pasa.2021.59
  * Plunkett et al. (2015) Plunkett, A. L., Arce, H. G., Corder, S. A., et al. 2015, ApJ, 803, 22, doi: 10.1088/0004-637X/803/1/22
  * Plunkett et al. (2013) —. 2013, ApJ, 774, 22, doi: 10.1088/0004-637X/774/1/22
  * Rau & Cornwell (2011) Rau, U., & Cornwell, T. J. 2011, A&A, 532, A71, doi: 10.1051/0004-6361/201117104
  * Rau et al. (2019) Rau, U., Naik, N., & Braun, T. 2019, AJ, 158, 3, doi: 10.3847/1538-3881/ab1aa7
  * Remijan et al. (2019) Remijan, A., Biggs, A., Cortes, P., et al. 2019, ALMA Technical Handbook Cycle 7 (ALMA Doc. 7.3, ver. 1.1)
  * Rodríguez-Fernández et al. (2008) Rodríguez-Fernández, N. J., Pety, J., & Gueth, F. 2008, Single-dish observation and processing to produce the short-spacing information for a millimeter interferometer, Tech. rep., IRAM, Grenoble
  * Sault et al. (1995) Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 433.  https://arxiv.org/abs/astro-ph/0612759
  * Stanimirovic (2002) Stanimirovic, S. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 278, Single-Dish Radio Astronomy: Techniques and Applications, ed. S. Stanimirovic, D. Altschuler, P. Goldsmith, & C. Salter, 375–396.  https://arxiv.org/abs/astro-ph/0205329
  * Tang et al. (2021a) Tang, Y., Wang, Q. D., & Wilson, G. W. 2021a, MNRAS, 505, 2377, doi: 10.1093/mnras/staa3230
  * Tang et al. (2021b) Tang, Y., Wang, Q. D., Wilson, G. W., et al. 2021b, MNRAS, 505, 2392, doi: 10.1093/mnras/stab1191
  * The CASA Team et al. (2022) The CASA Team, Bean, B., Bhatnagar, S., Castro, S., et al. 2022, PASP, 134, 114501, doi: 10.1088/1538-3873/ac9642
  * van Moorsel et al. (1996) van Moorsel, G., Kemball, A., & Greisen, E. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2020) Wang, Y., Beuther, H., Rugel, M. R., et al. 2020, A&A, 634, A83, doi: 10.1051/0004-6361/201937095
  * Williams et al. (2018) Williams, G. M., Peretto, N., Avison, A., Duarte-Cabral, A., & Fuller, G. A. 2018, A&A, 613, A11, doi: 10.1051/0004-6361/201731587
  * Wilner & Welch (1994) Wilner, D. J., & Welch, W. J. 1994, ApJ, 427, 898, doi: 10.1086/174195