# Neutral versus ionized gas kinematics at z 2.6: The AGN-host starburst galaxy PKS 0529-549

Federico Lelli, Carlos de Breuck, Theresa Falkendal, Filippo Fraternali,

Allison W. S. Man, Nicole P. H. Nesvadba, Matthew D. Lehnert

► To cite this version:

Federico Lelli, Carlos de Breuck, Theresa Falkendal, Filippo Fraternali, Allison W. S. Man, et al. Neutral versus ionized gas kinematics at z 2.6: The AGN-host starburst galaxy PKS 0529-549. Monthly Notices of the Royal Astronomical Society, 2018, 479 (4), pp.5440-5447. 10.1093/mnras/sty1795. hal-01851234

HAL Id: hal-01851234

https://hal.science/hal-01851234v1

Submitted on 10 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Neutral versus ionized gas kinematics at z ≃ 2.6: the AGN-host starburst galaxy PKS 0529-549

Federico Lelli,<fnref n="1" />†* Carlos De Breuck,<fnref n="1" /> Theresa Falkendal,<fnref n="1" />,<fnref n="2" /> Filippo Fraternali,<fnref n="3" />

Allison W. S. Man,<fnref n="1" /> Nicole P. H. Nesvadba<fnref n="4" /> and Matthew D. Lehnert<fnref n="2" />

Accepted 2018 June 28. in original form 2018 June 4

## ABSTRACT

We present a multiwavelength study of the AGN-host starburst galaxy PKS 0529-549 at $z \simeq 2.6$. We use (1) new ALMA observations of the dust continuum and of the [CⅠ] 370 µm line, tracing molecular gas, (2) SINFONI spectroscopy of the [O III] 5007 Å line, tracing ionized gas, and (3) ATCA radio continuum images, tracing synchrotron emission. Both [CⅠ] and [O III] show regular velocity gradients, but their systemic velocities and position angles differ by ~300 km s⁻<fnref n="1" /> and ~30°, respectively. The [CⅠ] is consistent with a rotating disc, aligned with the dust and stellar continuum, while the [O III] likely traces an outflow, aligned with two AGN-driven radio lobes. We model the [CⅠ] cube using 3D disc models, which give best-fitting rotation velocities $V_{rot} \simeq 310$ km s⁻<fnref n="1" /> and velocity dispersions $\sigma_v \lesssim 30$ km s⁻<fnref n="1" />. Hence, the [CⅠ] disc has $V_{rot}/\sigma_v \gtrsim 10$ and is not particularly turbulent, similar to local galaxy discs. The dynamical mass (~10¹<fnref n="1" /> M☉) is comparable to the baryonic mass within the errors. This suggests that baryons dominate the inner galaxy dynamics, similar to massive galaxies at $z \simeq 0$. Remarkably, PKS 0529-549 lies on the local baryonic Tully-Fisher relation, indicating that at least some massive galaxies are already in place and kinematically relaxed at $z \simeq 2.6$. This work highlights the potential of the [CⅠ] line to trace galaxy dynamics at high $z$, as well as the importance of multiwavelength data to interpret gas kinematics.

**Key words:** galaxies: active–galaxies: evolution–galaxies: individual: PKS 0529-549–galaxies: kinematics and dynamics–galaxies: starburst–(*cosmology*:) dark matter.

## 1 INTRODUCTION

Over the past decade much progress has been made in the study of gas dynamics at high redshifts. Integral field spectroscopy (IFS) in the near-infrared (NIR) has been widely used to trace the ionized gas kinematics in star-forming galaxies at $z \simeq 1–3$ (e.g. Förster Schreiber et al. 2009; Gnerucci et al. 2011; Wisnioski et al. 2015; Stott et al. 2016; Nesvadba et al. 2017). Radio and submillimetre observations with the Jansky Very Large Array and the Plateau de Bure Interferometer have proved to be powerful in tracing CO kinematics up to $z \simeq 4$ (e.g. Hodge et al. 2012; Ivison et al. 2012; Tacconi et al. 2013). Moreover, the Atacama Large Millimeter Array (ALMA) has opened a new window to study galaxy dynamics at even higher redshifts ($z \simeq 4–7$) using the [CⅡ] emission line (e.g.

De Breuck et al. 2014; Jones et al. 2017; Shao et al. 2017; Smit et al. 2018).

These different studies probe different gas phases. NIR spectroscopy probes emission lines like $H\alpha$ and [O III], tracing warm ionized gas (~$10^4$ K). Radio and submillimetre observations, instead, probe CO, [CⅠ], and [CⅡ] transitions, tracing cold neutral gas (<100 K) like molecular hydrogen ($H_2$).

In nearby galaxies, $H_2$ and HⅠ generally display regular disc rotation, whereas $H\alpha$ and [O III] kinematics may be more complex. For example, in galaxies hosting starbursts and/or active galactic nuclei (AGNs), $H\alpha$ and [O III] are often associated with galactic winds (e.g. Harrison et al. 2014; Arribas et al. 2014). In high-$z$ galaxies, it is unclear how these different gas phases relate to each other and whether they share the same kinematics, since only a few galaxies have been studied using multiple gas tracers (Chen et al. 2017; Übler et al. 2018). In this context, the redshift window $z \simeq 2–3$ is very interesting because (1) both ionized and neutral gas can be observed and spatially resolved with existing facilities, and (2) it

<footnote n="1">¹European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="2">²Sorbonne Université, CNRS UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3">³Kapteyn Astronomical Institute, University of Groningen, Postbus 800, NL-9700 AV, Groningen, the Netherlands</footnote>

<footnote n="4">⁴Institut d'Astrophysique Spatiale, UMR 8617, Université Paris-Sud, Bât. 121, F-91405, Orsay, France</footnote>

<footnote>* ESO fellow.</footnote>

<footnote>† E-mail: flelli@eso.org</footnote>
corresponds to the peak of the cosmic star formation history (SFH; Madau & Dickinson 2014), a key epoch to understand the formation of massive galaxies.

Here, we compare neutral and ionized gas kinematics in the star-burst galaxy PKS 0529-549 at $z \simeq 2.6$, combining new [C I] (2–1) data from ALMA with IFS data from SINFONI (Nesvadba et al. 2017). PKS 0529-549 is a well-studied radio galaxy with an exquisite set of ancillary data: optical spectroscopy, ground-based NIR imaging, and radio polarimetry (Broderick et al. 2007), Spitzer/IRAC imaging (De Breuck et al. 2010), and 1.1 mm data from AzTEC (Humphrey et al. 2011). PKS 0529-549 hosts a Type-II AGN and shows two radio lobes. The eastern lobe holds the record for the highest Faraday rotation measure to date, implying a strong magnetic field and/or a dense circumgalactic medium (Broderick et al. 2007). With only 5 min of ALMA on-source integration time, we spatially resolve the [C I] emission of PKS 0529-549. We show that [C I] and [O III] display different kinematics and discuss possible interpretations with the aid of 3D kinematic models.

Throughout this paper, we assume a flat $\Lambda$ cold dark matter (ΛCDM) cosmology with $H_0 = 67.8 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.308$, and $\Omega_\Lambda = 0.692$ (Planck Collaboration XIII 2016). In this cosmology, $1 \text{ arcsec} = 8.2 \text{ kpc at } z = 2.6$.

# 2 DATA ANALYSIS

## 2.1 ALMA data

PKS 0529-549 was observed on 2014 September 2 during ALMA Cycle 2 with 34 working antennas, giving maximum and minimum baselines of 1090 and 33 m, respectively. The on-source integration time was 5 min. The flux and bandpass calibrator was J0519-4546, while the phase calibrator was J0550-5732. We used four 1.875 GHz-wide spectral windows centred at 224.99, 226.87, 239.99, and 241.87 GHz. The second spectral window contains the redshifted [C I] (2–1) line, which corresponds to the $^3P_2 \rightarrow ^3P_1$ transition and has a rest-frame frequency of 809.3435 GHz.

The data reduction was performed using the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007). The uv data were flagged and calibrated using the standard CASA pipeline. The imaging was performed using the task *clean* with *Robust* = 2, which corresponds to natural weighting of the visibilities. This provides the highest sensitivity but the lowest spatial resolution. The resulting synthesized beam is $0.43 \text{ arcsec} \times 0.28 \text{ arcsec}$ with a position angle (PA) of 64.5°. A continuum image was obtained combining all four spectral windows, reaching a sensitivity of $0.055 \text{ mJy beam}^{-1}$. The [C I] cube was imaged using a channel width of $50 \text{ km s}^{-1}$, giving a sensitivity of $0.55 \text{ mJy beam}^{-1}$ per channel. The continuum was subtracted from the [C I] cube in the image plane, after fitting a first-order polynomial to line-free channels. We also attempted continuum subtraction in the uv plane (before the Fourier transform), but sanity checks revealed that the continuum was oversubtracted in this way. In any case, the continuum subtraction does not strongly affect our kinematic analysis. The [C I] flux is $2.1 \pm 0.4 \text{ Jy km s}^{-1}$, as estimated from fitting the global [C I] profile (presented in Man et al. in preparation).

After imaging, the [C I] cube was analysed using the Groningen Imaging Processing System (GIPSY, Vogelaar & Terlouw 2001). To build the total [C I] map (moment zero), we sum the signal within a Boolean mask, which follows the kinematic structure of the [C I] emission. Given the low signal-to-noise ratio (S/N), the mask is constructed interactively on a channel-by-channel basis by defining regions in which contiguous pixels have $S/N \gtrsim 2$. When using a

mask, the noise in the total map is not uniform because a different number of channels is summed at each spatial pixel. It is however possible to construct S/N maps and define a pseudo-$3\sigma$ contour by considering spatial pixels with S/N between 2.75 and 3.25 (see the appendix in Lelli, Verheijen & Fraternali 2014). We measure pseudo-$3\sigma \simeq 150 \text{ mJy beam}^{-1} \text{ km s}^{-1}$.

To build the velocity map (moment one), we estimate an intensity-weighted mean velocity for the pixels inside the Boolean mask. We stress that the velocity map is very uncertain due to beam smearing effects and the low S/N. This map merely provides a rough overview of the gas kinematics. We do not consider the dispersion map (moment two) because it is even more strongly affected by beam smearing: moment-two values are not indicative of the intrinsic gas velocity dispersion because the [C I] line is largely broadened by the low spatial and spectral resolutions. Our kinematic analysis uses the full 3D information and takes resolution effects into account: we build model cubes and compare them with the observed one using position-velocity (PV) diagrams and channel maps (Section 3.2).

## 2.2 SINFONI and ATCA data

Nesvadba et al. (2017) present IFS data of PKS 0529-549 using SINFONI at the ESO Very Large Telescope (VLT). Two different datacubes were obtained: one at low spatial resolution ($1.2 \text{ arcsec} \times 1.2 \text{ arcsec}$), where the stellar continuum and several emission lines are detected, and one at high spatial resolution ($0.7 \text{ arcsec} \times 0.6 \text{ arcsec}$), where only the strongest emission lines ([O III] $\lambda4958$, [O III] $\lambda5007$, and $H\alpha$) are recovered. The $H\alpha$ line is blended with the [N II] doublet, creating a broad ($\gtrsim 100 \text{ Å}$) and asymmetric line. We refer to Nesvadba et al. (2017) for further details.

We use the GIPSY package to analyse the high-resolution cube. To increase the S/N of the emission lines, the cube was box-averaged over 3 × 3 pixels. We constructed emission-line maps by summing the signal within narrow wavelength ranges, which were determined by visual inspection excluding channels contaminated by sky lines. No continuum subtraction was attempted, but we used the low-resolution datacube to estimate that the net continuum contribution on the emission-line maps is negligible. We constructed velocity fields by fitting a Gaussian function at each spatial pixel. This is preferable to intensity-weighted mean velocities given the relatively high S/N. Each emission line was fitted independently considering a narrow wavelength range and rejecting fit results with unphysically small amplitudes ($\lesssim 5$ counts) and/or dispersions ($\lesssim 5 \text{ Å}$). Given the overall data quality, it would be untrustworthy to fit multiple Gaussians to try de-blending the $H\alpha$ and [N II] lines. Hereafter, we consider only the [O III] $\lambda5007$ velocity field because the $H\alpha$+[N II] and [O III] $\lambda4958$ ones have lower quality due to the lower S/N. These emission lines, however, display similar kinematics.

The absolute positional accuracy of the SINFONI data relative to the World Coordinate System is $\sim 1 \text{ arcsec}$ (see Nesvadba et al. 2017 for details). Since the ALMA positional accuracy is much higher than 1 arcsec, we improve the SINFONI astrometry using the following approach. It is sensible to assume that the SINFONI stellar continuum and the ALMA dust continuum are physically related: the two maps indeed display a similar elongated morphology with a $\text{PA} \simeq 75^\circ$ (Fig. 2.2). Thus, we overlay the emission peak of the SINFONI stellar continuum with that of the ALMA dust continuum. The resulting spatial shift is smaller than 1 arcsec.

We also consider radio continuum images from the Australia Telescope Compact Array (ATCA), which were provided by
Broderick et al. (2007). These authors published natural-weighted images at various frequencies and resolutions. Here, we show the uniform-weighted image at 18.5 GHz because it has the highest spatial resolution (0.55 arcsec × 0.37 arcsec). This image is published here for the first time.

# 3 RESULTS

## 3.1 Multiwavelength overview

Fig. 2.2 provides a multiwavelength view of PKS 0529-549. The top panels present the ALMA data. The ALMA continuum (left-hand panel) is studied in detail in Falkendal et al. (in preparation) by modelling its spectral energy distribution (SED). It consists of two components: a central one due to dust emission and a secondary one to the east due to synchrotron emission. The secondary component, indeed, nicely overlaps with the eastern radio lobe from ATCA observations. The [C I] emission (middle panel) coincides with the dust emission and displays a similar elongated morphology with PA ~ 75°. Both dust and [C I] emissions are spatially resolved with ~2 beams along the major axis, so their morphology is not strongly driven by the elongated beam shape, which has a similar PA of ~64.5°. No [C I] emission is associated with the synchrotron lobes. The [C I] velocity field (right-hand panel) shows a velocity gradient along the major axis: this is consistent with a rotating disc as we show in Section 3.2.

The bottom panels of Fig. 2.2 present the SINFONI data. The stellar continuum (left-hand panel) is barely resolved but it seems to display roughly the same PA as the dust and [C I] emissions. The ionized gas (middle panel) appears more extended than any other galaxy component. The SINFONI point spread function (PSF) is significantly larger than the ALMA beam, but if we smooth the ALMA maps to the same spatial resolution (0.7 arcsec × 0.6 arcsec), they remain less extended than the ionized gas. The [O III] λ5007 velocity field (right-hand panel) displays a velocity gradient with PA~105°, so it is tilted by ~30° with respect to the [C I] distribution and kinematics. The [O III] velocity gradient, instead, aligns well with the orientation of the two synchrotron lobes. Strikingly, the [O III] velocities are systematically higher than the [C I] velocities by ~300 km s⁻¹.

The Hα displays similar kinematics as the [O III], but its velocity field (not shown) is very uncertain due to the low S/N and contamination from the [N II] lines. In Man et al. (in preparation), we present an X-Shooter spectrum of PKS 0529-549 revealing that the [O III] and Hα velocities are consistent with other emission lines from ionized gas like [OII] and [CII], while the [C I] velocities are consistent with photospheric absorption lines from young stars.

Given the evidence above, we interpret the [C I] emission as a rotating gas disc, having a similar PA as both the dust component and the stellar continuum. The [O III] and the Hα emission, instead, could trace a gas outflow, being aligned with the synchrotron lobes. Nesvadba et al. (2017) provide additional arguments for the outflow interpretation. For example, the [O III] λ5007 profiles are very broad, ranging from ~300 to ~900 km s⁻¹ (see their fig. A.5): this is likely due to the contribution of different velocity components along the line of sight as expected in a gas outflow. The velocity difference between [C I] and ionized gas may indicate that the outflow does not propagate exactly from the centre of the [C I] disc, but from its northeastern, redshifted side. Alternatively, we may be seeing only the redshifted side of the flow (further away from the observer) because the ionized bi-cone may have asymmetric density distribution, as observed in some local starburst galaxies (e.g. NGC 253, Boomsma

et al. 2005) and high-z ones (e.g. TN J1338-1942, Zirm et al. 2005). Interestingly, if only SINFONI data were available, it would be much more difficult to distinguish the [O III] velocity field of an outflow from that of a rotating disc. This highlights the importance of having multiwavelength observations to interpret ionized gas data. This point is further discussed in Section 4.1.

## 3.2 Disc modelling

To test the hypothesis of a rotating [C I] disc, we use two different codes: GALMOD in GIPSY (Sicking 1997) and 3DBAROLO (Di Teodoro & Fraternali 2015). Both codes build 3D kinematic models of rotating discs, taking projection and resolution effects into account. The main difference between the two is that GALMOD needs to be used iteratively by trial and error, while 3DBAROLO explores the full parameter space to find a best-fitting solution by χ² minimization. In practice, we use GALMOD to set the initial estimates of 3DBAROLO and to fix several geometric parameters that are otherwise degenerated.

The disc is built using a set of circular rings. Each ring is characterized by four projection parameters (centre, PA, inclination *i*, and systemic velocity $V_{\text{sys}}$) and four physical ones (surface density $\Sigma$, vertical scale height *h*, rotation velocity $V_{\text{rot}}$, and velocity dispersion $\sigma_{\text{v}}$). The disc is then projected onto the sky to generate model cubes that are smoothed to the same spatial and spectral resolutions of the observations. Since the [C I] distribution is resolved with only ~2 beams, we assume that the disc parameters are constant with radius (*R*). This is equivalent to using a single, large ring. Effectively, the physical parameters are intensity-weighted average values over the semimajor axis of the [C I] disc.

We explored two basic models: a disc with uniform surface density (the ring extends from $R=0$ to $R \simeq 0.45$ arcsec) and a disc with a central hole (the ring extends from $R \simeq 0.20$ arcsec to $R \simeq 0.45$ arcsec). Both models give comparable results because the ALMA beam (0.43 arcsec × 0.28 arcsec) smooth out any intrinsic structure in the gas density. This degeneracy only affects the determination of $\Sigma$, which cannot be robustly estimated and is not used in our analysis. For similar reasons, the model is insensitive to the vertical density distribution: we assume an exponential vertical profile with fixed $h = 0.04$ arcsec $\simeq 300$ pc.

We use GALMOD following Lelli et al. (2015): we inspect channel maps and PV diagrams to obtain initial estimates of the disc parameters, which are then refined until we find a good match between model and observations. The centre, PA, and $V_{\text{sys}}$ are easy to estimate and kept fixed during subsequent iterations. The remaining parameters (*i*, $V_{\text{rot}}$, and $\sigma_{\text{v}}$) are more challenging to determine due to degeneracies. The inclination is estimated by comparing model [C I] maps (after beam convolution) to the observed one. Fig. 2 (top panels) shows that the observed [C I] map is well reproduced by our model map: the residuals are smaller than ~3σ and less extended than the beam. These residuals are probably due to deviations from a pure axisymmetric gas distribution in the central regions. Clearly, the [C I] disc is not face-on ($i > 30°$), so the net correction on $V_{\text{rot}}$ is relatively small going as $1/\sin(i)$.

Finally, $V_{\text{rot}}$ and $\sigma_{\text{v}}$ are estimated by comparing PV diagrams extracted from both the model and observed cubes (Fig. 2, bottom panels). The major axis PV diagram is well reproduced by our axisymmetric kinematic model: the residuals are smaller than ~2σ. The width of the [C I] line profiles is largely driven by resolution effects, but it can be used to constrain the intrinsic gas velocity dispersion. Near the [C I] centre, the observed PV diagram shows asymmetric [C I] profiles with a tail of emission towards the sys-
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Multiwavelength overview of PKS 0529-549. **Top left:** ALMA continuum (blue contours and colour scale) overlaid with the ATCA continuum at 18.5 GHz (black contours). The central ALMA component traces dust, while the Eastern component traces synchrotron emission coinciding with the radio lobe. Blue contours are at 0.16, 0.32, 0.64 mJy beam⁻¹. Black contours are at 0.5, 1, 2, 4, 8, 16 mJy beam⁻¹. The ellipse shows the ATCA beam. **Top middle:** total [C I] map. Contours are at 150, 300, and 450 mJy beam⁻¹ km s⁻¹. The ellipse shows the ALMA beam. **Top right:** [C I] velocity field. The dashed line shows the [C I] kinematical PA. **Bottom left:** stellar continuum from SINFONI. The ellipse shows the PSF of the high-resolution SINFONI cube. **Bottom middle:** Hα+[N II] emission (red colour scale) overlaid with the [O III] λ5007 emission (black contours). **Bottom right:** [O III] λ5007 velocity field overlaid with the ATCA image at 18.5 GHz (black contours). The velocity scale is the same as in the top right-hand panel. The dashed line shows the [O III] λ5007 kinematical PA. In all panels, the cross shows the [C I] centre, while the bar to the bottom right corner corresponds to 2 kpc.</figcaption></figure>

temic velocity. These asymmetric profiles are well reproduced by our disc model, indicating that they are probably due to unresolved rotation. This is the typical effect of beam smearing in poorly resolved galaxy discs (e.g. Lelli, Fraternali & Sancisi 2010). Channel maps are also individually inspected (Fig. 3): the [C I] emission progresses along the major axis of the dust component, which is natural if both dust and [C I] lie in a rotating disc. In particular, the estimated values of centre, PA, and *i* are consistent with the dust continuum, providing an independent constraint to the model.

Given the relatively low-resolution and S/N, 3D BAROLO cannot be run blindly leaving all the parameters free. We use $V_{\text{rot}}$ and $\sigma_V$ as free parameters, but keep the others fixed to the values estimated with GALMOD. The fit is performed within the Boolean mask described in Section 2.1. The best-fitting values do not strongly depend on the mask definition, but they do depend on the initial estimates. The best-fitting values of $V_{\text{rot}}$ show a weak dependence and cluster around 310 km s⁻¹. The best-fitting values of $\sigma_V$ show a stronger dependence on the initial estimates, but they are always smaller than ~30 km s⁻¹. Considering that the [C I] cube has low S/N and independent channels with a width of 50 km s⁻¹, the minimum detectable velocity dispersion is ~50/2.35 ≃ 21 km s⁻¹, so we consider $\sigma_V$ ≲ 30 km s⁻¹ as an hard upper limit. These values of $V_{\text{rot}}$ and $\sigma_V$ are in agreement with the previous visual estimates from GALMOD. To estimate errors on each variable, we re-run 3D BAROLO leaving such variable free and keeping all the others

fixed to their best-fitting values. The final parameters are given in Table 1.

One may wonder whether the [C I] velocity gradient may be explained by a gas outflow or a galaxy merger. The outflow scenario implies that (1) the gas is ejected along the major axis of the dust/stellar component rather than the minor axis (see Fig. 2.2), and (2) the PA of the outflow is tilted with respect to the synchrotron lobes, which are generally tracing the AGN ionization cone (e.g. Drouart et al. 2012). Both facts seem unlikely.

The merger origin for the [C I] velocity gradient cannot be definitively ruled out, but we disfavour this scenario, given the high symmetry of the observed *PV* diagram. A kinematic merger model would have twice as many free parameters compared to a disc model and would need to be fine tuned: (1) two merging [C I] components should have roughly the same mass, size, and velocity dispersion, (2) they should be aligned along the major axis of the dust continuum and have the same projected distance from its centroid (see Fig. 3), and (3) they should be well separated in velocity by ~500 km s⁻¹ and lie very close on the sky (≈0.5 arcsec), but they should not overlap at any spatial position after beam convolution, otherwise we should observe double-peaked [C I] profiles (contrary to what is seen in Fig. 2). Naturally, PKS 0529-549 may still be a late-stage merger with an inner rotating disc, as observed in some starburst galaxies at $z \simeq 0$ (e.g. Kregel & Sancisi 2001; Weaver et al. 2018). Here, we simply point out that the observed [C I] velocity gradi-
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top panels: observed [C I] map (**left**), model [C I] map (**middle**), and residual map (**right**). Contours are at 150, 300, and 450 mJy beam⁻¹ km s⁻¹. In the left-hand panel, the dashed line shows the 3σ contour from the model map. In all panels, the cross, the dashed line, and the ellipse show the disc centre, the major axis, and the ALMA beam, respectively. Bottom panels: Major axis PV diagrams obtained from the observed (**left**), model (**middle**), and residual (**right**) cubes. Contours are at −2σ (dashed), 2σ (solid), 3σ, and 4σ where σ = 0.55 mJy beam⁻¹. The vertical and horizontal lines correspond to the dynamical centre and systemic velocity, respectively.</figcaption></figure>

ent is probably due to rotation, allowing the following dynamical analysis.

## 3.3 Mass budget

De Breuck et al. (2010) measured the rest-frame H-band luminosity ($L_H$) and stellar mass ($M_*$) of 70 radio galaxies. For PKS 0529-549, they give $L_H \simeq 4 \times 10^{11}$ L⊙ and $M_* \simeq 3 \times 10^{11}$ M⊙. The stellar mass is derived by fitting the *Spitzer* SED and assuming (i) a formation redshift of 10, (ii) an exponentially declining SFH, and (iii) a Kroupa (2001) initial mass function (IMF). These assumptions are sensible for radio galaxies in a statistical sense (Rocca-Volmerange et al. 2004). For individual galaxies, however, the SFH remains uncertain. PKS 0529-549 is indeed bursting at $z \simeq 2.6$ since the estimated star formation rate (SFR) is ~1000 M⊙ yr⁻¹ (Falkendal et al. in preparation). For reference, galaxies with $M_* \gtrsim 10^{11}$ M⊙ at $z \simeq 2-3$ have SFRs of the order of ~100 M⊙ yr⁻¹ (e.g. Bisigello et al. 2018): PKS 0529-549 lies well above the mean SFR-$M_*$ relation (the so-called star-forming main sequence).

The effect of starbursts in high-z radio galaxies has been studied by Drouart et al. (2016) in 11 objects by modelling their full SED, from ultraviolet to NIR wavelengths. Unfortunately, their sample does not include PKS 0529-549, but it does include 4C 41.17, which is analogous to PKS 0529-549 in several aspects. In particular, the rest-frame optical spectrum of 4C 41.17 is strikingly similar to the one of PKS 0529-549 (Man et al. in preparation). Drouart et al. (2016) finds that 4C 41.17 is experiencing a recent starburst (~30 Myr) that makes up ~92 per cent of the bolometric luminosity

and forms ~44 per cent of the total stellar mass. Using this analogy and considering the large uncertainties, we adopt a conservative errorbar on the stellar mass of PKS 0529-549: $M_* = (3.0 \pm 2.0) \times 10^{11}$ M⊙. This errorbar does not include systematic uncertainty due to the IMF: a Salpeter IMF would increase the stellar mass by ~60 per cent.

The gas mass is given by three main mass components: atomic hydrogen ($M_{H1}$), molecules ($M_{mol}$), and the warm ionized medium ($M_{WIM}$). Unfortunately, H I emission is currently not accessible at high z. However, it is likely that the central parts of massive starburst galaxies are dominated by H$_2$ rather than H I, as observed in nearby interacting galaxies (Kaneko et al. 2017). Considering the observed [C I](2–1) flux of 2.1 Jy km s⁻¹ and the theoretical relations from Papadopoulos, Thi & Viti (2004), we estimate $M_{mol} \simeq 4 \times 10^{10}$ M⊙. This [C I]–to–H$_2$ conversion assumes the [C/H$_2$] abundance of the local starburst galaxy M82 and is uncertain by about a factor of 2 (see also Man et al. in preparation). Finally, the mass of the WIM is clearly negligible: Nesvadba et al. (2017) estimate $2.4 \times 10^9$ M⊙ for a fiducial electron density of 500 cm⁻³. Hence, the total baryonic mass of PKS 0529-549 is $M_{bar} = (3.4 \pm 2.0) \times 10^{11}$ M⊙. The dynamical mass within a radius R is given by

$$M_R = \varepsilon \frac{R V_{\text{rot}}^2}{G}, \qquad (1)$$

where G is the gravitational constant and ε is a parameter of the order of unity that depends on the 3D geometry of the total mass distribution. For simplicity, we adopt ε = 1 (spherical geometry). For local disc galaxies, McGaugh (2005) estimates ε ≃ 0.8, so
<figure><img src="image_5.png" /><figcaption>**Figure 3.** [CI] channel maps from the observed (blue) and model (red) datacubes overlaid on the ALMA continuum (grey scale). Contours are at $-2\sigma$ (dashed), $2\sigma$ (solid), $3\sigma$, and $4\sigma$ where $\sigma = 0.55$ mJy beam$^{-1}$. In the top-left corner, we provide the line-of-sight velocity with respect to $V_{\text{sys}}$. The ALMA beam is shown in the bottom-left corner. The cross corresponds to the dynamical centre.</figcaption></figure>

<table><caption>**Table 1.** Results of the 3D dynamical analysis.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$x_0$ (J2000)</td><td>$05^{\circ}30^{\circ}25.44^{\circ} \pm 0.01^{\circ}$</td></tr><tr><td>$y_0$ (J2000)</td><td>$-54^{\circ}54^{\circ}23^{\prime}2 \pm 0.1^{\prime\prime}$</td></tr><tr><td>$z$ (from $V_{\text{sys}}$)</td><td>$2.570 \pm 0.001$</td></tr><tr><td>PA (°)</td><td>$75 \pm 12$</td></tr><tr><td>$i$ (°)</td><td>$50 \pm 4$</td></tr><tr><td>$V_{\text{rot}}(\text{km s}^{-1})$</td><td>$310 \pm 50$</td></tr><tr><td>$\sigma v(\text{km s}^{-1})$</td><td>$\lesssim 30$</td></tr><tr><td>$M_* (10^{11} M_{\odot})$</td><td>$3.0\text{yy} \pm 2.0^a$</td></tr><tr><td>$M_{\text{mol}} (10^{11} M_{\odot})$</td><td>$0.4 \pm 0.2^b$</td></tr><tr><td>$M_{\text{WIM}} (10^{11} M_{\odot})$</td><td>$0.02 \pm 0.02^c$</td></tr><tr><td>$M_{\text{bar}} (10^{11} M_{\odot})$</td><td>$3.4 \pm 2.0$</td></tr><tr><td>$M_{4\text{kpc}} (10^{11} M_{\odot})$</td><td>$0.9 \pm 0.3$</td></tr><tr><td>$M_{8\text{kpc}} (10^{11} M_{\odot})$</td><td>$1.8 \pm 0.6$</td></tr></tbody></table>

one would infer a slightly smaller dynamical mass. For PKS 0529-549, we measure $V_{\text{rot}} = 310 \pm 50$ km s$^{-1}$ within the deconvolved radius of the [CI] disc ($R \sim 0.45$ arcsec $\sim 4$ kpc), thus $M_{4\text{kpc}} = (0.9 \pm 0.3) \times 10^{11} M_{\odot}$.

The value of $M_{4\text{kpc}}$ is consistent with $M_{\text{bar}}$ within the errors albeit it appears somewhat smaller. It is likely, however, that the stellar component is more extended than the [CI] disc. Unfortunately, the available NIR images have poor spatial resolution, but they tentatively indicate that the stellar component may be twice as extended as the [CI] disc (cf. with fig. 5 of Broderick et al. 2007). If we extrapolate the measured rotation velocity out to 8 kpc, we get $M_{8\text{kpc}} = (1.8 \pm 0.6) \times 10^{11} M_{\odot}$ that is closer to $M_{\text{bar}}$.

Alter- natively, one may question the estimated inclination. If the [CI] disc were oriented more face-on ($i \sim 30^{\circ}$), the rotation velocity would increase to $\sim 475$ km s$^{-1}$ and give $M_{4\text{kpc}} \sim 2 \times 10^{11} M_{\odot}$. These values of $i$ and $V_{\text{rot}}$, however, are difficult to reconcile with the observed [CI] map (Fig. 2, top panels).

In conclusion, the baryonic and dynamical masses of PKS 0529-549 are consistent within the uncertainties. This suggests that PKS 0529-549 has little (if any) dark matter in its central parts. This is comparable to other massive galaxies at similar redshifts (Toft et al. 2012; Wuyts et al. 2016) as well as massive late-type galaxies (Noordermeer et al. 2007; Lelli, McGaugh & Schombert 2016a) and early-type galaxies (Cappellari 2016) at $z \sim 0$. The latter ones are commonly thought to be the descendants of high-$z$ radio galaxies like PKS 0529-549 (e.g. Pentericci et al. 2001).

# 4 DISCUSSION

## 4.1 Comparison with other high-$z$ galaxies

The first IFS surveys at $z \sim 1-3$ suggested that one-third of star-forming galaxies were rotation-dominated systems, whereas the remaining two-thirds were either merging systems or dispersion-dominated objects (Förster Schreiber et al. 2009; Gnerucci et al. 2011). Subsequently, it was acknowledged that beam-smearing effects can play a major role in assessing the intrinsic galaxy dynamics: the most recent IFS surveys now indicate that at least $\sim 80$ per cent of massive star-forming galaxies ($10^9 < M_* / M_{\odot} < 10^{11}$) host rotating gas discs at $z \sim 1-3$ (Wisnioski et al. 2015; Stott et al. 2016).

While rotating gas discs are now thought to be ubiquitous at $z \sim 1-3$, their detailed properties are still debated. The general picture

<footnote>a De Breuck et al. (2010), b Man et al. (in preparation), c Nesvadba et al. (2017).</footnote>
is that high-z discs are more turbulent than their local analogues (Förster Schreiber et al. 2009; Lehnert et al. 2009; Gnerucci et al. 2011): $\sigma_v$ is thought to increase systematically with z, while $V_{rot}/\sigma_v$ decreases (Wisnioski et al. 2015; Stott et al. 2016). Di Teodoro, Fraternali & Miller (2016), however, point out that high-z discs are poorly resolved and beam-smearing effects can lead to systematic overestimates of $\sigma_v$ and underestimates of $V_{rot}$, unless they are properly modelled in 3D. For PKS 0529-549, we find $V_{rot}/\sigma_v \gtrsim 10$, which is comparable to CO and H I discs in spiral galaxies at z = 0 (Leroy et al. 2008; Mogotsi et al. 2016). Thus, the [C I] disc of PKS 0529-549 is not particularly turbulent. This is remarkable considering that PKS 0529-549 is forming stars at $\sim 1000 \text{ M}_\odot \text{ yr}^{-1}$ and hosts a powerful radio-loud AGN, so a large amount of energy should be injected into the inter-stellar medium (e.g. Lehnert et al. 2009).

The gas phase may be important in assessing the degree of turbulence of a galaxy disc. In general, ionized gas has higher $\sigma_v$ than neutral gas due to thermal broadening, but is also more strongly affected by stellar and AGN feedback. Thus, it is possible that some high values of $\sigma_v$ in 'turbulent' discs are actually due to undetected non-circular motions and/or outflows in the ionized gas component. For example, the [O III] velocity field in Fig. 2.2 could be interpreted as a turbulent disc if radio images and/or [C I] kinematics were not available. It is possible, indeed, to fit the [O III] cube with a rotating disc model, following the same procedures described in Section 3.2 for the [C I] cube. This gives $V_{rot} \simeq 100 \text{ km s}^{-1}$ and $\sigma_v \simeq 300 \text{ km s}^{-1}$, leading to $V_{rot}/\sigma_v \simeq 0.3$. Large IFS surveys find that 'main-sequence' galaxies at $z \simeq 2-3$ have $V_{rot}/\sigma_v \simeq 1-4$ when using $H\alpha$ or [O III] lines (Wisnioski et al. 2015). The lower value of PKS 0529-549 corroborates the outflow interpretation for the [O III] line. However, it also raises the question of how much gas outflows may affect the mean value of $V_{rot}/\sigma_v$ in large IFS surveys. Clearly, a multiwavelength approach is needed to have a complete picture of the gas dynamics of high-z galaxies.

## 4.2 The baryonic Tully-Fisher relation

A key constraint for galaxy formation models is represented by the emergence and evolution of the Tully-Fisher (TF) relation (Tully & Fisher 1977). Over the past years, both the stellar-mass TF relation (using $M_*$) and the baryonic TF relation (using $M_{bar} = M_* + M_{gas}$) have been investigated up to $z \simeq 3$ with contradicting results. Some studies find significant evolution (Tiley et al. 2016; Übler et al. 2017), while others do not (Miller et al. 2011, 2012; Di Teodoro et al. 2016; Pelliccia et al. 2017). According to Turner et al. (2017), these discrepancies are due to selection effects since different authors use different criteria to define 'disc-dominated' galaxies. Data quality must also play a major role. Different authors use different techniques to derive rotation velocities, velocity dispersions, and stellar masses. The latter ones can be particularly uncertain when comparing objects at different redshifts due to systematic effects. For example, Turner et al. (2017) find that the normalization of the stellar-mass TF relation varies by $\sim 0.1$ dex from $z=0$ to $z \simeq 3$. This offset is comparable to the uncertainty in the absolute calibration of stellar masses at $z=0$ from Spitzer [3.6] photometry (e.g. Lelli et al. 2016a), so it is hard to tell whether the stellar-mass TF relation actually evolves. The situation is exacerbated when considering the baryonic TF relation: the gas fraction of galaxies ($M_{gas}/M_*$) likely evolves with redshift, but direct estimates of gas masses are rarely available.

In this context, PKS 0529-549 is particularly interesting because we have estimates of both $M_*$ and $M_{mol}$, other than $V_{rot}$ and $\sigma_v$.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The location of PKS 0529-549 on the local BTFR from the SPARC sample (Lelli et al. 2016a,b).</figcaption></figure>

Fig. 4 shows that PKS0529-549 lies on the local baryonic TF relation from the SPARC sample (Lelli et al. 2016a; Lelli, McGaugh & Schombert 2016b). This baryonic TF relation is consistent with previous calibrations (e.g. McGaugh 2005), but has the advantage of using homogeneous photometry at Spitzer [3.6], which greatly reduces the uncertainties in $M_*$. The stellar masses are calculated assuming $M_*/L_{[3.6]} = 0.5 M_\odot/L_\odot$ for all galaxies, as expected from stellar population synthesis models with a Kroupa (2001) IMF (same as for PKS 0529-549). The rotation velocities are measured along the flat part of the rotation curve ($V_{flat}$), which is probed by deep H I observations (see Lelli et al. 2016a,b for details). In the case of PKS 0529-549, $V_{rot}$ is an intensity-weighted estimate over the semimajor axis, since the [C I] emission is resolved with $\sim 2$ beams. Thus, one may wonder whether we are probing $V_{flat}$. Local galaxies with similar masses as PKS 0529-549 typically have rotation curves that peak at very small radii ($R \lesssim 1$ kpc) and decline by $\sim 20-30$ per cent before reaching $V_{flat}$ (Noordermeer et al. 2007). If the rotation curve of PKS 0529-549 has a similar shape, then our intensity-weighted value should not differ from $V_{flat}$ by more than $\sim 20-30$ per cent. In general, the adherence of PKS 0529-549 to the baryonic Tully-Fisher relation seems solid against systematic effects. Clearly, a single object cannot be used to infer general conclusions. However, Fig. 4 suggests that some massive galaxies can be in place and kinematically relaxed at $z \simeq 2.6$, when the Universe was only $\sim 2.5$ Gyr old.

## 5 CONCLUSIONS

We presented ALMA data of the [C I] (2–1) line for the starburst galaxy PKS 0529-549 at $z \simeq 2.6$. With only 5 min of on-source integration time, ALMA spatially resolves the [C I] distribution and
kinematics. We also analysed SINFONI data (Nesvadba et al. 2017), probing the ionized gas kinematics, and ATCA radio continuum images (Broderick et al. 2007), probing synchrotron emission from a central AGN. Our results can be summarized as follows:

(i) Both [C I] and [O III] display regular velocity gradients, but their systemic velocities and kinematic PAs differ by ~300 km s⁻¹ and ~30°, respectively. The [C I] is consistent with a rotating disc, being aligned with both the stellar and dust components, while the [O III] likely traces an outflow, being aligned with two AGN-driven radio lobes.

(ii) The [C I] cube is well reproduced by a 3D disc model with $V_{\text{rot}} \simeq 310 \text{ km s}^{-1}$ and $\sigma_v \lesssim 30 \text{ km s}^{-1}$. This gives $V_{\text{rot}}/\sigma_v \gtrsim 10$ similar to local spiral galaxies, indicating that the [C I] disc of PKS 0529-549 is not particularly turbulent.

(iii) The dynamical mass within 8 kpc is ~1.8 × 10¹¹ M☉. This is comparable to the baryonic mass within the errors, implying that baryons dominate over dark matter in the central parts. This is similar to massive galaxies at $z \simeq 0$.

(iv) PKS 0529-549 lies on the local baryonic Tully-Fisher relation, suggesting that some massive galaxies can be already in place and kinematically relaxed at $z \simeq 2.6$.

This study demonstrates the potential of the [C I] line to trace galaxy dynamics at high-z. It also highlights the importance of multiwavelength observations to properly interpret gas kinematics during the early stages of galaxy formation.

## ACKNOWLEDGEMENTS

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00521.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

## REFERENCES

Arribas S., Colina L., Bellocchi E., Maiolino R., Villar-Martín M., 2014, A&A, 568, A14

Bisigello L., Caputi K. I., Grogin N., Koekemoer A., 2018, A&A, 609, A82

Boomsma R., Oosterloo T. A., Fraternali F., van der Hulst J. M., Sancisi R., 2005, A&A, 431, 65

Broderick J. W., De Breuck C., Hunstead R. W., Seymour N., 2007, MNRAS, 375, 1059

Cappellari M., 2016, ARA&A, 54, 597

Chen C.-C. et al., 2017, ApJ, 846, 108

De Breuck C. et al., 2010, ApJ, 725, 36

De Breuck C. et al., 2014, A&A, 565, A59

Di Teodoro E. M., Fraternali F., 2015, MNRAS, 451, 3021

Di Teodoro E. M., Fraternali F., Miller S. H., 2016, A&A, 594, A77

Drouart G. et al., 2012, A&A, 548, A45

Drouart G., Rocca-Volmerange B., De Breuck C., Fioc M., Lehnert M., Seymour N., Stern D., Vernet J., 2016, A&A, 593, A109

Förster Schreiber N. M. et al., 2009, ApJ, 706, 1364

Gnerucci A. et al., 2011, A&A, 528, A88

Harrison C. M., Alexander D. M., Mullaney J. R., Swinbank A. M., 2014, MNRAS, 441, 3306

Hodge J. A., Carilli C. L., Walter F., de Blok W. J. G., Riechers D., Daddi E., Lentati L., 2012, ApJ, 760, 11

Humphrey A. et al., 2011, MNRAS, 418, 74

Ivison R. J. et al., 2012, MNRAS, 425, 1320

Jones G. C. et al., 2017, ApJ, 850, 180

Kaneko H., Kuno N., Iono D., Tamura Y., Tosaki T., Nakanishi K., Sawada T., 2017, PASJ, 69, 66

Kregel M., Sancisi R., 2001, A&A, 376, 59

Kroupa P., 2001, MNRAS, 322, 231

Lehnert M. D., Nesvadba N. P. H., Le Tiran L., Di Matteo P., van Driel W., Douglas L. S., Chemin L., Bournaud F., 2009, ApJ, 699, 1660

Lelli F. et al., 2015, A&A, 584, A113

Lelli F., Fraternali F., Sancisi R., 2010, A&A, 516, A11

Lelli F., Verheijen M., Fraternali F., 2014, MNRAS, 445, 1694

Lelli F., McGaugh S. S., Schombert J. M., 2016a, AJ, 152, 157

Lelli F., McGaugh S. S., Schombert J. M., 2016b, ApJ, 816, L14

Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2782

Madau P., Dickinson M., 2014, ARA&A, 52, 415

McGaugh S. S., 2005, ApJ, 632, 859

McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, ASP Conf. Ser. Vol. 376, Astronomical Data Analysis Software and Systems XVI. Astron. Soc. Pac., San Francisco, p. 127

Miller S. H., Bundy K., Sullivan M., Ellis R. S., Treu T., 2011, ApJ, 741, 115

Miller S. H., Ellis R. S., Sullivan M., Bundy K., Newman A. B., Treu T., 2012, ApJ, 753, 74

Mogotsi K. M., de Blok W. J. G., Caldú-Primo A., Walter F., Ianjamasi-manana R., Leroy A. K., 2016, AJ, 151, 15

Nesvadba N. P. H., De Breuck C., Lehnert M. D., Best P. N., Collet C., 2017, A&A, 599, A123

Noordermeer E., van der Hulst J. M., Sancisi R., Swaters R. S., van Albada T. S., 2007, MNRAS, 376, 1513

Papadopoulos P. P., Thi W.-F., Viti S., 2004, MNRAS, 351, 147

Pelliccia D., Tresse L., Epinat B., Ilbert O., Scoville N., Amram P., Lemaux B. C., Zamorani G., 2017, A&A, 599, A25

Pentericci L., McCarthy P. J., Röttgering H. J. A., Miley G. K., van Breugel W. J. M., Fosbury R., 2001, ApJS, 135, 63

Planck Collaboration XIII, 2016, A&A, 594, A13

Rocca-Volmerange B., Le Borgne D., De Breuck C., Fioc M., Moy E., 2004, A&A, 415, 931

Shao Y. et al., 2017, ApJ, 845, 138

Sicking F. J., 1997, PhD thesis, Univ. Groningen

Smit R. et al., 2018, Nature, 553, 178

Stott J. P. et al., 2016, MNRAS, 457, 1888

Tacconi L. J. et al., 2013, ApJ, 768, 74

Tiley A. L. et al., 2016, MNRAS, 460, 103

Toft S., Gallazzi A., Zirm A., Wold M., Zibetti S., Grillo C., Man A., 2012, ApJ, 754, 3

Tully R. B., Fisher J. R., 1977, A&A, 54, 661

Turner O. J., Harrison C. M., Cirasuolo M., McLure R. J., Dunlop J., Swin-bank A. M., Tiley A. L., 2017, preprint (arXiv:1711.03604)

Vogelaar M. G. R., Terlouw J. P., 2001, in Harnden F. R., Jr, Primini F. A., Payne H. E., eds, ASP Conf. Ser. Vol. 238, Astronomical Data Analysis Software and Systems X. Astron. Soc. Pac., San Francisco, p. 358

Weaver J. R. et al., 2018, A&A, 614, A32

Wisnioski E. et al., 2015, ApJ, 799, 209

Wuyts S. et al., 2016, ApJ, 831, 149

Zirm A. W. et al., 2005, ApJ, 630, 68

Übler H. et al., 2017, ApJ, 842, 121

Übler H. D. N. et al., 2018, ApJ, 854, L24

This paper has been typeset from a TEX/LATEX file prepared by the author.