# WALLABY Pilot Survey: Public release of HI kinematic models for more than
100 galaxies from phase 1 of ASKAP pilot observations

N. Deg  Department of Physics, Engineering Physics, and Astronomy, Queen’s
University, Kingston ON K7L 3N6, Canada  K. Spekkens  Department of Physics
and Space Science, Royal Military College of Canada, P.O. Box 17000, Station
Forces Kingston ON K7K 7B4, Canada  T. Westmeier  International Centre for
Radio Astronomy Research (ICRAR), The University of Western Australia, 35
Stirling Highway, Crawley WA 6009, Australia  T.N. Reynolds  International
Centre for Radio Astronomy Research (ICRAR), The University of Western
Australia, 35 Stirling Highway, Crawley WA 6009, Australia  P. Venkataraman
Dunlap Institute of Astronomy and Astrophysics, University of Toronto, 50 St.
George Street, Toronto, ON, M5S 3H4, Canada  S. Goliath  NRC Herzberg
Astronomy and Astrophysics Research Centre, 5071 W. Saanich Rd., Victoria, BC,
V9E 2E7, Canada  A. X. Shen  CSIRO Space and Astronomy, PO Box 1130, Bentley
WA 6102, Australia  R. Halloran  Department of Physics, Engineering Physics,
and Astronomy, Queen’s University, Kingston ON K7L 3N6, Canada  A. Bosma  Aix
Marseille Univ, CNRS, CNES, LAM, Marseille  B. Catinella  International Centre
for Radio Astronomy Research (ICRAR), The University of Western Australia, 35
Stirling Highway, Crawley WA 6009, Australia  W.J.G. de Blok  Netherlands
Institute for Radio Astronomy (ASTRON), Oude Hoogeveensedijk 4, 7991 PD
Dwingeloo, The Netherlands  H. Dénes  Netherlands Institute for Radio
Astronomy (ASTRON), Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
E. M. Di Teodoro  Department of Physics & Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA  A. Elagali  Telethon Kids Institute, Perth
Children’s Hospital, Perth, Australia  B.-Q. For  International Centre for
Radio Astronomy Research (ICRAR), The University of Western Australia, 35
Stirling Highway, Crawley WA 6009, Australia  C. Howlett  School of
Mathematics and Physics, The University of Queensland, Brisbane QLD 4072,
Australia  G. I. G. Józsa  Max-Planck-Institut für Radioastronomie, Auf dem
Hügel 69, D-53121 Bonn, Germany  P. Kamphuis  Ruhr University Bochum, Faculty
of Physics and Astronomy, Astronomical Institute, 44780 Bochum, Germany  D.
Kleiner  INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5,
09047 Selargius, CA, Italy  B. Koribalski  ATNF, CSIRO Space and Astronomy, PO
Box 76, Epping NSW 1710, Australia  K. Lee-Waddell  International Centre for
Radio Astronomy Research (ICRAR), The University of Western Australia, 35
Stirling Highway, Crawley WA 6009, Australia  F. Lelli  INAF - Arcetri
Astrophysical Observatory, Largo Enrico Fermi 5, 50125, Florence Italy  X. Lin
School of Physics, Peking University, Beijing 100871, People’s Republic of
China  C. Murugeshan  CSIRO Space and Astronomy, PO Box 1130, Bentley WA 6102,
Australia  S. Oh  Department of Astronomy and Space Science, Sejong
University, 209, Neungdong-ro, Gwangjin-gu, Seoul, Republic of Korea  J. Rhee
International Centre for Radio Astronomy Research (ICRAR), The University of
Western Australia, 35 Stirling Highway, Crawley WA 6009, Australia  T. C.
Scott  Instituto de Astrofísica e Ciências do Espaço (IA), Rua das Estrelas,
4150-762 Porto, Portugal  L. Staveley-Smith  International Centre for Radio
Astronomy Research (ICRAR), The University of Western Australia, 35 Stirling
Highway, Crawley WA 6009, Australia  J.M. van der Hulst  Kapteyn Astronomical
Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The
Netherlands  L. Verdes-Montenegro  Instituto de Astrofísica de Andalucía (IAA-
CSIC), Glorieta de la Astronomía, 18008 Granada, Spain  J. Wang  Kavli
Institute for Astronomy and Astrophysics, Peking University, Beijing 100871,
China  O. I. Wong  CSIRO Space and Astronomy, PO Box 1130, Bentley WA 6102,
Australia

(23 June 2022; 24 Aug 2022; 02 Sept 2022)

###### Abstract

We present the Widefield ASKAP L-band Legacy All-sky Blind surveY (WALLABY)
Pilot Phase I Hi kinematic models.  This first data release consists of Hi
observations of three fields in the direction of the Hydra and Norma clusters,
and the NGC 4636 galaxy group. In this paper, we describe how we generate and
publicly release flat-disk tilted-ring kinematic models for 109/592 unique Hi
detections in these fields. The modelling method adopted here – which we call
the WALLABY Kinematic Analysis Proto-Pipeline (WKAPP) and for which the
corresponding scripts are also publicly available – consists of combining
results from the homogeneous application of the FAT and 3DBarolo algorithms to
the subset of 209 detections with sufficient resolution and $S/N$ in order to
generate optimized model parameters and uncertainties. The 109 models
presented here tend to be gas rich detections resolved by at least 3–4
synthesized beams across their major axes, but there is no obvious
environmental bias in the modelling. The data release described here is the
first step towards the derivation of similar products for thousands of
spatially-resolved WALLABY detections via a dedicated kinematic pipeline. Such
a large publicly available and homogeneously analyzed dataset will be a
powerful legacy product that that will enable a wide range of scientific
studies.

###### doi:

article in press

ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions
(ASTRO 3D), Australia Australian SKA Regional Centre (AusSRC) ARC Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia
Department of Astronomy, University of Cape Town, Private Bag X3, Rondebosch
7701, South Africa Kapteyn Astronomical Institute, University of Groningen, PO
Box 800, 9700 AV Groningen, The Netherlands Space Telescope Science Institute,
3700 San Martin Drive, Baltimore, MD 21218, USA ARC Centre of Excellence for
All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia Department of
Physics and Electronics, Rhodes University, P.O. Box 94, Makhanda, 6140, South
Africa School of Science, Western Sydney University, Locked Bag 1797, Penrith
NSW 2751, Australia CSIRO Space and Astronomy, PO Box 1130, Bentley WA 6102,
Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions
(ASTRO 3D), Australia Department of Physics and Astronomy, Sejong University,
209, Neungdong-ro, Gwangjin-gu, Seoul, Republic of Korea ARC Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia ARC
Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
Australia International Centre for Radio Astronomy Research (ICRAR), The
University of Western Australia, 35 Stirling Highway, Crawley WA 6009,
Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions
(ASTRO 3D), Australia : currently in press

##  1 Introduction

The Widefield ASKAP L-band Legacy All-sky Blind surveY (WALLABY; Koribalski et
al. 2020) is one of the key science projects for the Australian SKA Pathfinder
(ASKAP; Hotan et al. 2021) telescope. It is an extragalactic survey expected
to detect the atomic hydrogen (Hi) gas content of $\sim 210,000$ galaxies out
to redshift $z\sim 0.1$. It is expected that thousands of these sources will
have sufficient spatial resolution for kinematic modelling. WALLABY has
completed the first phase of its pilot observations, consisting of the three
fields towards the Hydra and Norma clusters, and the NGC 4636 group. Westmeier
et al. (2022, hereafter W22) is the release paper for the pilot data release 1
(PDR1) and this paper describes the PDR1 rotating disk kinematic models.

The generation of reliable kinematic models for as many resolved galaxies as
possible is a key science driver for WALLABY. Most, but not all, sources with
significant Hi reservoirs are rotationally supported as the gas generally
settles into rotating disks due to the conservation of angular momentum. As
such, we attempt to model all the sufficiently resolved PDR1 detections using
‘rotating disk’ kinematic models. Such models provide important measurements
for galaxies that are useful for exploring a variety of questions. For
instance, the rotation curves generated from such models are key to answering
questions related to the mass distribution within galaxies. Such questions
include whether or not disks are maximal (van Albada et al., 1985; van Albada
& Sancisi, 1986; Lelli et al., 2016; Starkman et al., 2018), and, with a large
enough sample size, probing the core-cusp problem (de Blok, 2010).
Additionally, studies of the Tully-Fisher (TF) relation (Tully & Fisher, 1977)
are significantly improved by measurements of $v_{\rm{flat}}$ (which is
derived from the outer rotation curves; see Verheijen 2001; Ponomareva et al.
2016; Lelli et al. 2019). Getting a statistically significant sample of
$v_{\rm{flat}}$ measurements will be valuable for galaxy population studies
involving the TF relation and the baryonic TF relation (McGaugh et al., 2000;
Oh et al., 2011).

Another key use of these rotation models is the calculation of the resolved
velocity function down to low masses (Lewis, 2019). This, coupled with mass
modelling will help to address questions related to the dark matter (DM) halo-
to-Hi velocity relation (Papastergis & Shankar, 2016). For larger galaxies,
kinematic models can help to constrain their spin, warps, and angular
momentum. The DM spin and Hi warps of a galaxy are often connected to
environmental processes (Battaner et al., 1990; Stevens et al., 2016; Lagos et
al., 2018).

These are just a few of the questions that require robust Hi kinematic models.
There are a variety of different methods for generating such kinematic models
from interferometric observations. The most common is tilted-ring modelling
(Rogstad et al., 1974), initially applied in 2D to velocity moments of well-
resolved Hi detections (e.g. Bosma, 1978; Begeman, 1987; van der Hulst et al.,
1992). More recently, algorithms have been developed to apply tilted-ring
models directly to 3D datacubes (e.g. Józsa et al., 2007; Davis et al., 2013;
Kamphuis et al., 2015; Di Teodoro & Fraternali, 2015; Bekiaris et al., 2016).
A key advantage of 3D techniques relative to their 2D counterparts is the
reliability with which they can be applied to marginally spatially-resolved Hi
detections; thus making them particularly useful for homogeneous application
to large numbers of detections from blind widefield surveys such as WALLABY.

This paper describes the construction of the PDR1 kinematic models (for the
subset of galaxies that were successfully modelled) as well as the public
release of the resulting data products. In Sec. 2 we briefly describe the
WALLABY detections, but a full description is provided in the data release
paper (W22). Sec. 3 describes tilted ring modelling in general, while Sec. 4
describes the specific approach taken for the PDR1 observations. Sec. 5
provides the overall results of the PDR1 kinematic models. Sec. 6 describes
the population of kinematically modelled galaxies, and Sec. 7 provides the
conclusions and discusses the future of kinematic modelling for full WALLABY.

##  2 WALLABY PDR1 Detections

The WALLABY pilot phase 1 observations targeted three $60\,\mathrm{deg}^{2}$
fields that cover cluster and group environments at differing distances $D$:
the Hydra cluster ($D\sim 60\,$Mpc, Jørgensen et al. 1996; Reynolds et al.
2021), the Norma cluster ($D\sim 70\,$Mpc, Mutabazi 2021), and the NGC 4636
group ($D\sim 15\,$Mpc, Tully et al. 2013). A full description of the
observations, the data reduction and the application of the SoFiA source
finding code (the HI Source Finding Application; Serra et al. 2015; Westmeier
et al. 2021) to generate the PDR1 sample of 592 unique Hi detections divided
between Hydra Team Release 1 (TR1), Hydra TR2, Norma TR1 and NGC 4636 TR1 is
reported in W22.

The PDR1 detection cubelets (cutouts from the mosaiced cubes around each
detected source), imaged with a Gaussian restoring beam with a full-width at
half maximum of $30^{\prime\prime}$ (hereafter “the beam”) in $18.5\,$kHz-wide
spectral channels ( $=3.9\,\textrm{km s}^{-1}$ at $z=0$), are the starting
point for the kinematic analysis presented here. W22 also detail the
limitations of the data given the pilot nature of the observations; we discuss
the potential impacts of the those limitations on the kinematic models in
Section 4.6, which we expect to be mild.

Figure 1 plots the angular size as a function of the integrated signal-to-
noise ($S/N$) for the detected sources. These two properties strongly
influence whether a galaxy’s Hi content can be reliably kinematically modelled
(see Section 5). We define the size in Fig. 1 as the SoFiA-returned major axis
diameter ell_maj of an ellipse fitted to the source Moment 0 map (Westmeier et
al., 2021). As in W22, we compute the integrated $S/N$ via:

$S/N_{obs}=\frac{S_{mask}}{\sigma_{rms}\sqrt{N_{mask}\Omega}}~{},$ (1)

where $S_{mask}$ and $N_{mask}$ are the total flux and number of cells in the
SoFiA detection mask respectively, $\Omega$ is the beam area, and
$\sigma_{rms}$ is the root-mean-square noise of the detection-free cells in
the corners of the SoFiA-returned source cubelet. In Fig. 1 and throughout, we
plot the Hydra TR2 values for sources with both Hydra TR1 and Hydra TR2
entries.

Figure 1: Size (as estimated by ell_maj) as a function of integrated $S/N$
(given by Eq. 1) of PDR1 detections. Sources in the Hydra, Norma, and NGC 4636
fields are indicated by circles, stars, and triangles respectively. Coloured
points represent all detections with $\texttt{ell\\_maj}>2\,$beams or
$\log(S/N_{obs})>1.25$, for which kinematic models were attempted: successful
models are shown in blue, and failed models are shown in red (see Section 4).
Moment maps for the sources corresponding to the points outlined in larger
open symbols are shown in Fig. 8.

Fig. 1 shows a clear correlation between angular size and $S/N$ among the
detections: as expected, sources with larger angular sizes have higher
integrated $S/N$. As also shown in W22, the majority of the detections are
only marginally spatially resolved, with values of ell_maj that span only a
few beams. Moreover, most detections have relatively low $S/N$. As such, our
modelling must be tailored for the marginally resolved, low $S/N$ regime. We
discuss considerations that drive the adopted modelling approach in Sec. 3,
and describe the resulting procedure in Sec. 4.

##  3 Kinematic Modelling Considerations

Given that many science goals for WALLABY are enabled by statistical samples
of resolved source properties (see Sec. 1), two core principles underpin our
kinematic modelling approach:

  1. 1.

Models should be automatically and homogeneously applied to all suitable
detections.

  2. 2.

Model parameters should have robust estimates of their uncertainties.

These principles drive key choices in the modelling undertaken. First, we do
not tailor kinematic models to individual detections; rather, we apply the
same models using the same technique to all sources that meet our selection
criteria. Second, since available algorithms do not return statistical
uncertainties on all parameters, we apply different code implementations of
the same underlying model to a given source in order to estimate the
uncertainties for the returned parameters.

Given these principles and the properties of the spatially-resolved PDR1
detections described in Sec. 2, we discuss here the considerations that drive
our kinematic modelling procedure. Sec. 3.1 introduces tilted-ring modelling
and describes the Fully Automated TiRiFiC (FAT, where TiRiFiC itself stands
for Tilted Ring Fitting Code; Kamphuis et al. 2015; Józsa et al. 2007) and the
3D-Based Analysis of Rotating Objects From Line Observations (3DBarolo, Di
Teodoro & Fraternali 2015) algorithms that we use to generate the PDR1 models.
Sec. 3.2 then explores differences in how these two codes model the same
underlying observation, which is used to build and hone the modelling
procedure adopted in Sec. 4.

###  3.1 Tilted-Ring Modelling

Tilted-ring modelling, first introduced by Rogstad et al. (1974), is a widely-
used technique for generating kinematic models of a galaxy’s Hi disk. In this
procedure, a model galaxy is constructed from a series of concentric rings,
each with intrinsic properties such as a centre, rotation speed, surface
density and thickness, as well as quantities that arise from the ring’s sky
projection, like inclination and position angle. While the precise set of
parameters included in the models varies by implementation, the goal is to
generate mock observations of the ring ensemble and optimize the ring
parameters so that they resemble the observations.

Tilted-ring models were initially developed for application to 2D velocity
fields derived from 3D Hi datacubes, with rotcur in the gipsy package being an
early and widely-used implementation (Begeman, 1987; van der Hulst et al.,
1992). A suite of more recent 2D algorithms that also characterize non-
circular flows or complex disk geometries have since been developed and
publicly released (e.g. reswri, Schoenmakers et al. 1997; Kinemetry, Krajnović
et al. 2006; DiskFit, Spekkens & Sellwood 2007; 2DBAT; Oh et al. 2018). 2D
algorithms are relatively efficient, and reliably recover the intrinsic and
projected ring properties when the Hi disk is at intermediate inclination
(generally in the range $40^{\circ}-75^{\circ}$) and spatially resolved by
$\sim 8-10$ beams across the major axis (e.g. Bosma 1978, Kamphuis et al.
2015).

More recent tilted-ring codes have generalized the approach for application
directly to the 3D datacubes themselves (e.g. TiRiFiC, Józsa et al. 2007;
KinMS, Davis et al. 2013; 3DBarolo, Di Teodoro & Fraternali 2015; FAT,
Kamphuis et al. 2015; GBKFit, Bekiaris et al. 2016). 3D techniques have two
main advantages relative to 2D ones: first, they allow for more complicated
morphological and kinematic models to be applied to deep, high-resolution data
(e.g. Józsa et al., 2009; Khoperskov et al., 2014; Di Teodoro & Peek, 2021;
Józsa et al., 2021); and second, they can be robustly applied at lower spatial
resolutions and across a wider range of disk geometries than in 2D (e.g.
Kamphuis et al., 2015; Di Teodoro & Fraternali, 2015; Lewis, 2019; Jones et
al., 2021). Given the size distribution of sources implied by Fig. 1, it is
this latter property that makes 3D techniques most suitable for homogeneous
modelling of PDR1 detections.

We work with FAT and 3DBarolo, two publicly-available codes designed to
automatically apply 3D tilted-ring models to samples of Hi datacubes. Below,
we describe the salient properties of both algorithms in the context of the
PDR1 kinematic analysis.

####  3.1.1 FAT

FAT111https://github.com/PeterKamphuis/FAT (Kamphuis et al., 2015) automates
the application of TiRiFiC222https://gigjozsa.github.io/tirific/ (Józsa et
al., 2007), one of the first and most well-developed 3D tilted-ring codes.
TiRiFiC constructs models by populating rings with tracer particles,
projecting them into a 3D datacube, and convolving the result with a 3D kernel
to match the spatial and spectral resolution of the data to which the model is
compared. The model is then optimized by computing the channel-by-channel
goodness of fit using an implementation of the ‘golden section’ search
algorithm (Press et al., 1992).

The basic approach implemented in FAT is to automatically initialize TiRiFiC
using parameters determined from applying the SoFiA source finder to the input
datacube, and then to iteratively apply TiRiFiC, usually with increasing
complexity, until a satisfactory fit is achieved. FAT begins with a flat-disk
model in which the ring geometries are independent of galactocentric radius
$R$, and has the functionality to explore radial variations in subsequent
iterations or to fit flat-disk models. By design, FAT estimates an
axisymmetric rotation curve but computes the surface brightness profile on the
approaching and receding sides of the disk separately. Once a satisfactory fit
of the parameters is found, radial variations are smoothed by a polynomial fit
to avoid artificial fluctuations from the TiRiFiC fitting algorithm, with
differences between smoothed and unsmoothed curves returned as uncertainties
for some parameters.

In a series of validation tests on real and simulated data, Kamphuis et al.
(2015) show that FAT can reliably recover both the geometries and kinematics
of Hi disks with inclinations ranging from $20^{\circ}-90^{\circ}$ that are
spatially resolved by at least 8 beams across their major axes, while
extensive tests by Lewis (2019) imply that FAT can recover inclinations and
rotation curves for flat, axisymmetric disks resolved by as few as 3.5 beams
across their major axis diameters, $D_{HI}$, in the inclination range
$35^{\circ}-80^{\circ}$. We note that $D_{HI}$ differs from the SoFiA-returned
ell_maj shown in Fig. 1 (see Sec. 5). We work with FAT version 2.01.

####  3.1.2 3DBarolo

3DBarolo333https://editeodoro.github.io/Bbarolo/ (Di Teodoro & Fraternali,
2015) is a tilted-ring code that has been extensively used to apply 3D models
to Hi datasets in different resolution and $S/N$ regimes. Many elements of the
3DBarolo implementation are similar to those described for FAT above; below,
we highlight differences that are relevant to the PDR1 kinematic analysis.

Key 3DBarolo features that differ from FAT are parameter initialization, model
optimization and flux normalization. 3DBarolo can use a built-in source finder
based on DUCHAMP (Whiting, 2012) to initialize the models, or the user can
specify initial parameter estimates directly. Once the source(s) are found,
the model is optimized on a ring-by-ring basis using the Nelder-Mead algorithm
(Nelder & Mead, 1965), where beam effects are mimicked by convolving each
velocity channel with a 2D kernel. 3DBarolo can compute radial variations of
the geometric parameters using a number of different strategies such as
polynomial fits or Bezier interpolation, or can return median values if a
flat-disk model is specified. The model cube flux is normalized in 2D using
the observed moment 0 map, either on a pixel-by-pixel basis or on an azimuthal
ring basis. This approach increases the efficiency of the 3DBarolo
optimization relative to the channel-by-channel method adopted in TiRiFiC, but
limits the range of disk inclinations and surface density distributions that
can be robustly recovered. 3DBarolo implements a Monte Carlo approach to
estimate uncertainties for some parameters, where models are varied around the
best fit until the residuals increase by some factor (typically 5%).

In a series of validation tests on real data, Di Teodoro & Fraternali (2015)
show that 3DBarolo can efficiently recover the geometries and kinematics of
well-resolved and moderately-resolved Hi disks at intermediate inclinations
from the THINGS (Walter et al., 2008) and WHISP (Swaters et al., 2002)
surveys, respectively, while tests on both real data and galaxy mocks imply
that 3DBarolo can recover rotation curves and velocity dispersion profiles in
systems resolved by as few as 2 beams along the major axis when the
inclination is fixed and in the range $45^{\circ}-75^{\circ}$. We work with
3DBarolo version 1.6.

###  3.2 Application to PDR1 Detections

The key differences between FAT and 3DBarolo described above imply that the
same fitting options applied to the same dataset by each code may yield
different optimizations. These differences are typically small for spatially
well-resolved, high S/N detections, but may be significant in the low-
resolution, low-S/N regime in which most PDR1 sources lie (see Fig. 1;
Kamphuis et al. 2015; Di Teodoro & Fraternali 2015). Early in the pilot survey
phase, we therefore explored a suite of different FAT and 3DBarolo model
applications to over a dozen Hydra TR1 detections in order to develop the
technique we ultimately adopted. Because the sizes and S/N of most PDR1
sources pose challenges to tilted-ring modelling even with 3D applications, we
restricted the analysis to simple, flat-disk models where the disk geometry
does not vary with $R$.

Figure 2: Comparison between different flat-disk model outputs applied with
FAT and 3DBarolo to WALLABY J103915-301757 (ell_maj$=3.9~{}\textrm{beams}$,
$\log(S/N_{obs})=1.5$). The top row shows the moment 0 and moment 1 maps of
this source, with the cyan circle indicating the size of the beam. The panels
below show plots of the rotation curve (A), surface density profile (B),
inclination (C), position angle (D), kinematic centre relative to the PDR1
source centroid (E and F), systemic velocity (G) and velocity dispersion
profile (H) as a function of galactocentric radius $R$ for the flat-disk
models given in the legend (see text for details), evaluated at the locations
given by the points. The dashed black lines in panels C, D, E, F and G
indicate the PDR1 source parameters for those quantities from W22. The error
bars on some profiles in some panels are the final uncertainties returned by
either FAT or 3DBarolo for that model application.

This experimentation revealed that among possible flat-disk modelling choices
in FAT and 3DBarolo, a) the $S/N$ of the detected emission in each channel and
b) the model parameter initialization can both strongly influence the
optimizations returned in the PDR1 regime, with variations in other algorithm
switches having comparatively minor effects. The $S/N$ of the emission per
channel can impact the reliability of the built-in source finder which
initializes parameters, hence the importance of that modelling choice. The
parameter initialization, in turn, can impact the model outputs because
optimization schemes such as Golden Section (as in FAT) and Nelder-Mead (as in
3DBarolo) require robust initial guesses to converge in the complex, multi-
dimensional parameter spaces characteristic of 3D tilted-ring models (e.g.
Bekiaris et al. 2016).

We illustrate these trends in the PDR1 regime in Fig. 2, which shows the
output parameters for several flat-disk models applied to WALLABY
J103915-301757 (ell_maj$=3.9~{}\textrm{beams}$, $\log(S/N_{obs})=1.5$) with
FAT and 3DBarolo. This example is just one of the dozen galaxies tested with
different modelling options and illustrates well how different choices affect
the resulting models. The main differences between the different fitting
attempts are the spectral resolution of the cubelet to which the model is
applied (either the full-resolution cubelet, or a 3-channel Hanning-smoothed
cubelet), and the choice of geometric parameter initialization (either
initialized automatically by the code or initialized by the user to the PDR1
source values from W22):

  * •

Barolo_full_auto: 3DBarolo applied to the full-resolution cubelet, with
automated parameter initialization;

  * •

Barolo_smooth_auto: 3DBarolo applied to the spectrally-smoothed cubelet, with
automated parameter initialization;

  * •

Barolo_smooth_source: 3DBarolo applied to the spectrally-smoothed cubelet,
with geometric parameters initialized to the PDR1 source values;

  * •

FAT_full_auto: FAT applied to the full-resolution cubelet, with automated
parameter initialization;

  * •

FAT_smooth_auto: FAT applied to the spectrally-smoothed cubelet, with
automated parameter initialization;

  * •

Barolo_auto_smooth_vdisp: 3DBarolo applied to the spectrally-smoothed cubelet,
allowing the velocity dispersion to vary.

We note that since FAT does not allow the user to initialize parameters, there
no such model with this option in Fig. 2. We note also that since many PDR1
detections have no optical counterparts (particularly in Norma TR1, which is
close to the Galactic Plane), we do not attempt to initialize geometric
parameters with photometric values. We note that the Barolo_auto_smooth_vdisp
shown in Fig. 2 involves a different fitting mode than the other models, and
is discussed further below.

Comparing the optimized rotation curves (Fig. 2A), surface density profiles
(Fig. 2B) and disk geometries (Fig. 2C–G) across models for WALLABY
J103915-301757, the 3DBarolo application to the full-resolution cubelet
(Barolo_full_auto) differs markedly from the other outputs: its radial extent
is much smaller than that of the source plotted in the top row, and the disk
geometry (most notably the inclination) is discrepant with the source
morphology. This model failure stems from an incorrect source identification
and parameter initialization by the 3DBarolo source finder due to the low S/N
of the emission in each channel of the full-resolution cubelet.

Regardless of the model, the position angle and geometric center (Fig. 2D–F)
are recovered well. Given the pixel size is $\sim 6$′′, the kinematic center
is recovered within less than 2 pixels. The successful measurement of these
three geometric parameters is typical for kinematic modelling as they tend to
have fewer degeneracies with other parameters than the inclination or systemic
velocity. If there are large differences between the kinematic center and
position angle for various fits, it indicates a failure in one or more of
those fits. The systemic velocity itself (Fig. 2G) shows a larger variation,
but for the smoothed cubes, the differences are only of the order of 1-2
channels. The greatest outlier is the Barolo_smooth_source model, which also
is an outlier in terms of the spatial center.

We find that in general, models applied smoothed to cubelets are more stable
and converge faster than those applied to the full-resolution cubelets, with
little difference between the optimized values when both models succeed (e.g.
FAT_full_auto and FAT_smooth_auto in Fig. 2). This is not unexpected since the
modelled PDR1 detections are spectrally well-resolved in both the full-
resolution and smoothed cubelets, while the per-channel S/N of the emission is
$\sim$50% higher in the latter. We therefore elect to kinematically model PDR1
cubelets that have been spectrally smoothed by a 3-channel Hanning window.

Another trend that emerges among successful models in Fig. 2 is that the
outputs from different models applied using the same code (e.g. FAT_full_auto
vs. FAT_smooth_auto) are typically more similar than those from the same model
applied by different codes (e.g. Barolo_full_smooth vs. FAT_full_smooth), with
differences that can well exceed the returned uncertainties. We find this to
be generally the case for the rotation curve and surface brightness profiles
at radii within $\sim$1 beam of the kinematic centre (Fig. 2A and B): the FAT-
returned rotation curves tend to rise more steeply and surface brightness
distributions tend to exhibit greater central depressions than the 3DBarolo
counterparts, particularly for relatively high inclination and/or poorly
spatially resolved sources. These discrepancies may stem in part from the
different radial grid definitions adopted by the codes (3DBarolo returns model
values at the ring mid-points, whereas FAT uses the ring edges), but
differences in optimization methodology (see Sec. 3) likely play a stronger
role. Regardless of the cause, the key point is that the differences between
successful FAT and 3DBarolo fits are typically larger than the reported
uncertainties. Our PDR1 modelling approach therefore adopts an average of the
models returned by each code as the optimal model, and differences between
them as a measure of uncertainty.

Figure 3: Velocity dispersion profiles from models identical to
Barolo_auto_smooth_vdisp in Fig. 2, applied to all 36 PDR1 sources with
$2\leq\texttt{ell\\_maj}\leq 4$ and $1.25\leq\log(S/N_{obs})\leq 1.5$. The
profiles are coloured according to the model disk inclination.

The dashed horizontal lines in Fig. 2 plot the PDR1 source parameters that
best approximate the geometric parameters returned by the kinematic modelling:
$\cos^{-1}(\texttt{ell\\_min}/\texttt{ell\\_maj})$, kin_pa, ra, dec, and freq
(see table 3 of W22) converted to the appropriate units are shown in Fig. 2C-G
respectively. Model Barolo_smooth_source initializes the model geometric
parameters to these values; for most successful models the output parameters
are nearly identical to the inputs, and different from the outputs from runs
in which the geometric parameters are initialized automatically in either the
Barolo_smooth_auto model or in the FAT_smooth_auto model. We find this to be
generally true for the PDR1 models, and speculate that the tilted-ring
parameter space is sufficiently complex that the 3DBarolo optimizer is
unlikely to find a step that improves the goodness of fit during the runs as
configured. Since the PDR1 parameters only approximate the kinematic model
parameters in the first place, we elect to use automatic source initialization
in the kinematic analysis.

We now discuss model Barolo_auto_smooth_vdisp in Fig. 2, which is identical to
Barolo_auto_smooth except that the disk velocity dispersion is allowed to vary
with $R$. Save for small differences between the rotation curves at $R\sim
50$′′ and the very different velocity dispersion profiles, the returned
parameters are almost identical between the two models, with corresponding
lines overlapping completely in Fig. 2B–G.

We find the independence of the model outputs on the velocity dispersion
switch, as well as the large radial variations in velocity dispersion when it
is allowed to vary, to be general trends in our PDR1 models. The velocity
dispersion variations are further illustrated in Fig. 3, where models
identical to Barolo_auto_smooth_vdisp were applied to the 36 PDR1 sources in
Hydra TR2, Norma TR1, and NGC 4636 TR1 with $2\leq\texttt{ell\\_maj}\leq 4$
and $1.25\leq\log(S/N_{obs})\leq 1.5$. This figure shows that, independent of
disk inclination, there is a general trend of decreasing velocity dispersion
with increasing $R$, but also variations across profiles and between them that
well exceed the relatively tight range $8\,\textrm{km s}^{-1}\leq\sigma\leq
12\,\textrm{km s}^{-1}$ measured from high-resolution Hi maps (Tamburro et
al., 2009). This is perhaps not surprising given the $\sim\,12\textrm{km
s}^{-1}$ resolution of the Hanning-smoothed cubelets that we model. We
therefore keep the velocity dispersion fixed to $\sigma=10\,\textrm{km
s}^{-1}$ (intermediate to the range of values typically measured) in the PDR1
kinematic models.

##  4 The WALLABY Kinematic Analysis Proto-Pipeline

Having explored some key considerations for kinematic modelling of PDR1
detections in Sec. 3, we now describe the approach adopted to derive flat-disk
tilted-ring models by applying FAT and 3DBarolo to pre-processed PDR1 cubelets
and averaging successful fits. We call the procedure the WALLABY Kinematic
Analysis Proto-Pipeline (WKAPP), and the full set of driving scripts is
available from its distribution page444https://github.com/CIRADA-Tools/WKAPP.
Fig. 4 summarizes the modelling steps:

  1. 1.

Select detections on which kinematic modelling is attempted and pre-process
their PDR1 cubelets;

  2. 2.

Apply FAT and 3DBarolo models to pre-processed cubelets;

  3. 3.

Generate optimized models and uncertainties by averaging successful model
fits;

  4. 4.

Compute surface density profiles from PDR1 source Moment 0 maps using
optimized model geometries.

Figure 4: A schematic of the WKAPP modelling process. The blue parallelograms
indicate data products, the green diamonds indicate decision points, and
yellow boxes indicate automated code.

We describe each of these steps in Sec. 4.1–4.4, the model outputs in Sec.
4.5, and some limitations of the current approach in Sec. 4.6.

###  4.1 Detection selection and cubelet pre-processing

The first step of WKAPP is to select a set of PDR1 detections on which
kinematic modelling is attempted. Validation tests on FAT and 3DBarolo suggest
that the algorithms can be successfully applied to Hi disks with diameters
$D_{HI}$ that are resolved by as few as 2-3 beams depending in the $S/N$ (Di
Teodoro & Fraternali, 2015; Lewis, 2019). We use the PDR1 size measure ell_maj
in our selection, which is typically a factor of two smaller than
$D_{\text{H\sc{i}}}$ in our successful models (see Sec. 6 for a comparison of
$D_{\text{H\sc{i}}}$ to ell_maj). Therefore we attempt to model all detections
with $\texttt{ell\\_maj}\geq 2$ beams. Because ell_maj is not a direct measure
of disk size, we also attempt to model all detections with
$\log(S/N_{obs})\geq 1.25$, even if they are below the size threshold. These
selection cuts result in 209 unique PDR1 detections that we attempt to model,
shown by the red and blue points in Fig. 1.

Next, the PDR1 cubelets selected for modelling are pre-processed in two steps.
First, the spectral axis of the cubelets is converted to velocity units from
the frequency units provided in the PDR1 data release. Second, the cubelets
are Hanning-smoothed by three spectral channels (to a resolution of
$11.7\,\textrm{km s}^{-1}$ at $z=0$) using 3DBarolo. As discussed in Sec. 3.2,
the main driver of this choice is an increase in model stability for
essentially the same model fit quality. It also decreases the FAT and 3DBarolo
run time since there are fewer spectral channels.

###  4.2 Application of FAT and 3DBarolo models

For each of the PDR1 detections selected as described above, we automatically
apply flat-disk tilted-ring models to the pre-processed cubelets using
3DBarolo and FAT. As discussed in Sec. 3.2, we allow each code to
automatically initialize all parameters, and we fix the velocity dispersion to
$10\,\textrm{km s}^{-1}$ in the models.

For 3DBarolo, the ring widths are set to 2 rings/beam and we use the SEARCH
source-finding method and the azimuthal normalization method (in order to be
as similar to FAT as possible). FAT does not allow the ring size to be
specified, but it generally fits objects with 2 rings/beam as well. For
completeness, both the input and results files from the 3DBarolo and FAT
applications to each successfully-modelled detection are distributed with the
data release (see Sec. 5).

###  4.3 Fit Inspection and Optimal Model Geometry and Rotation Curve

Only a subset of selected sources are successfully modelled using either FAT
or 3DBarolo: some show complicated structures that are not well-described by
flat disks, some are actually pairs of galaxies, and many have too low of a
resolution or $S/N$ to be successfully modelled (see Sec. 5). The results of
the 3DBarolo and FAT fits for each source are therefore visually examined to
determine their success. If either code fails to produce a model, if the final
model for either code is non-physical (for example, Barolo_full_auto in Fig.
2; see Sec. 3.2), or if the models returned differ strongly (for example,
$\delta\sin(i)>0.2$ between the FAT and 3DBarolo results), then the source is
discarded from the kinematic modelling sample (see Fig. 4).

If both the 3DBarolo and FAT fits are successful, then the two fits are
averaged together in three distinct steps to generate an optimal kinematic
model. The first is to directly average the geometric parameters (center,
$V_{sys}$, inclination, and position angle), with the uncertainty set to half
the difference between them. Table 1 shows an example for WALLABY
J163924-565221 (ell_maj$=4.5$ beams and $\log(S/N_{\mathrm{obs}})=1.53$), a
relatively large and high $S/N$ PDR1 detection (see Fig. 1).

The averaged model geometry is then used to calculate the optimal rotation
curve from the outputs of the FAT and 3DBarolo models. Since these models
typically have different radial extents and are evaluated at different values
of $R$, a final radial grid must be constructed. The final grid has two points
per beam; the innermost point is set to the larger of the two smallest model
$R$, which also defines the grid values. To optimize the radial extent of the
models, the outermost rotation curve point is the largest $R$ on the grid at
which one model is interpolated and the other is extrapolated by no more than
half a beam. Figure 5 shows an example of the radial grid definition for the
rotation curve of WALLABY J163924-565221.

With the final geometric parameters calculated and the radial grid set, the
3DBarolo and FAT rotation curves are adjusted to the final inclination and
interpolated onto the final grid using a spline fit. As with the geometric
parameters, the uncertainty on each rotation curve point is set to half the
difference between the two interpolated curves at that $R$. We also propagate
the effect of the inclination uncertainty to the rotation curve, providing a
separate value for this source of error. It is recommended that these two
uncertainties be added in quadrature when working with the model rotation
curves. An example of the optimal rotation curve calculation is given in Fig.
5 for WALLABY J163924-565221.

We note that the FAT and 3DBarolo model rotation curves are generated with a
degree of internal consistency. But it is not guaranteed that our optimized
rotation curves will have similar levels of self-consistency as they are
generated by averaging the interpolatated, inclination corrected FAT and
3DBarolo outputs. However, the visual inspection of the different fits as well
as a final examination of the optimized model help to avoid any
inconsistencies, and in practice the best fitting disk centres and position
angles are typically very similar (see Sec. 3.2). We therefore judge that the
successful models have rotation curves that are internally consistent with the
rest of the model parameters.

###  4.4 Surface Density Profile Computation

In the final WKAPP step, the surface density profile is calculated from
ellipse fits to the PDR1 detection Moment 0 map and the average geometry. In
other words, the surface density profile is derived separately from the FAT
and 3DBarolo estimates of this parameter, but using the same disk geometry as
in the optimized model. This approach is similar to the 3DBarolo procedure for
calculating surface densities but differs strongly from the FAT approach (see
Sec. 3.1), where they are constrained directly from the cube; since the FAT
approach has not been vetted in the resolution and $S/N$ regime of the PDR1
detections, we use ellipse fits for this first public data release with the
goal of using cube fits in future ones (see Sec. 4.6).

The optimized surface density profile is computed using the same radial grid
values as the rotation curve, but with the extent determined by the PDR1 mask
width along the kinematic major axis of the Moment 0 map. In practice, this
implies that the surface density profile of a given model typically extends to
larger $R$ than its rotation curve; this choice implies that the majority of
the surface density profiles extend to the characteristic density
$\Sigma=1\,\mathrm{M_{\odot}\,pc^{-2}}$ at which disk radii are typically
defined (e.g. Wang et al., 2016), although this requires extrapolating the
disk geometry beyond the region used in the model fits. We adopt the standard
error on the mean as the uncertainty in the measured profiles, that is the
standard deviation of the pixels in each ring divided by the square root of
the number of beams in that ring.

In addition to providing the surface density profile directly measured from
the ellipse fits, we also provide a version to which a standard $\cos(i)$
correction has been applied to deproject the profiles to a face-on view. We
caution that this correction can strongly under-estimate the inner surface
densities of marginally-resolved HI disks, as is the case for many PDR1
detections (see Fig. 1). In addition, we do not attempt to correct the outer
surface density profiles for beam smearing effects. We discuss both effects in
Sec. 4.6, and recommend that their impact be considered when using the
corrected surface density profiles.

###  4.5 Model Outputs

Parameter | Units | FAT | 3DBarolo | Model | Uncertainty
---|---|---|---|---|---
$X$ | px | 40.5 | 39.5 | 40.0 | 0.5
$Y$ | px | 33.5 | 33.8 | 33.7 | 0.1
$V_{sys}$ | $\textrm{km s}^{-1}$ | 1466.9 | 1465.7 | 1466.3 | 0.6
Inc | deg | 49.4 | 37.0 | 43.2 | 6.2
PA | deg | 249.9 | 245.6 | 247.8 | 2.2
Table 1: Geometric parameter averaging for the WKAPP model to WALLABY
J163924-565221 (ell_maj$=4.5$ beams and $\log(S/N_{\mathrm{obs}})=1.53$). The
FAT and 3DBarolo columns show the results from the fits to the galaxy using
the respective codes. the Model and Uncertainty columns show the average
geometric parameters and rounded uncertainties adopted.

Every PDR1 source that is successfully modelled by WKAPP is characterized by a
set of model parameters as listed in Table 2. The geometric parameters and
associated uncertainties are single values, while the rotation curve and
surface density profiles and their associated uncertainties are arrays.

We note that among the geometric parameters provided for each model are PAs in
pixel coordinates (PA_model) and in global equatorial coordinates
(PA_model_g). For most PDR1 detections, there is a small but non-zero
rotational offset between those two coordinate systems that is defined by the
PDR1 cubelet header. This results in a small but systematic difference between
PA_model and PA_model_g (typically less than 2 degrees).

As described in Sec. 4.3, we provide estimates of uncertainty from two
different sources for each rotation curve: the first (e_Vrot_model) arises
from the FAT and 3DBarolo averaging process, and the second (e_Vrot_model_inc)
is the contribution to the uncertainty on the rotation curve obtained by
propagating the uncertainty on the inclination (e_Inc_model). Fig. 5 provides
an example of these two sources of uncertainty for WALLABY J163924-565221. We
recommend adding these sources in quadrature when using the rotation curve.

We also provide estimates of the uncertainty on the surface density profile
from two sources. The first (e_SD_model) is the standard error of the pixels
in each ring, which we recommend be adopted as the uncertainty in the
projected surface density profile (SD_model). We also provide an estimate of
the statistical uncertainty for the profile deprojection (e_SD_FO_model_inc)
obtained by propagating the uncertainty on the inclination (e_Inc_model). We
recommend adding these sources in quadrature when the deprojected surface
density profile (SD_FO_model) is used for scientific analysis, but caution
that for many PDR1 sources systematic errors in the standard $\cos(i)$
correction dominate. We discuss this further in Sec. 4.6 below.

We also provide a quality flag for each model:

  * •

QFlag_model = 0: No obvious issues.

  * •

QFlag_model = 1: $\texttt{Inc\\_model}\leq 20^{\circ}$ or
$\texttt{Inc\\_model}\geq 80^{\circ}$.

  * •

QFlag_model = 2: $\texttt{e\\_Vsys\\_model}\geq 15\,\textrm{km s}^{-1}$.

  * •

QFlag_model = 3: Both conditions 1 and 2 are met.

Figure 5: Example showing how the 3DBarolo and FAT rotation curve fits are combined into a single average model for WALLABY J163924-565221, where the geometric parameters are given in Table 1. The black line shows the optimal rotation curve, while the red and blue lines show the outputs from the automated 3DBarolo and FAT models fits respectively. The solid error bars show the uncertainty from averaging the interpolated inclination-adjusted rotation curves, while the dashed error bars show the uncertainty in the inclination propagated into the rotation curve. In this example, the latter uncertainty is much larger than the former for most points. Name |  Type |  Unit |  Description
---|---|---|---
X_model |  double |  px |  x-coordinate of the kinematic center†
e_X_model |  double |  px |  Uncertainty in X_model†
Y_model |  double |  px |  y-coordinate of the kinematic center.†
e_Y_model |  double |  px |  Uncertainty in Y_model†
RA_model |  double |  deg |  Right ascension (J2000) of the kinematic center
e_RA_model |  double |  deg |  Uncertainty in RA_model†
DEC_model |  double |  deg |  Declination (J2000) of the kinematic center
e_DEC_model |  double |  deg |  Uncertainty in DEC_model†
Vsys_model |  double |  $\textrm{km s}^{-1}$ |  Heliocentric systemic velocity
e_Vsys_model |  double |  $\textrm{km s}^{-1}$ |  Uncertainty in Vsys_model
Inc_model |  double |  deg |  Inclination
e_Inc_model |  double |  deg |  Uncertainty in Inc_model
PA_model |  double |  deg |  Position angle in pixel coordinates (counterclockwise from x=0)†
e_PA_model |  double |  deg |  Uncertainty in PA_model†
PA_model_g |  double |  deg |  Position angle in equatorial coordinates (East of North)
e_PA_model_g |  double |  deg |  Uncertainty in PA_model_g
Rad |  double array array |  arcsec |  Radial grid for Vrot_model
Vrot_model |  double array |  $\textrm{km s}^{-1}$ |  Rotation curve
e_Vrot_model |  double array |  $\textrm{km s}^{-1}$ |  Uncertainty in Vrot_model from the averaging process
e_Vrot_model_inc |  double array |  $\textrm{km s}^{-1}$ |  Uncertainty in Vrot_model due to e_Inc_model
Rad_SD |  double array |  arcsec |  Radial grid for SD_model and SD_FO_model
SD_model |  double array |  M⊙ pc-2 |  Projected surface density profile
e_SD_model |  double array |  M⊙ pc-2 |  Uncertainty in SD_model
SD_FO_model |  double array |  M⊙ pc-2 |  Deprojected surface density profile using a $\cos(\texttt{Inc\\_model})$ correction
e_SD_FO_model_inc |  double array |  M⊙ pc-2 |  The uncertainty in SD_FO_model due to e_Inc_model
QFlag_model |  integer |  |  Kinematic model quality flag
${}^{\dagger}\,$In pixel coodinates relative to the preprocessed cubelet,
which starts from the point (1,1).
Table 2: WKAPP model parameters.

We flag models with inclinations below $20^{\circ}$ and above $80^{\circ}$
(QFlag_model = 1) because, although we judge these fits to be successful,
these inclinations lie in the range where neither FAT nor 3DBarolo have been
vetted (Di Teodoro & Fraternali, 2015; Kamphuis et al., 2015; Lewis, 2019). We
similarly have judged fits with $\texttt{e\\_Vsys\\_model}\geq 15\,\textrm{km
s}^{-1}$ (QFlag_model = 2) to be successful, but they are strong outliers in
the distribution of this value (see Sec. 5) which may indicate a subtle
failure that is not obvious through visual inspection. $\sim 12\%$ (16/124) of
all modelled sources have been flagged: QFlag_model = 2 and QFlag_model = 3
have each been assigned once, with the remaining 14 sources having been
assigned QFlag_model = 1.

File suffix |  Type |  Description
---|---|---
_AvgMod.txt |  ascii file |  Model parameters
_DiagnosticPlot.png |  PNG file |  Model summary plot
_ProcData.fits |  FITS cube |  Pre-processed cubelet
_ModCube.fits |  FITS cube |  Model realization with pre-processed cubelet properties
_DiffCube.fits |  FITS cube |  Data - model cube with pre-processed cubelet properties
_ModRotCurve.fits |  FITS binary table |  Table containing the model rotation curve parameters
_ModSurfDens.fits |  FITS binary table |  Table containing the model surface density parameters
_ModGeo.fits |  FITS binary table |  Table containing the model geometry parameters
_FullResProcData.fits |  FITS cube |  Full spectral resolution cubelet with velocity units
_FullResModelCube.fits |  FITS cube |  Model realization with full resolution cubelet properties
_FATInput.txt |  ascii file |  The input file of the FAT run
_FATMod.txt |  ascii file |  The results of the FAT run
_BaroloInput.txt |  ascii file |  The input file of the 3DBarolo run
_BaroloMod.txt |  ascii file |  The geometry and rotation curve results of the 3DBarolo run
_BaroloSurfDens.txt |  ascii file |  The surface density results of the 3DBarolo run
Table 3: WKAPP data products available for each successfullly modelled PDR1
source.

In addition to the catalog of model parameters for all kinematically modelled
PDR1 sources, WKAPP also produces a number of data products for each source.
They are listed in Table 3. Several products serve to group model parameters
for individual sources into distinct files for ease of access: files with
suffix _AvgMod.txt contain all model parameters for the PDR1 source in the
prefix, while those with suffixes _ModRotCurve.fits, _ModSurfDens.fits, and
_ModGeo.fits store the rotation curve, surface density profile, and geometric
parameters as FITS binary tables respectively. Additionally, text files with
FAT or Barolo in the suffix provide the input and output files from the
automated FAT and 3DBarolo applications described in Sec. 4.2.

Several data and model cubelets are also provided as data products. The model
cubelets are realizations of the optimized models using the stand-alone
tilted-ring model generator MCGSuite code555https://github.com/CIRADA-
Tools/MCGSuite (Lewis 2019, Spekkens et al. in prep). The pre-processed PDR1
cubelets to which FAT and 3DBarolo are applied (see Sec. 4.2) are in files
with suffix _ProcData.fits. Realizations of the optimized models in cubelets
with the same properties as the pre-processed data cubelets as well as data –
model cubelets are in files with suffixes _ModeCube.fits and _DiffCube.fits,
respectively. For completeness, we also provide PDR1 cubelets at full spectral
resolution with the frequency axis in velocity units (suffix
_FullResProcData.fits), as well as model realizations with the properties of
those cubelets (suffix _FullResModelCube.fits).

Finally, a summary plot is provided for each modelled source as a PNG file
with suffix _DiagnosticPlot.png. As an example, Figure 6 shows the summary
plot for WALLABY J100426-282638 (ell_maj$=5.0$ beams and
$\log(S/N_{\mathrm{obs}})=1.83$), one of the largest and highest $S/N$ PDR1
detections (see Fig. 1). While the model cubelets and summary plots may be
useful for a variety of scientific applications, it is important to note that
the key data products are the WKAPP model parameters and uncertainties from
which the other data products are derived.

Figure 6:  Sample WKAPP model summary plot for WALLABY J100426-282638
(ell_maj$=5.0$ beams and $\log(S/N_{\mathrm{obs}})=1.83$). Similar summary
plots are included in the data release for each modelled PDR1 detection. The
upper left and right panels show the rotation curve and surface density
profile of the optimized model. The middle left panel shows the PDR1 Moment 0
map and the location of the model center marked with a black X. The middle
right panel shows the PDR1 Moment 1 map, along with the model velocity
contours (constructed from the MCGSuite cube realization), and the direction
of the model position angle marked by a black arrow. The bottom panels show
the major and minor axis position-velocity (PV) diagrams (left and right
panels respectively) along with the corresponding model PV diagrams (magenta
lines). The model contours are at 3 and 5$\sigma$ of the PV diagram noise. The
major axis PV diagram also shows the projected rotation profile
($\,=\,$Vrot_model$\times\sin[$Inc_model$]$).

The WKAPP data products are accessible via both CSIRO ASKAP Science Data
Archive (CASDA)666https://research.csiro.au/casda/ and the Canadian Astronomy
Data Centre (CADC)777https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/. A full
description of the data access options can be found in both W22 and through
the data access portal888https://wallaby-survey.org/data/.

###  4.6 Model Limitations

The procedure adopted here for producing kinematic models of the WALLABY PDR1
galaxies is a reasonable first effort. However, it is important to note that
there are limitations to both elements of the approach adopted as well as to
the underlying data; we discuss them below, in what we judge to be decreasing
order of importance from the perspective of using the WKAPP products. Many of
these issues will be solved in future releases through improved data analysis
and a custom kinematic pipeline that is optimitized for WALLABY detections.

Figure 7: Comparison of the deprojected surface densities (dotted and dashed
coloured lines) recovered from WKAPP of a mock input surface density (solid
black line) in PDR1-like sources that are resolved by $D_{HI}=4$ beams and
$D_{HI}=8$ beams across their major axes at disk inclinations of $20^{\circ}$,
$60^{\circ}$, and $80^{\circ}$.

The most obvious issue in the kinematic modelling approach is the deprojection
of the surface density measured from the Moment 0 maps (see Sec. 4.4). The
standard $\cos(i)$ correction adopted here is known to fail at high
inclinations due to the line-of-sight passing through multiple radii (e.g.
Sancisi & Allen 1979). The failure of the $\cos(i)$ correction in the PDR1
regime is illustrated very clearly in Fig. 7, which shows an input deprojected
surface density distribution and the distributions that are recovered for
PDR1-like sources generated using MCGSuite at different resolutions and
inclinations using the WKAPP method.

Fig. 7 illustrates the well-known result that as the galaxy becomes more
highly inclined and more poorly resolved, the deprojected surface density from
ellipse fits to Moment 0 maps becomes much less reliable in both total flux
and profile shape. In particular, the inner profile peak is strongly
underestimated as the inclination increases and the resolution decreases. We
caution the user against using the inner deprojected surface density profile
unless the impact of these biases is quantified for the particular application
at hand.

Conversely, the outer profile profile shape is similar to the input one except
at the poorest resolution and highest inclination shown in Fig. 7. However,
the profile extent is biased to larger $R$ due to beam smearing. We judge the
outer deprojected surface densities and Hi sizes to be robust enough for use
in many cases, although Hi sizes should be first corrected for beam smearing
(e.g. Wang et al., 2021; Reynolds et al., 2021). In the future, the custom
WALLABY pipeline will fit the surface density using the full 3D cube and so,
should be more accurate than the ellipse fitting adopted here.

A second limitation is the restriction of the kinematic analysis to flat-disk
models. As described in Sec. 3.2, the homogeneous application of models to all
suitable PDR1 detections is a key principle of our modelling approach, which
drives the flat-disk modelling choice: in the marginally-resolved regime, it
is often not possible to reliably explore warps, non-radial flows, and other
complicated features due to a lack of statistically independent constraints on
the underlying structure. Certainly, some of the more well-resolved galaxies
in the sample show evidence for these complicated structures; for example, the
slight offset in the minor axis PV diagram for J100426-282638 in Fig. 6
indicates some level of non-circular motions. More sophisticated modelling of
these objects is likely warranted, and well-suited to 2D analyses where non-
axisymmetric structures can also be explored (e.g. Oh et al., 2018; Sellwood
et al., 2021). This work is underway for PDR1 sources.

A related issue is that complicated structural features may be present but not
spatially resolved, biasing the flat-disk models constructed here. The
importance of this bias is not known at present, but can be constrained by
convolving mock or real galaxies that exhibit such features down to the
marginally-resolved regime and exploring how well flat-disk models recover
their structure. While such tests have not yet been performed for PDR1
sources, they will be investigated in future data releases.

As a result of the second key principle that underpins our modelling approach
and in light of the lack of statistical uncertainties returned by available
tilted-ring algorithms (Sec. 3.2), the uncertainties on the optimized model
are derived from the differences between the 3DBarolo and FAT applications to
the pre-processed PDR1 cubelets (Sec. 4.2). While we judge these uncertainties
to be reasonable estimates of the reliability of the model parameters returned
that can be used for scientific applications, they have not been vetted as
statistical representations of the dispersion in model properties for the PDR1
galaxy population. As such, we recommend that they be considered as lower
limits of the absolute uncertainty on the properties of the underlying Hi
disks. The custom WALLABY pipeline that will be implemented in future data
releases will include robustly determined statistical uncertainties through
either monte-carlo or bootstrap-resampling approaches (e.g. Davis et al.,
2013; Sellwood et al., 2021).

Finally, we note that we have used the output PDR1 source cubelets from as
inputs to WKAPP. W22 discuss a number of data quality issues that may affect
the PDR1 release (see their section 4). The most significant from a kinematic
modelling perspective is likely the adoption of a $30^{\prime\prime}$ circular
Gaussian beam even for sources that may not have been deconvolved during
calibration and imaging. This issue is likely to affect the poorly-resolved,
lower-$S/N$ detections, which may be better characterized by the dirty beam,
which has beam dimensions of $\sim 30"\times 27"$. Moreover, signatures of the
dirty beam in the form of negative sidelobes may still remain in the cubelets.
Considering the small difference in dimensions between the dirty beam and the
restored beam, the different beam sizes are unlikely to affect our kinematic
models. It is more plausible that the presence of negative sidelobes from the
dirty beam have biased the modelled disk morphologies, but since the first
ASKAP sidelobe peaks at $\sim 5\%$ of the main lobe response and since the
integrated systematic effect on the measured fluxes is only of order $\sim
20\%$ (W22), we expect the effect on the disk morphologies to be mild. On
balance, we conclude that both of these effects are likely to be insignificant
relative to the other limitations in the kinematic models discussed above.

##  5 Kinematic Model Catalog

We have successfully generated WKAPP kinematic models for 124 PDR1 sources;
since 15/19 modelled Hydra TR1 sources also have models in Hydra TR2, WKAPP
has produced kinematic models for a total of 109 unique PDR1 objects. Table 4
lists the number of sources in each field and team release, the number of
sources for which modelling was attempted, and the number for which successful
models were obtained. Considering that we attempted to model 209/592 ($35\%$)
unique sources, our model success rate is $\sim 60\%$. The coloured points in
Fig. 1 summarize these results in the source size-$\log(S/N)$ plane. The mean
uncertainties on the geometric model parameters are listed in Table 5: we
typically constrain the kinematic centre to a few arcsec and $\textrm{km
s}^{-1}$, and the disk inclination and position angle to better than $\sim
5^{\circ}$ and $\sim 2^{\circ}$, respectively.

We note that the differences between the 15 unique objects in the Hydra field
for which for both the TR1 and TR2 detections have successful kinematic
models, the differences between them are generally small: the rotation curve
differences are typically smaller than the uncertainties due to inclination.
We recommend using the Hydra TR2 models over the Hydra TR1 models when both
are available.

| Hydra TR1 | Hydra TR2 | Hydra Unique | Norma TR1 | NGC 4636 TR1
---|---|---|---|---|---
$N_{\rm{sources}}$ | 148 | 272 | 301 | 144 | 147
$N_{\rm{attempted}}$ | 37 | 74 | 79 | 63 | 67
$N_{\rm{success}}$ | 19 | 31 | 35 | 31 | 43
Table 4: Number of PDR1 sources in each field (first row), the number for which WKAPP modelling was attempted (second row), and the number of successful WKAPP models (third row). Parameter | Mean Uncertainty
---|---
V_sys | $2.2~{}\textrm{km s}^{-1}$
Inc_model | $4.3^{\circ}$
PA_model_g | $1.5^{\circ}$
RA_model |  $2.4$′′
DEC_model |  $2.0$′′
Table 5: Mean uncertainties for the geometric parameters of optimized models.

To illustrate the conditions under which WKAPP succeeds and fails, Figure 8
shows moment maps of PDR1 sources in both of these categories. These include
both the highly resolved, high $S/N$ (ell_maj $\geq~{}7$ beams and
$\log(S/N_{\mathrm{obs}})\geq 1.6$) and the marginally resolved, low $S/N$
regimes (ell_maj $\leq~{}2.5$ beams and $\log(S/N_{\mathrm{obs}})\leq 1.4$).
The galaxies in the A-C panel pairs have morphologies consistent with rotating
disks and do not show strong signatures of non-circular motions, asymmetries,
disturbances, or other such features (although panel pair B does show spiral
arms and small non-circular motions). Given that our modelling method treats
galaxies as ‘flat’ disks, it is unsurprising that we successfully model this
type of high $S/N$, high resolution detection. Nonetheless, each galaxy in the
top row is a candidates for more detailed modelling in the future as they have
sufficient resolution elements to identify warps, non-circular flows, measure
the velocity dispersion, etc.

The high resolution, high $S/N$ failures in Fig. 8 are all interesting, as
each galaxy fails for a different reason. Panel pair D shows a galaxy with a
very complicated velocity profile that may be related to infalling Hi from a
recent interaction. The E panel pair shows a galaxy with an extended tidal
tail; as explained in W22, this cubelet may also contain deconvolution
artifacts. This galaxy may be modelled in the future with slightly more
careful masking/modelling. Panel pair F show a pair of interacting galaxies
that SoFiA detected as a single object. As the masking improves within WKAPP,
both those objects may be modelled in the future, but such modelling will be
challenging due to their interaction.

Figure 8: Moment 0 and Moment 1 maps for a sample set of PDR1 sources where
the WKAPP modelling is a success (first and third rows) or a failure (second
and fourth rows). The top two rows show well-resolved and high-$S/N$ sources,
while the bottom two rows have low resolution and $S/N$ values. These sources
are shown in Fig. 1 with the outlined symbols. The open ellipse in the Moment
0 maps shows the beam FWHM. Figure 1 indicates each of the galaxies as open
blue squares (top row), open red squares (second row), open blue diamonds
(third row), and open red diamonds (bottom row).

The low resolution, low $S/N$ rows are also quite interesting. Unlike their
higher resolution counterparts, it is more difficult to identify the reasons
for the specific modelling successes or failures. On the whole, however, we
find that the most common cause of a failure in this regime is that the
default source-finding mode in FAT or 3DBarolo is unable to find the source in
the cubelet. Additionally, both 3DBarolo and FAT have a number of default
quality control flags, which a low resolution, low $S/N$ source may not
satisfy. Another situation where the codes may fail is when the automated
initial FAT or 3DBarolo estimate of the model parameters is poor, which then
results in a poor fit. It is again important to note here that both FAT and
3DBarolo can be tuned individually to overcome these issues and produce
accurate models for some of these cases, but that we elect to run both codes
automatically and homogeneously on PDR1 sources (Sec. 3).

Perhaps the most surprising result is the number of marginally resolved
objects that have been modelled with only a few beams of resolution. This is a
testament to the power of 3D tilted ring modelling. Contrasting the low
resolution, low $S/N$ successes to the failures suggests that it is the
combination of low $S/N$ and low resolution that leads to the modelling
efforts failing for apparently similar objects (based on a visual comparison
of the moment maps). This suggests that the size and $S/N$ cuts applied to the
PDR1 detection sample will be sufficient for future modelling efforts.

Figure 9 illustrates the diversity of the modelled galaxy rotation curves and
surface density profiles across PDR1 detections. When source distances are
calculated from barycentric redshifts under the assumption of a flat
$\Lambda$CDM cosmology with a Hubble parameter of
$H_{0}=70~{}\mathrm{km\,s^{-1}\,Mpc^{-1}}$ and a local matter density of
$\Omega_{\rm m}=0.3$, the measured sizes range from a few to tens of
kiloparsecs. The rotation velocity amplitudes range from $30-250~{}\textrm{km
s}^{-1}$. Such a range in size, velocity, and from inference, mass, means that
this sample of kinematic models will be valuable for many different studies
and science questions. This sample is of comparable size and covers a similar
mass range as SPARC (Spitzer Photometry and Accurate Rotation Curves; Lelli et
al. 2016), which contains 175 rotation curves. SPARC and the PDR1 kinematic
release are highly complementary across a range of scientific applications:
while SPARC galaxies are generally better resolved than PDR1 sources, the PDR1
selection function is well-defined (W22).

Figure 9: The full sample of optimized model rotation curves (top panel) and
deprojected surface density profiles (bottom panel) for all pilot fields. The
blue, red, and orange lines show galaxies from the Hydra, Norma, and NGC 4636
fields respectively. The radial sizes are calculated using redshift derived
distances. The horizontal dashed line shows the $1$M⊙ pc-2 surface density
that defines the Hi radius of a galaxy. The right-hand panels are normalized
by $R_{\text{H\sc{i}},c}$ as determined from the surface density profiles.

The deprojected surface densities in the lower panels of Fig. 9 suggest that
there is an Hi surface density saturation at $\sim
10\textrm{M}_{\odot}~{}\textrm{pc}^{-2}$, consistent with the results of
Bigiel et al. (2008). However, we caution against using the inner surface
densities for scientific applications without further modelling, given the
breakdown in the standard $\cos(i)$ correction used to derive the deprojected
surface densities in the marginally-resolved regime (see Section 4.6 and Fig.
7).

By contrast, the outer deprojected surface density profiles are reliably
recovered by WKAPP, modulo being radially smeared by the beam (see Section 4.6
and Fig. 7). The outer profiles in Fig. 9 show the characteristic exponential
outer decline noted in previous work (e.g. Wang et al., 2016). We plot in
Figure 10 the diameter $D_{HI,c}$ at which the deprojected surface crosses
$1\,\mathrm{M_{\odot}\,pc^{-2}}$ as a function of ell_maj recovered by SoFiA
(see Fig. 1). We note that since both parameters are estimates of disk size
from the PDR1 Moment 0 maps, they should be similarly beam smeared.

The best fitting line for the data (performed in linear space) of $m=1.97$
shown in Fig. 10 illustrates that, in general, $D_{HI,c}$ exceeds ell_maj by a
factor of $\sim$two. ell_maj is computed from the second spatial moment of the
Moment 0 map along the major axis (Serra et al., 2015), which for a Gaussian
profile approaches twice its standard deviation. The factor of $\sim$2
difference between $D_{HI,c}$ and ell_maj then arises naturally from the outer
Gaussian profile shape provided it peaks in the range
$6\,\mathrm{M_{\odot}\,pc^{-2}}-10\,\mathrm{M_{\odot}\,pc^{-2}}$, which is
generally the case for the PDR1 sources plotted in Fig. 9. This difference
between ell_maj and $D_{\text{H\sc{i}}}$ justifies our PDR1 selection
criterion of ell_maj $\,>\,2$ beams (see Sec. 4.1). Figs. 10 and 1 imply that
the majority of successful kinematic models have $D_{\text{H\sc{i}}}\geq 3.5$
beams, consistent with the modelling tests of Lewis (2019).

Figure 10: A comparison of the $R_{\text{H\sc{i}}}$ radius to the SoFiA
ell_maj parameter for all successfully modelled PDR1 detections. The black
dashed line shows the one-to-one line, the red dashed line shows the best fit
straight line to the data, while the circle, star, and triangle symbols
indicate galaxies in the Hydra, Norma, and NGC 4636 fields respectively. The
values for $m$ are the slopes of the one-to-one and best fit lines in linear
space, respectively. The open symbols indicate the fitted galaxies (rows 1 and
3) shown in Fig. 8.

##  6 The Population of Kinematically Modelled Detections

A key question to ask when producing a survey is what are the biases in a
particular sample? In this case, are there any biases/selection effects that
apply to the kinematically modelled sample of galaxies relative to the larger
WALLABY sample?

To investigate the possibility of an environmental selection, we used the
Reynolds et al. (2022) dataset of Hydra TR2 WALLABY detections with velocities
$cz<7000$ km/s. In that work, Reynolds et al. (2022) classifed the galaxy
environment as field, infalling, or cluster. The top panel of Figure 11 shows
these galaxies along with their measured stellar and Hi masses, while the
bottom panel shows successful and failed models in the star formation rate –
stellar mass plane. We find no qualitative evidence that galaxies in different
environments are more or less likely to be modellable, though the sample is
relatively small to subdivide by environment. Morevoer, such environmental
effects are difficult to discern using only detections in HI-blind, shallow
surveys due to selection effects (Cortese et al., 2021).

Figure 11: The population of kinematically modelled PDR1 detections in the
Hydra field in the context of the Hi mass - stellar mass relation (top panel)
and star-formation - stellar mass relation (bottom panel) from (Reynolds et
al., 2022). In both panels, the symbol shape denotes the environmental
designation from Reynolds et al. (2022). Galaxies that have been successfully
kinematically modelled are plotted in blue, while those for which modelling
was attempted but failed are plotted in red. In the top panel, the blue line
shows the predicted locus of points for galaxies that lie on the Hi-mass - Hi-
diameter relation (Wang et al., 2016) with Hi diameters between 3 and 4 beams
across at the at the Hydra cluster distance ($D=60\,\mathrm{Mpc}$, Jørgensen
et al. 1996). The galaxies that were successfully modelled tend to lie above
that region, indicating that one of the main drivers of modellability is
angular size. There is no qualitative correlation between environment, SFR and
galaxy modellability in the sample examined.

However, Fig. 11 suggests that kinematic models of sources with $M_{\star}\leq
10^{8.5}\,M_{\odot}$ tend to fail, models of sources with $M_{\star}\geq
10^{9.5}\,M_{\odot}$ tend to succeed, and successfully modelled sources with
$M_{\star}\sim 10^{9}\,M_{\odot}$ tend to be more gas-rich than sources where
the models failed. These trends are a consequence of the Hi mass - Hi size
relation (Wang et al., 2016), the locus of which for Hydra cluster galaxies
spatially resolved by 3–4 beams is shown by the shaded region in the top
panel. Model successes tend to lie above the shaded region, while failures
tend to lie below it. This threshold is broadly consistent with the benchmark
results of Lewis (2019) for FAT, and demonstrates that galaxy angular sizes
(and $S/N$ by virtue of the correlation in Fig. 1) are the strongest
predictors of modellability with WKAPP among PDR1 sources.

##  7 Conclusions

WALLABY will use the ASKAP telescope to detect the Hi content of $\sim
210,000$ galaxies out to redshift $z\sim 0.1$. The PDR1 observations target
three fields observed at the full resolution and sensitivity of the planned
survey. The source-finding analysis applied to these fields detected 592 Hi
sources (W22). Of those, we have kinematically modelled 109 galaxies.

In our modelling approach (WKAPP), we attempt to fit all detections with
$\texttt{ell\\_maj}\geq 2$ beams or $\log(S/N_{obs}\geq 1.25$ using both
3DBarolo and FAT. There are 209 unique galaxies that meet this criteria. Both
the 3DBarolo and FAT analyses are constrained to consider only pure flat disk
models in order to obtain a uniform and robust population of modelled
galaxies. The results of the individual applications are examined visually to
determine their plausibility.

The optimized PDR1 models are generated by first averaging the geometric
parameters derived from the FAT and 3DBarolo fits. Then the inclination
corrected, interpolated rotation curves are averaged together to generate the
optimized rotation curve. Finally, the surface density is extracted from the
SoFiA masked Moment 0 map via ellipse fitting using the optimized model
geometry.

The full set of kinematic models are publicly available at the WALLABY pilot
phase data portal. The modelled population tends to be gas-rich and tends to
have larger stellar masses than the non-modelled population. This is largely
expected from the Hi mass-size relation.

The WKAPP modelling success rate is roughly $20\%$ (109/592). This success
bodes well for the full WALLABY dataset. The $\sim 20\%$ modelling suggests
that we will generate kinematic models for $\sim 40,000$ galaxies over the
full survey. However, the three PDR1 fields were chosen for testing purposes
and contain galaxies that may be less distant than the average WALLABY galaxy.
Given the importance of size and $S/N$ the full WALLABY modelling success rate
may be somewhat lower than $20\%$, but it is still likely higher than $10\%$
(Koribalski et al., 2020).

The use of WKAPP on the PDR1 sources has been quite successful. While the
modelling success rate across the full sample is only $\sim$20%, for sources
with ell_maj$\geq 2$ beams it is $\sim 50\%$. Additionally, galaxies with
ell_maj$\geq 2$ beams are less resolved than prior estimates of the FAT and
3DBarolo resolution limits, allowing us to attempt kinematic models on a
greater number of sources than initially expected. Beyond the successes of
WKAPP for PDR1 sample, it is a critical step in developing the full, automatic
pipeline that will be deployed for the full WALLABY survey. In the meantime,
these kinematic models are useful for a large variety of science
investigations. Moreover, examining the population of galaxies where the
models failed is also informative, and has revealed many intriguing and
complicated objects. While the kinematic models presented here are, by
necessity, relatively simple, there are a number of candidates for more
detailed 2D and 3D modelling. Comparing the WKAPP models to existing models
for some of these candidates as well as exploring more detailed 2D and 3D
modelling efforts will help to understand the strengths and weaknesses of this
approach. Another important exercise will be testing WKAPP, the future full
pipeline, and other kinematic modelling software packages, using mock WALLABY
observations from large cosmological simulations.

The PDR1 kinematic models presented here are the first step towards full set
of WALLABY kinematic models. We plan to publically release the rotation
curves, surface density profiles, and other properties for the $10000-40000$
galaxies that we expect to model. This will form a large, homogeneous legacy
data set that will allow explorations of the velocity function, the TF
relation, investigations of galaxy mass distributions, and much more.

We would like to thank the referee for their useful comments and suggestions
for the paper.

The Australian SKA Pathfinder is part of the Australia Telescope National
Facility (https://ror.org/05qajvd42) which is managed by CSIRO. Operation of
ASKAP is funded by the Australian Government with support from the National
Collaborative Research Infrastructure Strategy. ASKAP uses the resources of
the Pawsey Supercomputing Centre. Establishment of ASKAP, the Murchison Radio-
astronomy Observatory and the Pawsey Supercomputing Centre are initiatives of
the Australian Government, with support from the Government of Western
Australia and the Science and Industry Endowment Fund. We acknowledge the
Wajarri Yamatji as the traditional owners of the Observatory site.

WALLABY acknowledges technical support from the Australian SKA Regional Centre
(AusSRC) and Astronomy Data And Computing Services (ADACS).

This research used the facilities of the Canadian Astronomy Data Centre
operated by the National Research Council of Canada with the support of the
Canadian Space Agency.

This paper includes archived data obtained through the CSIRO ASKAP Science
Data Archive, CASDA (http://data. csiro.au).

This paper uses resources from the Canadian Initiative for Radio Astronomy
Data Analysis (CIRADA), which is funded by a grant from the Canada Foundation
for Innovation 2017 Innovation Fund (Project 35999) and by the Provinces of
Ontario, British Columbia, Alberta, Manitoba and Quebec, in collaboration with
the National Research Council of Canada, the US National Radio Astronomy
Observatory and Australia’s Commonwealth Scientific and Industrial Research
Organisation.

Part of this research was conducted by the Australian Research Council Centre
of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through
project number CE170100013.

AB acknowledges support from the Centre National d’Etudes Spatiales (CNES),
France. EDT was supported by the US National Science Foundation under grant
1616177. JMvdH acknowledges support from the European Research Council under
the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant
Agreement nr. 291531 (HIStoryNU). KS acknowledges support from the Natural
Sciences and Engineering Research Council of Canada (NSERC). LVM acknowledges
financial support from the State Agency for Research of the Spanish Ministry
of Science, Innovation and Universities through the ”Center of Excellence
Severo Ochoa” awarded to the Instituto de Astrofísica de Andalucía
(SEV-2017-0709), from grant RTI2018-096228-B-C31 (MCIU/AEI/FEDER,UE), from the
grant IAA4SKA (Ref. R18-RT-3082) from the Economic Transformation, Industry,
Knowledge and Universities Council of the Regional Government of Andalusia and
the European Regional Development Fund from the European Union. PK is
partially supported by the BMBF project 05A17PC2 for D-MeerKAT. SHOH
acknowledges a support from the National Research Foundation of Korea (NRF)
grant funded by the Korea government (Ministry of Science and ICT: MSIT) (No.
NRF-2020R1A2C1008706). TS acknowledges support by Fundação para a Ciência e a
Tecnologia (FCT) through national funds (UID/FIS/04434/2013), FCT/MCTES
through national funds (PIDDAC) by this grant UID/FIS/04434/2019 and by FEDER
through COMPETE2020 (POCI-01-0145-FEDER-007672). TS also acknowledges the
support by the fellowship SFRH/BPD/103385/2014 funded by the FCT (Portugal)
and POPH/FSE (EC). TS additionally acknowledges support from DL
57/2016/CP1364/CT0009 from The Centro de Astrofísica da Universidade do Porto.

This research uses Astropy,999http://www.astropy.org a community-developed
core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018).
It also uses the Numpy (Harris et al., 2020), SciPy (Virtanen et al., 2020),
PANDAS (Reback et al., 2020), and MatPlotLib (Hunter, 2007) libraries.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Battaner et al. (1990) Battaner, E., Florido, E., & Sanchez-Saavedra, M. L. 1990, A&A, 236, 1
  * Begeman (1987) Begeman, K. G. 1987, PhD thesis, University of Groningen, Netherlands
  * Bekiaris et al. (2016) Bekiaris, G., Glazebrook, K., Fluke, C. J., & Abraham, R. 2016, MNRAS, 455, 754
  * Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846
  * Bosma (1978) Bosma, A. 1978, PhD thesis, University of Groningen, Netherlands
  * Cortese et al. (2021) Cortese, L., Catinella, B., & Smith, R. 2021, PASA, 38, e035
  * Davis et al. (2013) Davis, T. A., Alatalo, K., Bureau, M., et al. 2013, MNRAS, 429, 534
  * de Blok (2010) de Blok, W. J. G. 2010, Advances in Astronomy, 2010, 789293
  * Di Teodoro & Fraternali (2015) Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021
  * Di Teodoro & Peek (2021) Di Teodoro, E. M., & Peek, J. E. G. 2021, ApJ, 923, 220
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hotan et al. (2021) Hotan, A. W., Bunton, J. D., Chippendale, A. P., et al. 2021, PASA, 38, e009
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Jones et al. (2021) Jones, G. C., Vergani, D., Romano, M., et al. 2021, MNRAS, 507, 3540
  * Jørgensen et al. (1996) Jørgensen, I., Franx, M., & Kjærgaard, P. 1996, MNRAS, 280, 167
  * Józsa et al. (2007) Józsa, G. I. G., Kenn, F., Klein, U., & Oosterloo, T. A. 2007, A&A, 468, 731
  * Józsa et al. (2009) Józsa, G. I. G., Oosterloo, T. A., Morganti, R., Klein, U., & Erben, T. 2009, A&A, 494, 489
  * Józsa et al. (2021) Józsa, G. I. G., Thorat, K., Kamphuis, P., et al. 2021, MNRAS, 501, 2704
  * Kamphuis et al. (2015) Kamphuis, P., Józsa, G. I. G., Oh, S. . H., et al. 2015, MNRAS, 452, 3139
  * Khoperskov et al. (2014) Khoperskov, S. A., Moiseev, A. V., Khoperskov, A. V., & Saburova, A. S. 2014, MNRAS, 441, 2650
  * Koribalski et al. (2020) Koribalski, B. S., Staveley-Smith, L., Westmeier, T., et al. 2020, Ap&SS, 365, 118
  * Krajnović et al. (2006) Krajnović, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, MNRAS, 366, 787
  * Lagos et al. (2018) Lagos, C. d. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573
  * Lelli et al. (2016) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157
  * Lelli et al. (2019) Lelli, F., McGaugh, S. S., Schombert, J. M., Desmond, H., & Katz, H. 2019, MNRAS, 484, 3267
  * Lewis (2019) Lewis, C. 2019, PhD thesis, Queen’s University at Kingston, Canada
  * McGaugh et al. (2000) McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99
  * Mutabazi (2021) Mutabazi, T. 2021, ApJ, 911, 16
  * Nelder & Mead (1965) Nelder, J. A., & Mead, R. 1965, Computer Journal, 7, 308
  * Oh et al. (2011) Oh, S.-H., de Blok, W. J. G., Brinks, E., Walter, F., & Kennicutt, Robert C., J. 2011, AJ, 141, 193
  * Oh et al. (2018) Oh, S.-H., Staveley-Smith, L., Spekkens, K., Kamphuis, P., & Koribalski, B. S. 2018, MNRAS, 473, 3256
  * Papastergis & Shankar (2016) Papastergis, E., & Shankar, F. 2016, A&A, 591, A58
  * Ponomareva et al. (2016) Ponomareva, A. A., Verheijen, M. A. W., & Bosma, A. 2016, MNRAS, 463, 4052
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C, 2nd edn. (Cambridge, USA: Cambridge University Press)
  * Reback et al. (2020) Reback, J., McKinney, W., jbrockmendel, et al. 2020, pandas-dev/pandas: Pandas 1.0.3, doi:10.5281/zenodo.3715232
  * Reynolds et al. (2021) Reynolds, T. N., Westmeier, T., Elagali, A., et al. 2021, MNRAS, 505, 1891
  * Reynolds et al. (2022) Reynolds, T. N., Catinella, B., Cortese, L., et al. 2022, MNRAS, 510, 1716
  * Rogstad et al. (1974) Rogstad, D. H., Lockhart, I. A., & Wright, M. C. H. 1974, ApJ, 193, 309
  * Sancisi & Allen (1979) Sancisi, R., & Allen, R. J. 1979, A&A, 74, 73
  * Schoenmakers et al. (1997) Schoenmakers, R. H. M., Franx, M., & de Zeeuw, P. T. 1997, MNRAS, 292, 349
  * Sellwood et al. (2021) Sellwood, J. A., Spekkens, K., & Eckel, C. S. 2021, MNRAS, 502, 3843
  * Serra et al. (2015) Serra, P., Westmeier, T., Giese, N., et al. 2015, MNRAS, 448, 1922
  * Spekkens & Sellwood (2007) Spekkens, K., & Sellwood, J. A. 2007, ApJ, 664, 204
  * Starkman et al. (2018) Starkman, N., Lelli, F., McGaugh, S., & Schombert, J. 2018, MNRAS, 480, 2292
  * Stevens et al. (2016) Stevens, A. R. H., Croton, D. J., & Mutch, S. J. 2016, MNRAS, 461, 859
  * Swaters et al. (2002) Swaters, R. A., van Albada, T. S., van der Hulst, J. M., & Sancisi, R. 2002, A&A, 390, 829
  * Tamburro et al. (2009) Tamburro, D., Rix, H. W., Leroy, A. K., et al. 2009, AJ, 137, 4424
  * Tully & Fisher (1977) Tully, R. B., & Fisher, J. R. 1977, A&A, 54, 661
  * Tully et al. (2013) Tully, R. B., Courtois, H. M., Dolphin, A. E., et al. 2013, AJ, 146, 86
  * van Albada et al. (1985) van Albada, T. S., Bahcall, J. N., Begeman, K., & Sancisi, R. 1985, ApJ, 295, 305
  * van Albada & Sancisi (1986) van Albada, T. S., & Sancisi, R. 1986, Philosophical Transactions of the Royal Society of London Series A, 320, 447
  * van der Hulst et al. (1992) van der Hulst, J. M., Terlouw, J. P., Begeman, K. G., Zwitser, W., & Roelfsema, P. R. 1992, in Astronomical Society of the Pacific Conference Series, Vol. 25, Astronomical Data Analysis Software and Systems I, ed. D. M. Worrall, C. Biemesderfer, & J. Barnes, 131
  * Verheijen (2001) Verheijen, M. A. W. 2001, ApJ, 563, 694
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Walter et al. (2008) Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563
  * Wang et al. (2016) Wang, J., Koribalski, B. S., Serra, P., et al. 2016, MNRAS, 460, 2143
  * Wang et al. (2021) Wang, J., Staveley-Smith, L., Westmeier, T., et al. 2021, ApJ, 915, 70
  * Westmeier et al. (2021) Westmeier, T., Kitaeff, S., Pallot, D., et al. 2021, MNRAS, 506, 3962
  * Westmeier et al. (2022) Westmeier, T., Deg, N., Spekkens, K., et al. 2022, PASA, doi:doi.org/10.1017/pasa.2022.50
  * Whiting (2012) Whiting, M. T. 2012, MNRAS, 421, 3242