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

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 optimised model parameters and uncertainties. The 109 models presented here tend to be gas rich detections resolved by at least 3–4 synthesised 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 analysed dataset will be a powerful legacy product that that will enable a wide range of scientific studies.

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} 210000$ 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, McGaugh, & Schombert 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, Verheijen, & Bosma 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, Florido, & Sanchez-Saavedra 1990; Stevens, Croton, & Mutch 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, Lockhart, & Wright 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 Section 2 we briefly describe the WALLABY detections, but a full description is provided in the data release paper W22. Section 3 describes tilted-ring modelling in general, while Section 4 describes the specific approach taken for the PDR1 observations. Section 5 provides the overall results of the PDR1 kinematic models. Section 6 describes the population of kinematically modelled galaxies, and Section 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, Franx & Kjærgaard 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 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^′′ (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 Figure 1 as the SoFiA-returned major axis diameter ell_maj of an ellipse fitted to the source Moment 0 map (Westmeier 2021). As in W22, we compute the integrated $S/N$ via:

(1) \begin{equation} S/N_{obs}=\frac{S_{mask}}{\sigma_{rms} \sqrt{N_{mask} \Omega}},\end{equation}

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 Figure 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 Equation (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 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 Figure 8.

Figure 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 Section 3, and describe the resulting procedure in Section 4.

3. Kinematic modelling considerations

Given that many science goals for WALLABY are enabled by statistical samples of resolved source properties (see Section 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 Section 2, we discuss here the considerations that drive our kinematic modelling procedure. Section 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. Section 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 Section 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 optimise 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 characterise non-circular flows or complex disk geometries have since been developed and publicly released (e.g. reswri, Schoenmakers, Franx, & de Zeeuw 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 generalised 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 Figure 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

FAT a (Kamphuis et al. 2015) automates the application of TiRiFiC b (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 optimised 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 initialise 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$ to $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 Figure 1 (see Section 5). We work with FAT version 2.01.

3.1.2. 3DBarolo

3DBarolo c (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 initialisation, model optimisation, and flux normalisation. 3DBarolo can use a built-in source finder based on DUCHAMP (Whiting 2012) to initialise the models, or the user can specify initial parameter estimates directly. Once the source(s) are found, the model is optimised 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 normalised 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 optimisation 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 optimisations. 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 Figure 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.

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 initialisation can both strongly influence the optimisations 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 initialises parameters, hence the importance of that modelling choice. The parameter initialisation, in turn, can impact the model outputs because optimisation 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 Figure 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 initialisation (either initialised automatically by the code or initialised by the user to the PDR1 source values from W22):

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

• Barolo_smooth_auto: 3DBarolo applied to the spectrally smoothed cubelet, with automated parameter initialisation;

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

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

• FAT_smooth_auto: FAT applied to the spectrally smoothed cubelet, with automated parameter initialisation;

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

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.

We note that since FAT does not allow the user to initialise parameters, there no such model with this option in Figure 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 initialise geometric parameters with photometric values. We note that the Barolo_auto_smooth_vdisp shown in Figure 2 involves a different fitting mode than the other models, and is discussed further below.

Comparing the optimised rotation curves (Figure 2A), surface density profiles (Figure 2B) and disk geometries (Figure 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 initialisation 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 centre (Figure 2D–F) are recovered well. Given the pixel size is ${\sim} 6 ''$ , the kinematic centre 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 centre and position angle for various fits, it indicates a failure in one or more of those fits. The systemic velocity itself (Figure 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 centre.

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 optimised values when both models succeed (e.g. FAT_full_auto and FAT_smooth_auto in Figure 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 Figure 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 (Figure 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 optimisation methodology (see Section 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.

The dashed horizontal lines in Figure 2 plot the PDR1 source parameters that best approximate the geometric parameters returned by the kinematic modelling: $\cos^{-1}$ (ell_min/ell_maj), kin_pa, ra, dec, and freq (see Table 3 of W22) converted to the appropriate units are shown in Figure 2C-G respectively. Model Barolo_smooth_source initialises 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 initialised 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 optimiser 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 initialisation in the kinematic analysis.

We now discuss model Barolo_auto_smooth_vdisp in Figure 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 Figure 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 Figure 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 \le$ ell_maj $\le 4 $ and $1.25\le \log(S/N_{obs})\le 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} \le \sigma \le 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.

Figure 3. Velocity dispersion profiles from models identical to Barolo_ auto_smooth_vdisp in Figure 2, applied to all 36 PDR1 sources with $2\le$ ell_maj $\le 4 $ and $1.25\le \log(S/N_{obs})\le 1.5$ . The profiles are coloured according to the model disk inclination.

4. The WALLABY kinematic analysis proto-pipeline

Having explored some key considerations for kinematic modelling of PDR1 detections in Section 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 page.d Figure 4 summarises 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 optimised models and uncertainties by averaging successful model fits;

4. 4. Compute surface density profiles from PDR1 source Moment 0 maps using optimised 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 Sections 4.1–4.4, the model outputs in Section 4.5, and some limitations of the current approach in Section 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{Hi}}$ in our successful models (see Section 6 for a comparison of $D_{\text{Hi}}$ to ell_maj). Therefore we attempt to model all detections with ell_maj $\ge 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})\ge 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 Figure 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 Section 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 Section 3.2, we allow each code to automatically initialise 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 $\mathrm{beam}^{-1}$ and we use the SEARCH source-finding method and the azimuthal normalisation 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 $\mathrm{beam}^{-1}$ 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 Section 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 Section 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 Figure 2; see Section 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 Figure 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 (centre, $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 Figure 1).

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.

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 optimise 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.

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.

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 Figure 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 optimised rotation curves will have similar levels of self-consistency as they are generated by averaging the interpolated, inclination corrected FAT and 3DBarolo outputs. However, the visual inspection of the different fits as well as a final examination of the optimised model help to avoid any inconsistencies, and in practice the best fitting disk centres and position angles are typically very similar (see Section 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 optimised model. This approach is similar to the 3DBarolo procedure for calculating surface densities but differs strongly from the FAT approach (see Section 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 Section 4.6).

The optimised 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 Figure 1). In addition, we do not attempt to correct the outer surface density profiles for beam smearing effects. We discuss both effects in Section 4.6, and recommend that their impact be considered when using the corrected surface density profiles.

4.5. Model outputs

Every PDR1 source that is successfully modelled by WKAPP is characterised 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.

Table 2. WKAPP model parameters.

${}^{\rm a}\,$ In pixel coordinates relative to the pre-processed cubelet, which starts from the point (1,1).

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 Section 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). Figure 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 Section 4.6 below.

We also provide a quality flag for each model:

• QFlag_model = 0: No obvious issues.

• QFlag_model = 1: Inc_model $\leq 20^{\circ}$ or Inc_model $\geq 80^{\circ}$ .

• QFlag_model = 2: e_Vsys_model $\geq 15\,\textrm{km s}^{-1}$ .

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

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 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 Section 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.

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 Section 4.2. Several data and model cubelets are also provided as data products. The model cubelets are realisations of the optimised models using the stand-alone tilted-ring model generator MCGSuite codee (Lewis 2019, Spekkens et al. in prep). The pre-processed PDR1 cubelets to which FAT and 3DBarolo are applied (see Section 4.2) are in files with suffix _ProcData.fits. Realisations of the optimised 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 realisations with the properties of those cubelets (suffix _FullResModelCube.fits).

Table 3. WKAPP data products available for each successfully modelled PDR1 source.

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 Figure 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 optimised model. The middle left panel shows the PDR1 Moment 0 map and the location of the model centre 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 realisation), 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)f and the Canadian Astronomy Data Centre (CADC).g A full description of the data access options can be found in both W22 and through the data access portal.h

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 optimised for WALLABY detections.

The most obvious issue in the kinematic modelling approach is the deprojection of the surface density measured from the Moment 0 maps (see Section 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 Figure 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.

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$ .

Figure 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 Figure 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 Section 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 Figure 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 (Section 3.2), the uncertainties on the optimised model are derived from the differences between the 3DBarolo and FAT applications to the pre-processed PDR1 cubelets (Section 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^′′ 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 characterised 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 Figure 1 summarise 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.

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).

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.

Table 5. Mean uncertainties for the geometric parameters of optimised 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 ${\ge}7$ beams and $\log(S/N_{\mathrm{obs}})\ge 1.6$ ) and the marginally resolved, low $S/N$ regimes (ell_maj ${\le}2.5$ beams and $\log(S/N_{\mathrm{obs}})\le 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.

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 Figure 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 high resolution, high $S/N$ failures in Figure 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 artefacts. 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.

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 (Section 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 optimised 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 \textit{M}_{\odot}\ \mathrm{pc}^{-2}$ surface density that defines the Hi radius of a galaxy. The right-hand panels are normalised by $R_{\text{Hi},c}$ as determined from the surface density profiles.

The deprojected surface densities in the lower panels of Figure 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 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 Figure 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 Figure 7). The outer profiles in Figure 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 Figure 1). We note that since both parameters are estimates of disk size from the PDR1 Moment 0 maps, they should be similarly beam smeared.

Figure 10. A comparison of the $R_{\text{Hi}}$ 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 Figure 8.

The best fitting line for the data (performed in linear space) of $m=1.97$ shown in Figure 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– $10\,\mathrm{M_\odot\,pc^{-2}}$ , which is generally the case for the PDR1 sources plotted in Figure 9. This difference between ell_maj and $D_{\text{Hi}}$ justifies our PDR1 selection criterion of ell_maj $\,>\,2$ beams (see Section 4.1). Figures 10 and 1 imply that the majority of successful kinematic models have $D_{\text{Hi}} \ge 3.5$ beams, consistent with the modelling tests of Lewis (2019).

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\ \mathrm{km s}^{-1}$ . In that work, Reynolds et al. (2022) classified 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. Moreover, such environmental effects are difficult to discern using only detections in HI-blind, shallow surveys due to selection effects (Cortese, Catinella, & Smith 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, Figure 11 suggests that kinematic models of sources with $M_\star \le 10^{8.5}\,{\rm M}_\odot$ tend to fail, models of sources with $M_\star \ge 10^{9.5}\,{\rm M}_\odot$ tend to succeed, and successfully modelled sources with $M_{\star} \sim 10^9\,{\rm 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 Figure 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} 210000$ 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 ell_maj $\ge 2$ beams or $\log(S/N_{obs}\ge 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 optimised 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 optimised rotation curve. Finally, the surface density is extracted from the SoFiA masked Moment 0 map via ellipse fitting using the optimised 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} 40000$ 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 ${\ge} 2$ beams it is ${\sim} 50\%$ . Additionally, galaxies with ell_maj ${\ge} 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.