© The Authors 2024

1. Introduction

Galaxies are gravitationally bound objects made of luminous and non-luminous matter, and dynamical modelling enables us to map their spatial distributions (Binney & Tremaine 2008). Understanding the intrinsic shape of galaxies is essential in order to prevent biased conclusions when interpreting results derived from such modelling. The shape is quantified by the axis ratios of a 3D ellipsoid: p = b/a and q = c/a, where a > b > c are semi-major, intermediate, and minor axes; a > b = c is considered prolate, a = b > c is considered oblate and, otherwise, the shape is triaxial. Observations of the luminous components of early-type galaxies (ETGs) have revealed that stars show oblate, prolate, or even triaxial distributions (Weijmans et al. 2014; Foster et al. 2016, 2017; Pulsoni et al. 2018; Ene et al. 2018). These studies used photometry and kinematics to constrain the shapes and found a link between the amount of rotational support and the degree of triaxiality, where more rotation favours oblate distributions, but does not exclude triaxial shapes. Furthermore, an oblate-axisymmetric distribution of stars can be embedded within a triaxial dark matter (DM) halo (e.g. see the results from the TNG simulations by Pulsoni et al. 2020). The Milky Way and its satellites are the only galaxies for which we have detailed 6D phase space information for a sufficiently large sample of tracers with which we can reconstruct the shape of the total mass distribution of both luminous and dark matter (e.g. Sesar et al. 2011; Debattista et al. 2013; Posti & Helmi 2019; Iorio & Belokurov 2019; Naidu et al. 2021; Han et al. 2022). Beyond the Local Group, we rely on projected quantities only. Because of the random orientation of galaxies, we observe these intrinsic distributions projected along a line of sight (LOS), which makes the recovery of intrinsic shape very challenging, especially for ETGs with non-regular rotation (Krajnović et al. 2005; van den Bosch & van de Ven 2009).

While great advances have been made in recent years in our understanding of the intrinsic shapes of the central stellar components of ETGs with dynamical modelling (Jin et al. 2020; Santucci et al. 2022; Thater et al. 2023), much less is known about the haloes, where non-luminous matter dominates. In this domain, cosmological simulations have been the source of insight into the intrinsic shapes of the dark matter haloes (e.g. Cole & Lacey 1996; Allgood et al. 2006; Prada et al. 2019; Chua et al. 2019). The most extensively used simulations are those based on the Λ cold dark matter (ΛCDM) cosmology, which agrees best with the observed properties of galaxies. N-body simulations in the ΛCDM framework have shown strong radial dependence on the halo shape, with a preference for prolate haloes in the inner regions r < 0.1 r_vir, and increasingly triaxial haloes in the outer regions (e.g. Bailin & Steinmetz 2005; Hayashi et al. 2007). The inclusion of baryonic physics acts to make the haloes rounder at all radii and the exact prescription of the baryonic feedback impacts the triaxiality of galaxies (e.g. Kazantzidis et al. 2004; Bryan et al. 2013; Tenneti et al. 2015; Chua et al. 2019). Studies of N-body dissipationless simulations found a strong dependence of the halo shape on different cosmological parameters and the dark matter fraction (Allgood et al. 2006). For example, haloes in warm DM (WDM) cosmology are found to be rounder than cold DM (CDM) haloes (e.g. Bullock 2002; Mayer et al. 2002). Directly observed measurements of the dark matter halo shape can thus provide additional constraints for cosmological simulations.

The inner regions of luminous galaxies are typically dominated by baryons, while in their outer regions, dark matter becomes the dominant contributor to the mass distribution (e.g. Deason et al. 2012; Cappellari et al. 2013). In the case of ETGs, the DM fraction at 5 R_e has been found to be 30−80% (Gerhard 2013, and references therein). Harris et al. (2020a) found a median DM fraction of 80−90% for galaxies with M_⋆ > 10¹⁰ M_⊙. For the recovery of galaxy DM distribution and shapes, we need luminous tracers in the regions where DM dominates. Specific features such as stellar streams and polar rings have been used to put observational constraints on the 3D distribution of the dark matter (e.g. Iodice et al. 2003; Varghese et al. 2011; Khoperskov et al. 2014; Pearson et al. 2015; Leung et al. 2021). Weak lensing provides constraints on the projected flattening of DM haloes on a statistical basis (e.g. Robison et al. 2023).

Recovery of the mass distribution and deprojection of the galaxy orientation of gas-poor ETGs is most efficiently achieved through dynamical modelling of the luminous tracers that are distributed throughout the galaxy (e.g. Cappellari 2008). In the inner few R_e, stars can be used as kinematic tracers; but further out, other halo tracers are needed such as globular clusters (GCs) and planetary nebulae (PNe). Globular clusters are old, dense stellar systems populating the haloes of every large galaxy (Brodie & Strader 2006, and references therein). There is empirical evidence that the total mass in GCs is a better tracer of the total halo mass of its host galaxy than its stellar mass (Blakeslee 1997; Spitler & Forbes 2009; Georgiev et al. 2010; Hudson et al. 2014; Harris et al. 2017). This makes GCs suitable tracers of the overall mass of the galaxies. Thanks to their compact sizes and relatively high luminosities, the photometry of individual GCs can be observed in galaxies up to a few 100 Mpc away from us with the Hubble Space Telescope (HST; e.g. Jordán et al. 2009, 2015; Harris 2023) and in nearby groups and clusters from the ground-based surveys (e.g. Taylor et al. 2016). However, spectroscopic observations needed for radial velocities (RVs) are limited to a few nearby groups and clusters (Brodie et al. 2014; Chaturvedi et al. 2022).

NGC 5128 is the closest ETG, it is classified as giant elliptical or S0 (Harris 2010) and has clear evidence of a recent accretion history or possibly a major merger (Wang et al. 2020, and references therein). It has a stellar mass of 1.6 × 10¹¹ M_⊙ (Hui et al. 1995) and is estimated to host up to 2000 GCs (e.g Harris et al. 2006). NGC 5128 hosts a radio source, Centaurus A, and for simplicity we refer to the galaxy as Cen A. Its proximity has enabled detailed studies of the resolved red giant branch (RGB) stars out to more than 25 R_e (140 kpc) (Rejkuba et al. 2022), providing constraints on the density distribution of the stellar halo. Spectroscopic observations of a large population of PNe (Peng et al. 2004a; Walsh et al. 2015) and GCs (Peng et al. 2004b; Woodley et al. 2010a; Hughes et al. 2023) over a vast halo area provide a necessary sample of kinematic tracers for dynamical modelling; we exploit these in the present work. Several literature studies took advantage of these kinematic samples to measure the spherical mass distribution of this galaxy (Peng et al. 2004a; Woodley et al. 2007, 2010a; Hughes et al. 2021a; Dumont et al. 2024). However, observational evidence of the PNe kinematics (Hui et al. 1995; Peng et al. 2004a; Pulsoni et al. 2018) supports strong triaxiality of the Cen A stellar component and the gravitational potential. In this work, we aim to test the validity of the assumption of spherical symmetry of the Cen A DM halo with axisymmetric dynamical modelling.

We use the large photometric and kinematics datasets of GCs recently compiled to determine the shape of the DM halo of NGC 5128. Throughout the paper, we adopt the following properties: distance D = 3.8 Mpc (Harris et al. 2010)¹, stellar effective radius of R_e = 5 arcmin, position angle (PA) of 35° (Dufour et al. 1979), and systemic velocity v_sys = 541 km s⁻¹ (Hui et al. 1995). In Sect. 2 we introduce the photometric and kinematic datasets used in this work. In Sect. 3 we present a careful analysis of both datasets. The aim of the analysis is to determine the population of GCs that is consistent with being dynamically relaxed in the gravitational potential of Cen A. We describe the dynamical modelling tool used in this work in Sect. 4 and the statistical approach to find the best-fit parameters in Sect. 5. We present the results of the modelling in Sect. 6 and discuss our findings in Sect. 7.

2. Data

2.1. Photometric dataset

In this work, we use the photometric catalogue from Taylor et al. (2017, hereafter T17) with multi-band u′g′r′i′z′ photometry based on the observational programme aiming to characterise the baryon structures around Cen A with the Dark Energy Camera (DECam, Taylor et al. 2016). T17 used 500 000 detected sources within ∼140 kpc of Cen A and covering ∼21 deg². To remove fore- and background contamination, the authors compared the distribution of their candidates with the kinematically confirmed GCs from Woodley et al. (2010b) in the (u′−r′)₀–(r′−z′)₀ plane. The GCs are well separated from Milky Way stars at the red end of the colour–colour plane, however, the confusion increases for blue GCs (see their Fig. 2). The final dataset of ∼3200 candidate GCs is therefore still affected by contamination, which we discuss in more detail in Sect. 3.

We excluded GCs that did not have photometry in all of the photometric bands and selected those within 5 < R < 130 arcmin, to remove the impact of the inner polar dust ring and observational footprint. Additionally, we corrected the online catalogue for extinction. We use the online tool IRSA² to obtain the reddening correction E(B − V)_SFD³ at the location of each GC. We combine these with A_b/E(B − V)_SFD values from Schlafly & Finkbeiner (2011, see Table 6) for R_V = 3.1 in relevant SDSS bandpasses (b) to compute the colour excess A_b. Throughout the paper we also correct all the photometry for a distance modulus value μ = 27.88 (Harris et al. 2010) to bring it to the absolute plane. The reddening and distance corrected magnitude is then $b_0=b-\mu -A_{\mathrm{b}}^{\mathrm{SFD}}$.

2.2. Kinematic catalogue

In this work, we use a heterogeneous compilation of literature kinematic measurements of GCs as presented in Hughes et al. (2021b, hereafter H21; see their Sect. 3 and Table 4). Together with the positions and velocities the authors also report their g and r-band photometry. Out of the whole sample of 557 GCs only 482 have both photometric and kinematic information and we use this subset of GCs in the following analysis and dynamical modelling. This is the largest sample of kinematic measurements for Cen A GCs in the literature and the authors carefully matched the different datasets. The 482 GCs span a magnitude range −11 < r₀ < −6 mag and have a median radial extent of 8 arcmin (9 kpc, 1.6 R_e) with 95% of the GCs within 25 arcmin (27 kpc, 5 R_e). The median velocity uncertainty is 33 km s⁻¹.

The spatial distributions of both the kinematic and photometric catalogues used in this work are shown in Fig. 1. The photometric dataset extends much further than the kinematic. The photometric dataset at R > 130 arcmin clearly shows the observational footprint. The kinematic sample is concentrated in the inner 30 arcmin (∼6 R_e), with a few GCs located in the outer halo. It also shows an elongation in the NE–SW direction, with a larger number of GCs in the II. and IV. quadrants compared with I. and III. in the inner halo. The same trend is also seen in the outer halo, but the number of GCs there is significantly lower.

Fig. 1. Spatial distribution of the photometric GC candidates from T17 in small grey points and of the kinematic dataset from H21 shown with larger blue dots. The dashed and solid blue circles mark the radius of 30 arcmin and 130 arcmin respectively. The blue lines delimit the quadrants used in Sect. 3.1 to determine the magnitude-selected subsample of GC candidates from T17 that are consistent with a relaxed axisymmetric distribution. The Roman numerals indicate the number of the respective quadrant and the Arabic numbers below the number of GCs with kinematic measurements in the corresponding quadrant.

3. Selection of the relaxed tracer component

The dynamical models used in this work (see Sect. 4.1) assume that the kinematic tracers are relaxed and fully phase-mixed in the gravitational potential. Additionally, they require an accurate characterisation of their spatial distribution, the tracer density profile.

The presence of shells, streams and a young population of halo stars all support a recent merger scenario of Cen A (e.g. Malin et al. 1983; Rejkuba et al. 2005, 2011; Crnojević et al. 2016). Recently Wang et al. (2020) used an N-body simulation that can reproduce several observed features of the galaxy. The major merger scenario dates the merger back to ∼6 Gyr, with the final fusion of the progenitors ∼2 Gyr ago. Woodley et al. (2010b) found that Cen A’s globular cluster system (GCS) is dominated by old GCs (τ > 8 Gyr) and has evidence of intermediate age (5 < τ < 8 Gyr) and young (τ < 5 Gyr) populations. The metallicities and ages of the GCs and halo stars suggest that few minor mergers took place after the bulk of the galaxy was assembled ∼10 Gyr ago (e.g. Beasley et al. 2008; Woodley et al. 2010b). Such timescales suggest that a population of GCs may not be completely phase-mixed in the Cen A potential and necessitates a careful photometric and kinematic analysis to determine the relaxed population of GCs for the modelling.

In this section, we present our analysis of the photometric GC candidates from T17 to identify the relaxed population. Furthermore, we analysed the kinematic GC catalogue from H21 and used the smooth component of the PN velocity field as determined by Pulsoni et al. (2018) to investigate the possible presence of prominent non-relaxed features. Lastly, we present how we build the tracer density profile used for dynamical modelling.

3.1. Analysis of the photometric dataset

To build up the tracer density profile of the relaxed GC population we followed the two-step approach presented in Veršič et al. (2024). In the first step, we select a subsample of GCs that is complete up to a limiting magnitude and radius and in the second step we fit an analytic function to the complete subsample. A hallmark of a clean and dynamically relaxed sample of GCs is that the magnitude distribution should closely follow a lognormal globular cluster luminosity function (GCLF) up to the completeness limit of the survey and the spatial distribution should be close to point symmetric.

3.1.1. Comparison between observed and theoretical GCLF

Deep observations of GCs in the nearby Universe have revealed that the magnitude distribution of GCs around galaxies follows a nearly universal shape (e.g. Jordán et al. 2006; Rejkuba 2012; Harris et al. 2014). The Gaussian shape of the GCLF is characterised by a peak turnover magnitude (TOM_r) and a spread (σ_r). The values of these parameters mildly depend on the mass of the host galaxy (Villegas et al. 2010) and, for Cen A, we adopt literature values TOM_r = −7.36 and σ_r = 1.2. In the first step of the analysis, we leverage the universality of the GCLF to further remove contaminants from the T17 catalogue of GC candidates.

To remove the effect of the central polar dust ring and the observational footprint of DECAM (as seen in Fig. 1), we additionally limited the sample to 4′≤R ≤ 130 arcmin (4.4 kpc ≤ R ≤ 143 kpc). By comparing the theoretical GCLF and the observed magnitude distribution, we identified 5 different photometric subsamples, which we labelled A–E. The observed GCLF corrected for distance and extinction is compared with the theoretical GCLF. Figure 2 shows the two GCLFs together with the delineation of the different magnitude subsamples. The limiting magnitudes that correspond to the different subsamples are shown in Table 1. Column 2 shows the number of candidate GCs in the corresponding magnitude range and in Col. 3 we computed the estimated number of GCs based on the lognormal GCLF and a total of N_CG = 1450 ± 160 as estimated by H21. Based on the strong deviations from the theoretical GCLF, we refer to all sources with r₀ > −7.9 (D, E) as faint and r₀ ≤ −7.9 (A–C) as bright in the following discussion. The brightest GC bin A is expected to have the lowest contamination fraction, with the contamination expected to increase in fainter magnitude bins; a trend that is observed in subsamples D and E. Subsample D shows the largest excess of observed GCs with respect to the theoretical GCLF, and the completeness limit of the survey is evident in subsample E, with the number of GCs dropping rapidly towards the fainter magnitudes. The faint sources are progressively affected by the survey’s incompleteness and contamination, and therefore we excluded them from the subsequent analysis.

Fig. 2. Comparison between the observed GCLF from the T17 candidate GCs in grey and theoretical GCLF in red. The top panel shows the GCLF and the middle panel the cumulative GCLF, which allows a more precise determination of the magnitude subsamples. Both theoretical GCLFs are scaled, to match the total observed number of GCs in regions A–C. The bottom panel is based on the cumulative GCLF from the middle panel and it shows the difference between the theoretical and observed distributions. The coloured vertical lines delineate different magnitude subsamples labelled with A–E in corresponding colours. The light blue histogram on the top panel shows the magnitude distribution of the kinematic H21 dataset and is arbitrarily scaled for clarity.

On the bright end, we observe less strong deviations between the observed and theoretical GCLF. The brightest magnitude bin A has potential contributions from ultracompact dwarfs (UCDs; Dumont et al. 2022), which have distinct origins from GCs (Haşegan et al. 2005; Mieske et al. 2008; Voggel et al. 2020), and therefore could bias the GC tracer density profile. There is an unexpected overdensity of sources in subsample B that deviates from the theoretical GCLF. Identifying whether this is a real feature or the result of the GC candidate selection process is beyond the scope of this work.

Among all of the samples, subsample C shows the best agreement between the theoretical and observed GC numbers. The values in Cols. 4 and 5 of Table 1 should be taken as guidelines and not as exact contamination fractions, since different studies estimate the total number of GCs in Cen A from 1000 to 2000 GCs (Harris et al. 2006, 2010), which would alter the exact values of the fraction, although the general trends would remain.

3.1.2. Point symmetry of magnitude-selected subsamples

To further investigate whether or not the bright GC subsamples are consistent with being dynamically relaxed, we investigated the azimuthal variation of GCs in eight quadrants as shown in Fig. 1. Guided by the findings of Rejkuba et al. (2022), who reported a transition between the inner and outer halo of RGB stars at R ∼ 30 arcmin, we separately inspected the number of sources in 4 quadrants of the inner (5 < R ≤ 30 arcmin) and outer (30 < R ≤ 130 arcmin) regions. These regions are shown in Fig. 1 as solid and dashed lines for the outer and inner halo regions, respectively.

The azimuthal variation of the fraction of GCs for the subsamples A, B and C in each of the eight quadrants is shown in Fig. 3. A flattened point-symmetric distribution of tracers will show an increased fraction along the major axis of the flattened system. For comparison, we also sampled positions of mock GCs from an elliptical distribution aligned with the stellar PA. The azimuthal variation of such a point-symmetric distribution is shown with a grey dashed line. The coloured lines show the fraction of GCs compared with the subsample as f_GC = N_i, quad./N_i [%], where i stands for A, B, and C, coloured yellow, light blue and purple, respectively. The shaded regions highlight the corresponding 1σ uncertainty region. The orientation of the stellar major axis is along the NE–SW direction and the magnitude-selected subsample C agrees best with the same alignment both in the inner and outer halo regions. This subsample exhibits a higher than 1σ difference in the fraction of the sources along the major and minor axis.

Fig. 3. Fraction of candidate GCs in each of the eight angular bins as shown in Fig. 1. We focus on the subsamples A, B, and C in yellow, light blue, and purple, respectively as defined in Fig. 2 and Table 1. The top panel shows the azimuthal variation of the fraction of GCs in the inner halo (5 < R < 30 arcmin) and the bottom panel the outer halo (30 < R < 130 arcmin). The shaded regions show the associated 1σ uncertainties. The stellar position angle is along NE–SW. To guide the eye, the grey dashed lines show an azimuthal variation for a mock ellipsoidal distribution of points aligned along the Cen A stellar PA.

Subsample A shows deviations from point symmetry both in the inner and outer halo indicating a strong presence of non-relaxed features. Subsample B does show signs of point symmetry in the inner regions, however in the outer regions the behaviour mimics that of subsample A. There has been recent observational evidence supporting the presence of a younger and brighter population of GCs that might not be fully relaxed, which might explain the trends we observe (Voggel et al. 2020; Dumont et al. 2022). However, as is evident from Table 1, there might still be substantial contamination and it is beyond the scope of this work to identify whether the source of the deviation from point symmetry in subsamples A and B is intrinsic or observational. Therefore, we adopt subsample C as the most relaxed sample of candidate GCs from T17 and use it to determine the tracer density profile.

3.2. Planetary nebula versus globular cluster velocity field

Radial velocity measurements in the H21 dataset offer the possibility to investigate whether this sample of kinematic tracers is phase-mixed in the galactic potential. Non-relaxed structures, originating in a recent interaction, appear as overdensities in phase space or perturbations of a point-symmetric velocity field (e.g. Coccato et al. 2013). In this section, we investigated whether the kinematic sample of GCs is consistent with a smooth, point-symmetric velocity field or whether there are substructures that sample a broader range of magnitudes. We investigate both the whole kinematic sample as well as the magnitude selected subsample C. The magnitude distribution of the kinematic sample is shown in Fig. 2 as a light blue histogram on the top panel. Compared with the photometric dataset it samples a narrower range in magnitude, due to a brighter limiting magnitude of the spectroscopic observations.

To investigate the presence of such features, we take advantage of the results from Pulsoni et al. (2018). The authors used the PN catalogue from Walsh et al. (2015) as tracers of the stellar halo kinematics to determine the minimum relaxed velocity and velocity dispersion fields. Non-relaxed structures are higher-level perturbations of the relaxed field and the authors found evidence for deviation from a point-symmetric velocity field at all radii. In this work, we consider the point-symmetric model fitted to the PNe as a good proxy for the mean velocity field of the relaxed population and thus use it for comparison with the GCs.

We carry out this analysis by using the entire sample of GCs with velocity measurements as well as separating them by colour. A more detailed discussion on the colour distribution follows in the next section and here we adopt a colour cut of (g − r)₀ = 0.74 to separate the blue and red GCs. We used this combination of photometric bands because it was available for the whole kinematic sample and a two-component Gaussian mixture model to determine the colour cut.

To compare the GC velocities with the PN velocity field, we used the radial variations of the harmonic expansion coefficients to generate the moments of the smooth PN LOS velocity field at the position of each GC. The coefficients were computed in fixed annular bins and we used a linear interpolator to compute the values of these parameters for arbitrary radii within the extent of the literature measurements. Using higher-order interpolators resulted in unrealistic radial variations of the coefficients.

The point-symmetric PN velocity field from Pulsoni et al. (2018) is plotted in Fig. 4 in the background and the individual GC velocity measurements used in this work are shown with coloured points. The stellar photometric PA and the kinematic PA show clear signs of misalignment which is a signature of triaxiality. However the amplitude of rotation ∼50 km s⁻¹ is much smaller than the velocity dispersion ∼150 km s⁻¹, which supports the axisymmetry assumption in the dynamical modelling approach we use in this work (Sect. 4.1).

Fig. 4. Mean velocity field of the PNe in the background obtained with harmonic expansion (Pulsoni et al. 2018) in comparison with the individual kinematic measurements of GCs. The photometric position angle of the stellar component is shown with a black solid line (Dufour et al. 1979).

For each GC, we use the interpolated moments of the velocity field to compute the local LOS velocity distribution (LOSVD). Then, we determine the probability of observing a GC with a velocity as extreme as measured assuming a local Gaussian LOSVD. The latter assumption is reasonable for discrete halo tracers and has been used in the literature (e.g. Zhu et al. 2016). By investigating the spatial distribution of outliers greater than 1, 2 and 3σ from the PN velocity field we find no discernible overdensity or trend with magnitude or colour out to R ∼ 40 kpc. We find that the individual GCs are consistent with being randomly drawn from a local Gaussian LOSVD, where the mean and σ are informed by the PN velocity field and convolved with the measurement error. The fractions of the GCs with velocities more extreme than 1, 2 and 3σ from the local LOSVD are given in Table 2. We repeated the analysis for the whole kinematic sample and again for the kinematic sample in the magnitude range of subsample C. In both cases, we redid the analysis for the whole subsample as well as the colour-selected subsample. Given that the local Gaussian LOSVD is informed by the point-symmetric velocity and dispersion field our investigation justifies dynamical modelling of the whole kinematics sample of GC based on the assumption of dynamical equilibrium.

These results expand the findings of Peng et al. (2004b) to R ∼ 40 kpc. These authors determined a significant correlation between the GC and PN velocity field, finding that the red GCs correlate more strongly with the velocity field of the PN. They also observed a mild rotational signature both around the major and minor axis for the GCs. A rotation signature was also identified in Hughes et al. (2023), where the velocity field was investigated with a constant kinematic position angle and a constant rotation. The authors similarly find a stronger rotation signature for red GCs than blue GCs further motivating the modelling of the blue and red populations separately.

3.3. Summary of our input observables and the use of GC system subsample C for the dynamical modelling

As illustrated in the previous sections we use photometric subsample C to build the tracer density profile to use in dynamical modelling throughout this work. Among the bright GCs only the subsample C shows a point-symmetric spatial distribution which is expected for a relaxed system. While all bright GCs are expected to be photometrically complete and clean, subsamples A and B nonetheless show deviation from point symmetry. This could be due to the observational biases and presence of non-GC contaminants, such as Milky Way stars, or other observational-related effects, which we cannot ascertain here. Concerning the kinematics, such magnitude selection is not necessary. In Table 2 we demonstrated that the entire kinematic dataset of GCs is consistent with the relaxed point symmetric velocity field from PNe, as velocity measurements are independent of photometric biases. This holds both for the entire sample as well as for the analogously selected kinematic subsample C. Therefore, in the next steps, we use the spatial distributions of T17 subsample C, to constrain the spatial distribution of our tracers, and the kinematics available for the H21 GCs sample.

3.4. Tracer density profile

In the dynamical modelling used in this work (see Sect. 4.1), GCs are used as massless tracers of the gravitational potential. An essential input for the modelling is an accurate description of their spatial distribution. In this section, we build the tracer density profile of the relaxed GCs subsample (the subsample C) in the form of multi-Gaussian expansion (MGE, Emsellem et al. 1994). We also considered separately GCs according to their colour as independent tracer populations as described in the following section.

3.4.1. Colour distribution

The GCSs of all large galaxies show bimodal or multi-modal colour distributions (Brodie & Strader 2006). The different peaks in the colour distribution correspond to different populations of GCs with distinct origins, and therefore to different spatial distributions and kinematics. The latter has been used in the literature in the dynamical modelling of massive ETGs (e.g. Romanowsky et al. 2009; Zhu et al. 2016) and can introduce biases in the modelling results if not accounted for (Veršič et al. 2024). Here we determine the colour separation for the magnitude-selected subsample C.

To separate the GCs into red and blue populations we use Gaussian mixture model algorithm implemented in Python package sklearn. As already pointed out in T17 the largest separation is in (u − z)₀ colour, owing to the large wavelength baseline, and therefore we use this combination of magnitudes. We determine the separation of the 2 components as the colour where the probability of belonging to either component is equal. We use bootstrap with replacements to provide uncertainties and determine the separation to be ${(u-z)}_0=2.{69}_{-0.03}^{+0.04}$. Reassuringly, a similar separating colour was found by Hughes et al. (2021b).

3.4.2. Position angle and flattening

In the following sections, we investigate the position angle and projected flattening (q′ = 1 − e = b/a)⁴ of the red and blue GCs. To measure the PA and flattening, we used the algorithm from mgefit⁵ as summarised in Veršič et al. (2024).

We binned the GCs in radial annuli and computed the PA and q′ in each. We use the different bins to compute the mean parameters and their uncertainties. We explored different binning strategies and found that they have a negligible effect on the final result. For different magnitude regions, we compared the recovered flattening in the 4 radial regions informed by the overdensity in GC candidates observed in T17 (R/R_e ≤ 9 and 9 < R/R_e ≤ 17), as well as the radial regions noted in the RGBs by Rejkuba et al. (2022) (R ≤ 30 kpc and R > 30 kpc).

For subsample C, we observe that the flattening and PA are consistent within all the radial regions. A comparison of the red and blue GCs shows that the red GCs prefer a more flattened distribution with $q_{\mathrm{GCS},\mathrm{red}}^{\prime }=0.{72}_{-0.08}^{+0.08}$ and blue a more spherical distribution with $q_{\mathrm{GCS},\mathrm{blue}}^{\prime }=0.{82}_{-0.02}^{+0.01}$. Within the uncertainties, the PAs of the two components agree with the photometric position angle of the main body stellar component of the galaxy at 35°: ${\mathrm{PA}}_{\mathrm{GCS},\mathrm{red}}=46.0_{-10.5}^{+10.5}$ and ${\mathrm{PA}}_{\mathrm{GCS},\mathrm{blue}}=30.2_{-8.0}^{+8.0}$.

To correct for the PA and rotate the sample to a common PA, we use the literature results from the main stellar component of the galaxy PA = 35°. This ensures the potential generated by the stars and the kinematics of the GCs are aligned with the same PA as required by the axisymmetric dynamical models used in this work.

3.4.3. Fitting an analytic tracer density profile

Based on the colour, and as discussed in Sect. 3.4.1, we split the GCs in subsample C into blue and red. To build the tracer density profile that will be an input to the dynamical modelling, we first fit the observed number density profile of the GCs with an analytic function and then decompose it in a series of Gaussians utilising mgefit. Due to the low number of sources in red and blue populations, we find that the most robust way to determine the MGE components is to construct the binned profile, fit an analytic function to that binned profile, and then decompose the fitted function into a series of Gaussians.

Assuming the PA of 35°, we align the major axis of the GCS with that of the main body of the galaxy. To account for the flattened spatial distribution of GCs, we assume the density is constant on confocal ellipses $R_{\text{ell}}^2=x^2+{(y/q^{\prime })}^2$ and use this to make a binned profile with 16 logarithmically spaced bins. Bins with fewer than 3 GCs are removed because of low statistical significance. We explored different binning strategies to ensure the best-fit profile is minimally affected by it.

We found that a sum of a power-law (Σ_P) and Sérsic (Σ_S) profile (Sersic 1968) can best represent the data:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\mathrm{P}}(R_{\mathrm{ell}})& =A_i\left[{\left(\frac{R_{\mathrm{ell}}}{1\mathrm{arcsec}}\right)}^{-{\alpha }_i}\right],\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\mathrm{S}}(R_{\mathrm{ell}})& ={\mathrm{\Sigma }}_{\mathrm{e},i}exp[-b_n[{\left(\frac{R_{\mathrm{ell}}}{R_{\mathrm{e},i}}\right)}^{1/n_i}-1]],\hfill \end{array}$$(2)

where R_e, i is the effective radius of the ith subpopulation, Σ_e, i is the corresponding number density, n_i the Sérsic index, and we use the approximation of b_n = 2n − 1/3 + 4/405 × 1/n + 46/25515 × 1/n² (Ciotti & Bertin 1999). The power-law profile was computed in units of arcsec, A_i is its scale and α_i is the slope. Σ_S and Σ_P were computed in units of [N_GC arcsec⁻²].

To find the best-fit parameters we use a Gaussian likelihood and uninformative (flat) priors on the modelled parameters (Σ_e, i,  R_e, i,  n_i,  A_i and α_i), the priors are given in Appendix A. The slope and Sérsic index were sampled in linear space and Σ_e, i,  R_e, i, and A_i in logarithmic space as they can span multiple orders of magnitude and because this transformation results in more efficient sampling of the parameter space.

We sample the posterior, find the maximum likelihood region and explore the parameter space with EMCEE, an affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013). The posteriors on each of the 5 free parameters are unimodal with the weakest constraints on the Sérsic index. The chains were allowed to evolve for 12 000 steps to ensure convergence of the best-fit parameters; the first 1000 were discarded as burn-in.

The best-fit parameters for different subpopulations are shown in Table 3. The Sérsic profile can better capture the behaviour of the density profile in the inner regions and the single power-law in the outer regions. Because of the choice of the fitting function, the interpretation of the individual best-fit parameters is not straightforward, especially in reference to literature studies. We note that the blue population shows a slightly shallower Sérsic and power-law slope and a larger R_e compared to the red population, which is consistent with the literature.

The resulting profile for the blue and red C subpopulations are shown as a solid black line in Fig. 5, with dark and light shaded regions corresponding to 1 and 3σ uncertainties. The top panel shows the blue GC data as blue points and the associated best-fit model. The bottom shows the same for the red population. Arbitrarily scaled dashed and dotted profiles show the literature results.

Fig. 5. Tracer density distribution for subsample C of GCs. The top figure shows the blue and the bottom the red population. The blue and red points with associated uncertainties indicate the binned profile computed for the respective subpopulation GCs. A solid black line shows the best-fit model to the data, and dark and light grey shadings show 1 and 3σ uncertainties computed as the 16th and 99.7th percentiles from the posterior. The dotted blue and red lines show the best-fit double power-law from T17 for each population. H21 used different datasets to determine three different stages of the clean photometric datasets. They refer to the best sample as golden (Au) and the next best sample as the silver (Ag) sample. The dashed yellow line shows the best-fit model of the golden sample from H21 and the dashed grey line the silver and golden combined. The two literature results are arbitrarily scaled for comparison.

The dashed H21 results were computed without the colour selection. Our profile in the inner regions (R ≤ 30 arcmin) is steeper than their silver sample (Au+Ag) for both the blue and red populations. The blue population profile shows a slope more consistent with their higher purity (gold) sample, which is shown as a dashed orange line. In the outer regions (R > 30 arcmin) our best-fit profile drops off faster than the literature results. The latter might not be surprising as they impose a constant background.

The dotted lines show the broken double power-law fits presented in T17. In the inner region, we find a consistent slope for the red GCs and a steeper slope for the blue GCs. While in the outer region, our best-fit profile shows a steeper slope compared with T17 profiles.

In general, our profiles lie within the literature results and show the main characteristics that have been observed for the GCS of Cen A. The differences stem from the different methods of selecting a high-purity and high-completeness sample of GCs and different fitting functions.

3.4.4. Multi Gaussian expansion decomposition

In the last step, we constructed an MGE representation of the number density profiles separately for the blue and red GCs of subsample C (as defined in Fig. 2). Such a representation is necessary for the dynamical modelling used in this work. We drew samples from the posterior of the parameters from the tracer density fit and used them to evaluate the profile at each radius. We used the mgefit implementation of the MGE algorithm (Cappellari 2002) on such profiles that mimic observational noise following the procedure in Veršič et al. (2024). Lastly, we also accounted for the flattening of both the blue and red GC systems. MGE components for blue and red GCs used in this work are given in Appendix A.

4. Dynamical modelling with GCs

To model the GCs, we use the axisymmetric Jeans Anisotropic MGE (JAM) formalism introduced by Cappellari (2008). In this section, we introduce the models and then describe how we use them to reproduce the GC kinematics and spatial distribution, and thus derive the mass distribution of Cen A. The JAM formalism has been extensively used to model ETGs with stars, gas, PNe, and GCs as tracers (e.g. Cappellari et al. 2013; Poci et al. 2017; Bellstedt et al. 2018).

4.1. C implementation of the JAM models

The JAM models are available assuming both spherical symmetry and cylindrical (R, z, ϕ) axisymmetry. As we need to reproduce the flattening of the tracer density distribution and the stellar and dark matter distributions, we use the axisymmetric models here. We use the C implementation of the JAM models (CJAM, Watkins et al. 2013).

In principle, these models have two free angular parameters that describe the orientation of the galaxy: the PA and the inclination angle i. The PA defines the location of the photometric major axis. In Sect. 3.4.2, we discussed the PA, and we henceforth keep this fixed. The main rotation of the galaxy is along the photometric major axis, which supports our choice of symmetry axis coinciding with the minor axis of the stellar and GC distribution.

The inclination angle is used to connect the intrinsic 3D quantities from the model with the observables (density profiles and velocity components) projected on the sky, where i = 0° for a face-on system and i = 90° for an edge-on system. We keep this parameter free in some of our models.

Both the tracer density and mass density distributions are incorporated into the model through MGE components. By scaling the stellar surface brightness with mass-to-light ratio M_⋆/L_B⁶ we derive the stellar mass density and add it together with the assumed dark matter density to define the gravitational potential:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{tot}}=(M_{\star }/L_B)\times {\nu }_B+{\rho }_{\mathrm{DM}},\end{array}$$(3)

where ν_B is the MGE representation of the deprojected luminosity density and ρ_DM MGE representation of the DM component.

The orbital properties of the kinematic tracers are incorporated through the anisotropy parameter β_z and rotation parameter κ. CJAM assumes a cylindrically aligned velocity ellipsoid and that the anisotropy is defined in the meridional plane as:

$$\begin{array}{c}\hfill {\beta }_z=1-\frac{\overline{v_z^2}}{\overline{v_R^2}},\end{array}$$(4)

where $\overline{v_z^2}$ is the second velocity moment along the z-direction and $\overline{v_R^2}$ along the cylindrical radial R direction. The rotation parameter modulates the ratio of ordered to disordered motion:

$$\begin{array}{c}\hfill \kappa =\frac{\overline{v_{\varphi }}}{\sqrt{\overline{v_{\varphi }^2}-\overline{v_R^2}}},\end{array}$$(5)

where $\overline{v_{\varphi }}$, $\overline{v_{\varphi }^2}$ are first and second order moments in the tangential coordinate ϕ.

4.2. Stellar surface density

To construct the stellar surface density profile, we combine the stellar surface brightness in the inner 1 R_e from Dufour et al. (1979) with the appropriately scaled halo RGB number counts from Rejkuba et al. (2022).

Accounting for the inner dust disc the stellar surface brightness in the range 70 arcsec ≤ R ≤ 255 arcsec (0.2−0.8 R_e) is best represented by the de Vaucouleurs profile in the B-band magnitude (Dufour et al. 1979, Eq. (12)):

$$\begin{array}{c}\hfill I_B=8.32[{\left(\frac{R}{R_{\mathrm{e},B}}\right)}^{1/4}-1]+22.94[{\mathrm{mag}\mathrm{arcsec}}^{-2}],\end{array}$$(6)

where I_B is in mag arcsec⁻² and R_e, B = 305 arcsec is the effective radius in the B band. We adopt a position angle of 35° East of North and a flattening of q_{⋆,  main} = 0.93 from Dufour et al. (1979).

Assuming the distance modulus of μ = 27.88 mag and absolute solar magnitude M_⊙, B = 5.48 (Binney & Tremaine 2008), the B-band magnitude profile is converted to surface brightness profile [L_⊙ arcsec⁻²]:

$$\begin{array}{c}\hfill L_B={10}^{-0.4(I_B-(m-M)-M_{\odot ,B})}[L_{\odot }{\mathrm{arcsec}}^{-2}].\end{array}$$(7)

Rejkuba et al. (2022) used HST images to study the properties of resolved stellar halo of Cen A out to 30 R_e. Constructing deep colour–magnitude diagrams (CMDs) that contained individual bright RGB stars enabled the authors to remove fore- and background contamination and construct the stellar number density profile from 28 individual HST pointings across the galaxy. Using the RGB number counts in the range 0.23 < R/R_e < 30 they find the tightest fit to the number density profile of RGB stars by fitting both a single power-law and for de Vaucouleurs profile with an increased flattening in the outer stellar halo corresponding to $q_{\star ,\hspace{0.17em}\text{halo}}^{\prime }=0.54\pm 0.02$. We adopt the best fit de Vaucouleurs profile from Rejkuba et al. (2022):

$$\begin{array}{c}\hfill log{\varphi }_{\mathrm{RGB}}=7.685-3.008{(a/R_{\mathrm{e}})}^{1/4},\end{array}$$(8)

where a is the semimajor axis position. The authors convert the number counts to luminosity and provide a shift to match the inner stellar surface brightness profile from Dufour et al. (1979), thereby extending the inner stellar surface brightness out to 30 R_e. By integrating the combined functions, we recover the total luminosity of L_B ∼ 0.3 × 10¹¹ L_⊙ corresponding to M_B = −20.7 mag, consistent with B_T = 7.84 mag (M_B = −20.6 mag), the measurements in the Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al. 1991).

To obtain the stellar surface brightness we combine the profile from Dufour et al. (1979) and Rejkuba et al. (2022). To account for the difference in the flattening when making the MGE of the galaxy’s baryonic component, we generate a 2D image using Eqs. (7) and (8) in their respective radial regions. We do not have information on the radial variation of the flattening, so we kept q′ fixed to 0.93 within 0.8 R_e (the inner region) and 0.54 in the outer region beyond 1.5 R_e. In the regions defined by both profiles, we assign the mean luminosity between both models. We decompose the image with the mgefit and the resulting components are shown in Table 4. The profiles along the major and minor axes together with the MGE components are shown in Fig. 6. The blue points show the observed stellar surface density along the major and minor axis in the left and right panels, respectively. The individual Gaussian components are shown with faint coloured lines and the summed profile is shown with a solid red line. The discontinuities along the minor axis reflect the transition region where both profiles are evaluated as a result of the different flattenings of the main stellar body and the halo RGBs.

Fig. 6. Combined stellar surface brightness profile of Cen A from Dufour et al. (1979) and Rejkuba et al. (2022) and its MGE decomposition. The left panel shows the profile along the semimajor axis and the right along the semiminor axis. The blue points show the observed 2D profile we generated with Eq. (6) in the region highlighted by the grey solid line (D+79 fitting region) and Eq. (8) in the outer regions highlighted with the grey dashed line (R+22 fitting region). The red line shows the resulting MGE fit to the profile and the faint coloured lines are the individual Gaussian components. The blue vertical bars at the bottom of the figure show the positions of the kinematic tracers.

4.3. Dark matter

For the dark matter potential, we assume a Navarro, Frenk and White profile (NFW, Navarro et al. 1997) in elliptical coordinates m: $m^2=x^2+y^2+z^2/q_{\text{DM}}^2$, where q_DM is the intrinsic flattening of the dark matter profile. It is important to emphasise that we fix the axis of symmetry to the stellar component. Therefore, the oblate and prolate shapes of dark matter haloes have the same symmetry axis in our modelling approach. Then the DM density is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}(m)=\frac{{\rho }_{\mathrm{s}}}{\frac{m}{r_{\mathrm{s}}}{(1+\frac{m}{r_{\mathrm{s}}})}^2},\end{array}$$(9)

where ρ_s and r_s are the scale density and radius. We parameterise the profile with the total mass M₂₀₀ = (4π/3) $200{\rho }_{\text{crit}}\hspace{0.17em}{r}_{200}^3$, defined as the mass within a sphere of radius r₂₀₀ and an average density 200 times the critical density of the Universe (ρ_crit), and concentration c₂₀₀ = r₂₀₀/r_s. We use the latest results from Planck Collaboration VI (2020) for ${\rho }_{\mathrm{crit}}=(3H_0^2)/(8\pi G)$, where H₀ = 67.4 km s⁻¹ Mpc⁻¹ is the Hubble constant and G is the gravitational constant.

Similar to the stellar density profile, the dark matter profile is parametrised as an MGE to generate the input for CJAM. We decompose the 1D radial profile and account for the intrinsic flattening of the dark matter. Finally, q_DM is projected based on the inclination $q_{\mathrm{DM}}^{\prime }=\sqrt{q_{\mathrm{DM}}^2{sin}^{2}i+{cos}^{2}i}$ as well as the central surface density of the Gaussian components.

5. Discrete likelihood analysis

To compare our models with the data, we follow the discrete likelihood analysis approach described in Watkins et al. (2013), which has been applied to different studies of discrete halo tracers of ETGs and simulated galaxies (e.g. Hughes et al. 2021a; Veršič et al. 2024). CJAM provides projected velocity components of the tracer population given the gravitational potential and orbital parameters (β_z, κ). Through the assumption of a local Gaussian LOSVD we compare the velocity of each of the tracers with model predictions. We use this assumption to construct a likelihood function and constrain best-fit parameters. Given the size of the dataset and its precision, a Gaussian LOSVD is a reasonable assumption. For a tracer k with LOS velocity v_k, the likelihood is given by

$$\begin{array}{c}\hfill ln{\mathcal{L}}_k=-\frac{1}{2}ln\left(2\pi {\sigma }_{\mathrm{LOS},k}^2\right)-\frac{{(v_k-\overline{v_k})}^2}{2{\sigma }_{\mathrm{LOS},k}^2},\end{array}$$(10)

where $\overline{v_k}$ and σ_LOS, k are the mean LOS velocity and LOS velocity dispersion predicted by the model at the position of the kth tracer. The total likelihood of the whole sample is then given by

$$\begin{array}{c}\hfill ln\mathcal{L}={\displaystyle \sum _k}ln{\mathcal{L}}_k.\end{array}$$(11)

Since the measurement uncertainties are Gaussian there is a simple correlation between the likelihood and χ²,

$$\begin{array}{c}\hfill {\chi }^2=-2ln\mathcal{L}.\end{array}$$(12)

In our modelling, we use seven free parameters: three to describe the dark matter (M₂₀₀, c₂₀₀ and q_DM), one to scale the stellar surface density M_⋆/L_B, one to characterise the inclination (inc), and two to describe the orbital properties of the tracer population (β_z and κ).

This seven-dimensional parameter space is explored by evaluating CJAM models on a fixed grid. Each grid point represents a combination of model parameters. Due to the computational time intensity of the CJAM models, we computed the likelihood and χ² statistics on this grid instead of full posterior sampling. We carried out the analysis in two steps. First, to find the rough location of the global minimum, we varied only one parameter, keeping the others fixed. Second, to find the best-fit parameters as well as to investigate their covariances and quantify uncertainties, we varied two parameters at a time keeping the others fixed. We carried out the analysis separately for the blue and red GCs to get independent constraints on each of the parameters.

We found the best-fit solution through a series of iterations. Not all of the seven parameters were explored in each of the iterations. For each parameter (or parameter pair) that is explored, we evaluate the likelihood on a grid; while one of the parameters (or two) was varied the rest were kept fixed to the initial guess at the corresponding iteration. The initial guess of each successive iteration is determined from the output of the previous iteration based on the likelihood that was also used to determine the convergence of the best-fit solution. As the last step of the analysis, we derive both the best-fit parameters and the uncertainties by combining the constraints from the individual 2D grids.

5.1. Finding an approximate location of the global minimum

At each iteration step n, we compute the χ² surface for each of the parameters θ_i that is varied while keeping the other parameters ${\theta }_j^{(n,\mathrm{fix})}$ (for j ≠ i) fixed. We adopted initial estimates for the two components of the gravitational potential from previous investigations: log M₂₀₀ = 12⁷, c₂₀₀ = 10 (Peng et al. 2004a), and the dynamical (M/L)_B = 3.9 (Hui et al. 1995). Our initial estimate for the inclination was i = 75°, and for the orbital parameters we used β_z = 0 and κ = 0.

We used 20 grid points for each parameter as a compromise between the computational time and the coarseness of the grid. For each free parameter θ_i, we evaluated the CJAM model, calculated the χ²(θ_i), and determined the minimum ${({\chi }^2)}_i^{(n)}$ that corresponds to the solution ${\theta }_i^{(n)}$. We used two criteria to inform the location of the global minimum as shown in Fig. 7. One is related to the difference of the minimum χ² of the two subsequent iterations. The second is related to the difference in the value of the solution ${\theta }_i^{(n)}$ between two subsequent iterations ($\mathrm{\Delta }{\theta }_i={\theta }_i^{(n+1)}-{\theta }_i^{(n)}$). Figure 7 shows a sketch of χ² for two iteration steps and demonstrates these two criteria. Between the steps (n) and (n + 1) we modified parameters ${\theta }_j^{n+1,\mathrm{fix}};(j\ne i)$ and the resulting minimum of χ² statistic is lower than the minimum at the previous step: ${\chi }^2({\theta }_i^{(n+1)})<{\chi }^2({\theta }_i^{(n)})$. This negative difference in the minimal χ² value indicates that parameters ${\theta }_j^{n+1,\mathrm{fix}}$ are closer to the best-fit parameters than ${\theta }_j^{n,\mathrm{fix}}$. Additionally, the value of the parameter at the minimum is different: $|{\theta }_i^{(n+1)}-{\theta }_i^{(n)}|\gg \mathrm{\Delta }{\theta }_i$, meaning that the ${\theta }_i^{(n+1)}$ is closer to the best-fit parameter than ${\theta }_i^{(n)}$. Therefore we considered the solution has converged once Δθ_i was a few grid step sizes and ${({\chi }^2)}_i^{(n+1)}\approx {({\chi }^2)}_i^{(n)}$. Because of the correlations in 7D parameter space, ${\theta }_i^{(n)}$ will depend on the values of all other parameters ${\theta }_j^{(n)}$, j ≠ i.

Fig. 7. Sketch of a variation of χ2 statistics as a function of parameter θi for an iteration step (n) in blue and (n + 1) in orange. The points demonstrate the grid points in θi, with a grid step size of Δθi, where the χ2 is computed. The coloured arrows highlight the value of the parameter at the minimum of the χ2 statistics on the x axis and the absolute value of the statistics on the y axis.

As the result of the above iterative procedure, we identify the following best-fit parameters that are close to the global minimum and use them for the initial estimates for the next step of the analysis. These were: log M₂₀₀ = 12, c₂₀₀ = 10, q_DM = 1.2, i = 85°, (M_⋆/L_B) = 4, ${\beta }_z^{\text{GC},\hspace{0.17em}\text{blue}}=-0.2$, ${\beta }_z^{\text{GC},\hspace{0.17em}\text{red}}=0.1$, κ^{GC,  blue} = 0, κ^GC,  red = 0.

5.2. Two-parameter grid model

Next, we explored the covariances by covarying two parameters (θ_i, θ_j) and keeping the others fixed ({θ}_k, where i ≠ j ≠ k). The summary of the ${\left\{\theta \right\}}_k^{(m,\text{fix})}$, the ranges of parameters (i, j) and number of grid points in each parameter is presented in Table 5. The first column shows the iteration step – m, Cols. (2)–(10) contain the values of the ${\theta }_k^{m,\text{fix}}$ in the top row and grid range for (θ_i, θ_j) in the bottom row of each iteration. The last column shows the number of linearly spaced grid points for each θ_i and θ_j. If the covariance is not explored for a given parameter at a given iteration step, then the range row is left empty (–).

In iterations 1−6, we set the rotation of both red and blue GCs to 0 (κ = 0) and we introduce rotation in iterations 7 and 8. This is motivated by an increased computational time when κ^GC ≠ 0 and a recent study that showed that β_z and κ have negligible impact on the enclosed mass measurements from GCS with similar methodology (Veršič et al. 2024).

In iteration step 5, we introduce the literature mass-concentration relation (Correa et al. 2015) using the Python package COMMAH to assign c₂₀₀ given M₂₀₀ and the cosmological parameters from Planck Collaboration XIII (2016). This choice was driven by two observations of the 2D χ² surfaces in the iterations 1−4. First, the strong degeneracy of c₂₀₀ with both log M₂₀₀ and M_⋆/L_B and second, the minimum χ²c₂₀₀ is consistent with the literature studies of mass-concentration relation as explained in Appendix B. In iteration step 5, we additionally introduce a constant inclination, i = 90° (Hui et al. 1995). This choice was motivated because the inclination angle was not further constrained within the range 60°–90° by our data in iterations 1−4 (see 2nd column of 1st–4th row in Figs. C.1 and C.2). We explored the potential correlations between the inclination and other modelling parameters and no significant correlations were found with other parameters.

We set two criteria for the convergence of the solution. First, ${\theta }_i^{(m)}\approx {\theta }_i^{(m,\mathrm{fix})}$ for a given parameter i in all two-dimensional slices of the χ² surface. Second, ${\chi }^2({\theta }_i^{(m)},{\theta }_j^{(m)})\approx {\chi }^2({\theta }_k^{(m)},{\theta }_l^{(m)})$ for all (i, j), (k, l) parameter pairs as shown in Fig. 8. Once both criteria were met, we are confident to have found the global minimum. We used only the blue GCs to locate the minimum in this step. The constraints from the red GCs were used only to inform the population-specific parameters: velocity anisotropy β_z and rotation κ. We chose this approach because, in this way, the convergence of the parameters from the red GCs provides an independent verification of our method and the location of the global minimum.

Fig. 8. Evolution of the minimal χ2 and associated values of the parameter at the minimum, when different parameters were covarying for log M200. Each colour highlights the iteration step and the symbols correspond to the different parameters covarying in the 2D grid. In iteration 7, log M200 was kept constant as seen in Table 5. The top panel shows the evolution of the minimal χ2 value for the individual 2D covariation compared with the global minimum of χ2 for all grids and iterations. The bottom panel shows the difference between the value of the parameter at the minimum with 1σ uncertainties and the value of the parameter shown in Table 5 at a given iteration. The grey horizontal line is shown as a guide.

The iterative process for the blue GCs is demonstrated in Fig. 8 with log M₂₀₀ being a representative θ_i component (all of the parameters are shown in Appendix C for both blue and red GCs). The top panel shows the difference between the minimum value of χ² statistics at the iteration step (m) and the global minimum value. The bottom panel shows the difference in the best-fit value of log M₂₀₀ at a given step and the fixed value $log{M}_{200}^{(m,\mathrm{fix})}$. The θ_j components are shown with different symbols and each colour corresponds to a given iteration. On the top panel the convergence trend is evident, with the minimum value of χ² decreasing with each subsequent iteration. Additionally, the spread of the minimum of the χ² is decreasing, suggesting that all the parameter pairs are converging to the same global minimum. A similar trend is seen in the bottom panel as the $log{M}_{200}^{(m)}$ becomes more consistent among the different (i, j) pairs and gets closer to the value assumed at a given iteration $log{M}_{200}^{(m,\mathrm{fix})}$. A comparison of the initial run marked with blue symbols and the final set-up with black symbols shows improved consistency between parameter pairs. Based on this, we conclude that we identified the global minimum in the 5D parameter space.

5.3. Best-fit parameters and uncertainties

For the quantification of the uncertainties, we start by taking the 2D likelihood surface and, by marginalising over one of the dimensions, we come to the marginalised posterior, which is equivalent to the probability. For each of the θ_i, we obtain four independent constraints on the parameter from the blue GCs and another four from red GCs. We use the joint posterior to quantify the best-fit solution and the 16th and 84th percentile uncertainties.

6. Results

In this section, we present the results from the final iteration – iteration 8. The convergence in both the minimum χ² and the values of the parameter at the minimum support the conclusion that the iteration 8 has converged. This is therefore our final iteration. We first focus on the Δχ² surface followed by the combined constraints from independent modelling of the red and blue GCs and finally, we present the results of the modelling.

6.1. Δχ² surface

The Δχ² surface for the five parameters explored in iteration 8 (see Table 5) is shown in Fig. 9 for the constraints from the blue GCs and constraints from red GCs show similar trends. Each panel shows the smoothed Δχ² surface on a grid of parameter pairs. The black contours highlight the 68.3th, 90th, and 95th percentiles, where the 68.3th percentile approximately corresponds to the 1σ region. The blue cross highlights the location of the minimum χ².

Fig. 9. Smoothed Δχ2 surfaces for the different combinations of parameters from iteration 8 from the blue GCs. The black lines highlight the 68th, 90th, and 95th percentile corresponding to Δχ2 levels of 2.3, 4.6, and 6.3, respectively. The minimal χ2 solution is highlighted with the blue cross.

This figure provides useful insight into the correlations between the different parameters. However, combining the constraints on each of the parameters coming from the Δχ² surface is not trivial and instead, we take advantage of the evaluated likelihoods. We show this figure for demonstrative purposes and the black contours help to guide the eye.

There are several degeneracies between the parameters. The parameter best constrained by the data is log M₂₀₀, which has tight constraints and shows minimal degeneracy with the parameters of the velocity distribution of the GCs – β_z and κ. The total dark matter mass is strongly anticorrelated with the stellar M_⋆/L_B, in such a way that models with lower log M₂₀₀ prefer higher M_⋆/L_B. The flattening of the dark matter halo is positively correlated with log M₂₀₀, with more massive haloes preferring stronger elongation of the isodensity along the minor axis of the stellar and GC distribution of NGC 5128.

The flattening of the DM halo depends very weakly on both the virial dark matter mass and M_⋆/L_B. The strongest constraints on this parameter come from the properties of the GC velocity distribution. As can be seen in the second column, κ provides the strongest constraints on the flattening q_DM, disfavouring an oblate or spherical halo. The velocity anisotropy provides somewhat weaker constraints. The parameters of the velocity distribution are not correlated.

6.2. Combining constraints from red and blue GCs

Following the methodology outlined in Sect. 5.3, we compute the marginalised probability for each of the parameters θ_i. Combining the probabilities we arrive at the joint posterior distribution for all free parameters. The final results are shown in Fig. 10. The red and blue lines correspond to the constraints from the corresponding population of GCs, and the black lines show the joint posterior. The faint lines show the probability distributions coming from the individual parameters.

Fig. 10. Posterior distributions of the individual parameters for the final iteration. Red and blue coloured lines correspond to the constraints from the corresponding GC subpopulation and black shows the joint posterior. Thin lines are the marginalised posteriors from the different parameter pairs shown in Fig. 9 and the thick lines are the joint probabilities from all of the parameter pairs. The parameters in the bottom row are intrinsic to the two GC populations and we do not combine them.

The top row shows the parameters of the stellar and dark matter distributions, which should be consistent between both GC populations. The bottom row shows population-specific parameters. The median, 16th and 84th percentiles of the bold joint probabilities are presented in Table 6.

The total dark matter mass is best constrained by the data, with very consistent posteriors coming from both GC populations. The flattening of the dark matter halo in the top-middle panel shows a broader posterior, with the red population slightly favouring a more prolate halo. The top-right panel shows that the stellar M_⋆/L_B results differ slightly in that the red population favours lower values than the blue population. However, these differences are within the 1σ uncertainties.

The bottom-left panel shows that the red population prefers an isotropic to mild radial velocity anisotropy, while the blue population strongly disfavours isotropy, showing a preference for negative cylindrical anisotropy. The bottom-right panel shows strong consistency in the rotation parameters of both GC populations, with the red GCs preferring a slightly stronger rotation.

6.3. Best-fit model

Now that we have established the best-fit parameters, we can compare the resulting velocity fields of the GCs with that of the best-fitting model. Figure 11 shows the velocity field of the blue GCs in the left column and the velocity dispersion field in the right column. The first row shows the smooth velocity and dispersion fields of the observed GCs, the second row shows the model predictions, the third shows the residuals, and the final row the histogram of the residuals. The smoothing parameters are the same as discussed in Sect. 3.2.

Fig. 11. Comparison between the smooth velocity field of blue GCs and the model. The left column shows the mean velocity and the right column velocity dispersion. The first row shows the data and the second row the CJAM model evaluation for the best-fit values as given in Table 6. The third row shows the residuals between the smooth data and the model and the histogram of the residuals is shown in the bottom row.

In the top-left panel, both major and minor axis rotation is strongly visible, indicative of a triaxial potential. Axisymmetric modelling cannot reproduce this feature as can be seen in the panel below. Instead, the model is sensitive to the average rotation between the major and minor axis rotation and compensates by lowering the rotational amplitude. This can be seen as the negative residuals close to the minor axis, where the axisymmetric rotation has the opposite sign, and underestimation close to the major axis, where the residuals are positive.

The right column shows the velocity dispersion maps. The smoothed data shows signs of a drop in the velocity dispersion along the minor axis, which is even more prominently visible in the model panel below. The velocity dispersion residuals, which are of the order of ∼25%, do not show coherent patterns, unlike in the left panel, and their distribution is much lower compared with the velocity field residuals.

Both top panels show signs of clustering, which are partially driven by the smoothing procedure adopted in this work. As the smooth velocity fields are only used for demonstration purposes and were not used in the model, a full investigation is beyond the scope of this work.

7. Discussion

7.1. Enclosed mass

Based on the best-fit model, the parameters of the NFW distribution are: $M_{200}=1.{86}_{-0.69}^{+1.61}\times {10}^{12}{M}_{\odot }$, $c_{200}=8.{58}_{-0.41}^{+0.35}$, $r_{\mathrm{s}}=29.{94}_{-5.21}^{+6.95}$ kpc, and $r_{200}=257.{11}_{-36.23}^{+44.10}$ kpc. Figure 12 shows the comparison between the constraints for the enclosed mass from this work and literature studies. The lines show results from this work with the stellar mass in grey, dark matter mass in red and combined baryonic and non-baryonic in black. The shaded regions span the 16th to 84th percentiles.

Fig. 12. Enclosed mass measurements from this work and literature studies. Measurements from this work are shown with lines and literature studies with symbols. The grey line shows the stellar enclosed mass, with M⋆/LB determined in this work. The red line with the associated shaded 1σ uncertainty region shows the enclosed dark matter mass integrated in spherical shells with rs the scale radius of an NFW halo. The black line shows the total mass. Circles show literature measurements of the total mass with GCs as dynamical tracers, stars show PN-based measurements, and the square shows measurements from X-ray gas from Kraft et al. (2003). The green point shows the measurement with the whole sample of GCs from Woodley et al. (2007), the yellow point shows the most robust measurement from Woodley et al. (2010a), and the light blue point shows the measurement from Hughes et al. (2023). Blue and red circles show the enclosed mass estimates from blue and red GCs respectively from Hughes et al. (2023). The stars show constraints from PNe, with yellow showing estimate from Hui et al. (1995), light blue from Peng et al. (2004a), with correction from Woodley et al. (2007), and the measurement from the latter study is shown in purple. The location of the red and blue GCs used in this work is shown in the bottom panel.

Our enclosed mass measurement is on average lower than the literature estimates from other dynamical modelling of GCs. The origin of the discrepancy could be the difference in the slope of the tracer density profile. Woodley et al. (2007) reanalysed the enclosed mass from Peng et al. (2004a) and found a factor of 2 higher mass, when using the deprojected 3D tracer density profile instead of 2D, used in the original study. The slope used in Woodley et al. (2007) is steeper than that used in our work.

Woodley et al. (2007) modelled the pressure- and rotationally supported masses of PNe and GCs separately, making no separation based on GC colour. To avoid biases in the inner regions obscured by the dust they excluded R < 5 arcmin. They adopted a distance of 3.9 Mpc to Cen A and we account for the different distances between our and the literature study. Later, Woodley et al. (2010a) applied a tracer mass estimator (TME, Evans et al. 2003) to the whole sample of GCs to estimate the enclosed mass within 20 arcmin. They found the total mass of (5.9 ± 2.0)×10¹¹ M_⊙, which agrees within 2σ with our result at 20 arcmin. They excluded eight outliers based on their projected velocity and radius (v²R), which could inflate the mass estimate and their robust measurements exclude GCs within 1 R_e (R < 5 arcmin). Relaxing the latter criteria results in increased mass estimate at all radii almost by a factor of 2, hinting at the sensitivity of this method to outliers.

Hughes et al. (2023) used the most spatially extended sample of GCs to provide enclosed mass estimates separately for the blue and red GCs. Our measurements agree within 2σ, but are systematically lower. They also used a TME (corrected for rotational support) to estimate the total enclosed mass within a given radius assuming spherical symmetry in the total mass distribution. In their modelling, GCs are spherically distributed with a single power-law profile with an isotropic distribution of orbits. Their measurements are in good agreement with mass estimates from GCs, but are higher than the estimates using PNe. They postulate that the simplified assumption of isotropy in their estimates could be one of the reasons for this, as also commented in Woodley et al. (2007), but the mass-anisotropy degeneracy in spherical symmetry (Binney & Tremaine 2008) prevents them from making a conclusive determination.

Dumont et al. (2024) used literature datasets of PNe, GCs and satellite galaxies to trace the total mass distribution. They employed a similar methodology, with separate luminous and dark matter components and used CJAM. A spherical generalised NFW halo was used to characterise the dark matter halo and in their modelling parameterised with log(ρ_s) and $log({\rho }_{\mathrm{s}}{r}_{\mathrm{s}}^3)$. To characterise the spatial distribution of red and blue GCs they used single power power-law fits from Hughes et al. (2023). In their fiducial model, they used Gaussian priors on the c_vir, concentration at the virial radius⁸. Their fiducial model results in a dark matter mass of $M_{200}=3.6_{-1.0}^{+1.7}\times {10}^{12}{M}_{\odot }$. This value agrees within uncertainty with our measurements, however, it is systematically higher. Some of the differences in our results could arise from the different tracer density profiles adopted, as well as the different data used. The authors carried out robustness tests by constraining the parameters of the gravitational potential with individual tracer groups. This analysis results in very different values of M_vir (see their Table 3) when different tracers are used. The constraints from GCs and satellites agree within the uncertainties, however, when modelling PNe kinematics only they recover a virial mass 6 times lower than their Fiducial model. The implication of this result on the Fiducial model is not discussed.

Pearson et al. (2022) used simulations to model stream formation due to a disruption of a dwarf galaxy. They sampled different masses of the host galaxy and found a lower limit on the mass of M₂₀₀ > 4.7 × 10¹² M_⊙. They compared the visual appearance of the stream observed by Crnojević et al. (2016) and a single LOS velocity along the stream with their simulation results to identify the gravitational potential needed for the observed stream morphology and kinematics. Because of covariances between the mass of the disrupted galaxy, the infall velocity and the line-of-sight position, more LOS velocity measurements along the stream are needed to provide stronger constraints on the mass.

Müller et al. (2022) used the satellite galaxies of the Cen A to determine the dynamical mass within R_max ∼ 800 kpc (the maximal extent of their tracers). Using the mean circular velocity of 27 dwarfs, they estimated a spherical enclosed mass and determined M₂₀₀ = (5.3 ± 3.5)×10¹² M_⊙. This agrees with our measurements within the uncertainties. The dynamical timescales at such large distances, where the satellites are distributed, are longer than the Hubble time and non-virialised motion could inflate the enclosed mass estimate.

The best-fit M_⋆/L_B from this work is lower than that of Hui et al. (1995) and Mathieu et al. (1996), which could be attributed to the different treatment of the DM component in the modelling. Cappellari et al. (2009) modelled the central region of Cen A to determine the mass of the central BH and found M_⋆/L_K = 0.65 ± 0.15.

7.2. Limitations of the modelling assumptions

In the context of CJAM, the flattening of q_DM > 1 that we recover corresponds to the elongation of the dark matter distribution along the minor axis of the stellar and GC distribution. Literature investigations (Hui et al. 1995; Mathieu et al. 1996; Peng et al. 2004a) have identified clear signs of a triaxial halo and, in this section, we compare the modelling results from a prolate and oblate halo to identify the features in the observing plane that favour q_DM > 1.

Figure 13 compares the projected velocity moments from the minimum χ² solution for our fiducial result q_DM = 1.45 in blue and an oblate case with q_DM = 0.7 in orange. For the latter, we recomputed the other parameters from the 2D χ² grids shown in Fig. 9. We chose a value of q_DM = 0.7 as an intermediate flattening between the outer stellar halo flattening of 0.54 and the flattening of the GCS (0.72 for red and 0.82 for blue GCs). Here we present simplified velocity moments along the major and minor axis only. A more detailed analysis of the 2D velocity moments is presented in Appendix D.

Fig. 13. Comparing best-fit prolate (blue) and oblate (orange) CJAM model results. The top panel shows the mean LOS velocity along the major axis. The bottom panel shows the velocity dispersion along the major axis as solid lines and along the minor axis as dashed lines. The LOS velocity is shown symmetrically along the major and minor axes, to highlight the asymmetric rotational signature. However, the velocity dispersion on the bottom panel is shown for positive values of the axes due to the symmetry of σLOS.

The mean velocity along the major axis is higher in the oblate solution, however, the difference is smaller than the typical measurement uncertainty. Beyond 10 arcmin the velocity dispersion profiles of the oblate and prolate DM haloes are very similar for both major and minor axes, where the majority of the kinematic tracers lie. To match the velocity dispersion at these radii, the resulting M₂₀₀ and M_⋆/L_B estimates for the oblate halo are lower compared with the prolate halo. This anti-correlation between the flattening and mass can be seen in Fig. 9. A lower total mass results in a central value of the velocity dispersion that is lower for the oblate solution (∼150 km s⁻¹), compared with the prolate solution (∼165 km s⁻¹). The latter value is in better agreement with the observed central stellar root-mean-square velocity of ∼170 km s⁻¹ in the inner 0.2 arcsec of the Cen A centre from Cappellari et al. (2009).

A drop in the velocity dispersion along the minor axis is also seen in the smooth velocity field in Fig. 11. Such a drop is consistent with the projected dispersion map of a triaxial system with short- and long-axis tube orbits as seen in van den Bosch & van de Ven (2009). The full details of how measured flattening constrained with axisymmetric modelling can be mapped to the intrinsic shape and orientation of the DM halo is a topic for a future study.

7.3. Extremely misaligned dark matter haloes in cosmological simulations

Within the axisymmetric modelling approach used in this work, the dark matter halo shape with q_DM > 1 corresponds to its elongation along the minor axis of the stellar distribution. This raises a question of how frequent are such misalignments between the minor axes of the stellar and DM profiles. For example, in the Milky Way evidence disfavours a configuration in which the minor axis of the disc and DM halo are parallel (e.g. Law et al. 2009; Debattista et al. 2013; Han et al. 2023). In the Auriga analogues of the Milky Way, Prada et al. (2019, see their Fig. 6) found that the minor axes of the DM haloes and angular momentum vector of the stellar discs are aligned for most of their 30 haloes, with 2 haloes showing misalignment > 70° at 1/16 R₂₀₀.

In this section, we conduct a qualitative investigation to examine the frequency of highly misaligned distributions in two cosmological simulations. This analysis aims to enhance our understanding of the results by employing radial regimes analogous to the sensitivity observed in the Cen A study. For this purpose, we use IllustrisTNG (Nelson et al. 2019, and references therein) and Magneticum Pathfinder⁹ simulations to find analogues of Cen A with a central oblate stellar component hosted in a prolate DM halo elongated along the stellar minor axis. We follow a procedure similar to Pulsoni et al. (2021) and Valenzuela et al. (2024) to measure the shapes of the simulated galaxies, corresponding to method S1 from Zemp et al. (2011). This is not carried out through dynamical modelling, but instead by computing the eigenvalues of the inertial tensor calculated using the positions of the simulated particles and their masses within ellipsoidal shells.

In both simulations, we selected galaxies in the mass range of Cen A as a first step: M_vir = 3 × 10¹² − 3 × 10¹³ M_⊙ and M_⋆ = 1 − 5 × 10¹¹ M_⊙. Secondly, we selected those that have oblate stellar distributions embedded in prolate haloes. We then looked at the angle between the stellar minor axis and DM major axis and inspected the frequency of small misalignments between these two axes. We use the definition for triaxiality: T ≡ (1 − q²)/(1 − p²), where q = b/a and p = c/a are the axis ratios between the major (a), intermediate (b) and minor (c) axis. A completely prolate halo has T = 1 (a > b = c), and perfectly oblate halo has T = 0 with (a = b > c). We first selected galaxies with T_⋆ < 1/3 and that are hosted in a halo with T_DM > 2/3. To mimic the configuration of stellar and dark matter in line with Cen A as presented in this work, we compute the shape of the stellar component at 3r_⋆, half (stellar half-mass radius) and the shape of the dark matter halo at 5r_⋆, half. For those, we compute the angle between the minor axis of the stellar distribution and the major axis of the dark matter.

Magneticum. From Box4 (uhr) of the Magneticum Pathfinder simulation suite, which has side length of 68 Mpc, DM particle mass of m_DM = 5.1 × 10⁷ M_⊙ and average stellar particle mass of m_* = 1.8 × 10⁶ M_⊙ (Teklu et al. 2015; Valenzuela & Remus 2024), we find 111 galaxies that match the mass criteria, of which 94 have a modest star formation rate of SFR < 3 M_⊙ yr⁻¹. Of these, 12 galaxies meet the shape criteria as laid out above, of which 10 have modest star formation rates. Only one of those galaxies has a small angle of 17° between the stellar minor axis and the DM major axis, whereas the rest all have such angles above 80°. That galaxy is at the lower end of the mass selection limits with M_⋆ = 1.0 × 10¹¹ M_⊙ and M_vir = 4.6 × 10¹² M_⊙. It has no star formation and shows signs of recent interaction in the stellar halo. However, it should be noted that the stellar component has a relatively spherical shape, such that the orientations of its axes are not very well defined.

IllustrisTNG. We considered the TNG100 simulation run, with a simulation box of 110.7 Mpc on a side, DM particle mass of 7.5 × 10⁶ M_⊙ and mean baryonic particle mass of 1.4 × 10⁶ M_⊙ (Nelson et al. 2019). From the redshift z = 0 snapshot, we identified 455 galaxies in the mass range around Cen A. Of these, 121 have a prolate DM halo (∼27%), but only 17 have an oblate stellar component within 3r_⋆, half. Most of the 121 prolate DM haloes have the major axis close to the plane defined by the intermediate and major axis of the stellar component within 20°. Out of the 17 galaxies selected, only in 1 case there is significant misalignment. This is a clear merger, as shown by the presence of substructures or strong asymmetries in the stellar distribution and velocities. This galaxy has a moderate SFR of 1.8 M_⊙ yr⁻¹. The similarities between the simulated galaxy in the Magneticum simulation and Cen A as well as the consistent merger imprints and low probability of such extreme misalignment found in both simulation suite suggest that it is likely a transient phenomenon. This work opens up possibilities for further studies of how the dark matter halo shape is influenced by the merger history of its host and the environment.

7.4. Connection with the Cen A satellite distribution

Recent studies of Cen A satellites have found evidence of a corotating plane of satellites around Cen A (Müller et al. 2018). In their latest study, Müller et al. (2021) showed that 21 out of 28 dwarf galaxies, which have both distance and radial velocity measurements that make them satellites of Cen A, share a coherent motion within a flattened planar structure. This plane of satellites is almost perpendicular to the dust disc and lies very close to the line connecting Cen A and its nearest massive neighbour M83. Müller et al. (2019) used the accurate distances of Cen A satellites using the tip of the RGB to determine the intrinsic 3D shape of the satellite system. They find the $\mathrm{PA}=60{}_{-23}^{+21}°$ and determine the axis ratios b/a = 0.71 and c/a = 0.41, which results in a triaxiality of T = 0.6, close to prolate shape. They find that the long axis is aligned along the line-of-sight, in contrast with the stellar shape measurements traced by PNe within 25 kpc by Hui et al. (1995) and within 80 kpc by Peng et al. (2004a). The former found T = 0.4 (triaxial, but more oblate) and determined the viewing angles, such that the intermediate axis lies along the line-of-sight. Using these viewing angles, Peng et al. (2004a) found that the kinematics of the PNe at larger distances favour a significantly more prolate shape of T = 0.98.

Müller et al. (2021) explored the frequency of similar flattened structures in TNG100 (Nelson et al. 2019). Such flattened and coherently moving structures are found in only 0.3% of Cen A analogues in dark-matter-only simulations and in 0.2% cases in hydrodynamical simulations. More generally, the distribution of satellite galaxies in the ΛCDM is not expected to be completely random, i.e. isotropic (e.g. Zentner et al. 2005), due to the anisotropic accretion from the cosmic filaments. In this framework, the DM halo of Cen A is expected to be elongated along the major axis of the PNe and the satellite system, which is in contrast with the results of our modelling under axisymmetric assumptions. On the other hand, the triaxial shape of the satellite system with T = 0.6 (Müller et al. 2019) and of PNe with T = 0.98 (Peng et al. 2004a) strongly disfavour spherical or oblate distributions, which is in line with our findings.

8. Summary and conclusions

We carried out axisymmetric dynamical modelling of the blue and red GCs in the Cen A galaxy to measure the total mass and the intrinsic shape of its dark matter halo. The axisymmetric JAM modelling used in this work requires a relaxed system. Given that Cen A shows signatures of past mergers, special care was taken to identify a sample of GCs that show a point symmetry indicative of a relaxed distribution when building the tracer density profile. For the kinematic analysis, we use the entire available sample; we verified that it follows the relaxed point-symmetric velocity field from PNe.

The dataset that constrains the model is composed of the photometric dataset from T17 and the radial velocity sample from H21, which extend out to 115 kpc and include 225 red and 257 blue GCs. We note that 90% of the kinematic dataset is within 40 kpc.

We used axisymmetric JAM modelling together with a discrete Gaussian likelihood analysis to determine the best-fitting parameters of the gravitational potential and orbital properties of the GCs. The gravitational potential is composed of an axisymmetric stellar surface brightness out to 140 kpc (25 R_e), a mass-to-light ratio M_⋆/L_B, and an axisymmetric dark matter distribution. We assumed an NFW form for the dark matter profile, parametrised by M₂₀₀, and used the latest literature mass–concentration relation to determine c₂₀₀. By varying q_DM, we explored the deviation of the DM halo from spherical symmetry. To find the best-fit model, we maximised the Gaussian likelihood. The best-fit parameters were determined in an iterative approach, accounting for the covariances in the modelled parameters. The only parameter that was fixed in the final iteration was inclination, which was taken from Hui et al. (1995), and we verified that it does not correlate with other parameters. We modelled the blue and red GCs separately and used them as two independent samples to verify our approach and determine the highest likelihood solution for the parameters of the gravitational potential. The need to treat them separately was further motivated by the intrinsically different spatial distribution and orbital properties.

We derive a dark matter mass of $M_{200}=1.{86}_{-0.69}^{1.61}\times {10}^{12}{M}_{\odot }$ and $M_{\star }/L_B=2.{98}_{-0.78}^{+0.96}$, resulting in a total stellar mass of M_⋆ ≈ 0.9 × 10¹¹ M_⊙. The total halo mass is lower than values found in previous investigations using GC kinematics, but is in agreement with measurements from PNe as discussed in Sect. 7.1. Our flattening of $q_{\mathrm{DM}}=1.{45}_{-0.53}^{+0.78}$ disfavours an oblate or spherical DM halo. We explored the implications of q_DM > 1 under the axisymmetry assumption of the JAM modelling by examining the frequency of prolate DM haloes that host an oblate stellar distribution in the Magneticum and TNG simulations. We find that in both simulation suites, it is very rare to find DM elongated along the minor axis of the stellar distribution and that this is only seen in galaxies with strong signatures of active or recent mergers. Comparing the smooth velocity field of the GCs, we identified a difference in the velocity dispersion along the major and minor axes that can also arise in intrinsically triaxial potentials (Sect. 7.3). This suggests that GCs, and the modelling carried out in this work, are sensitive to deviations from spherical and cylindrical symmetry.

We find that the orbital properties of blue and red GCs are different. While both populations show mild rotation, we find that the rotation parameter κ is larger for red GCs than for blue GCs. The latter population also shows negative orbital anisotropy in cylindrical coordinates while the red population is consistent with being isotropic. Both results highlight the need for incorporating rotation and anisotropy in future dynamical modelling.

Future work should include more complex modelling that can account for triaxial symmetry and non-relaxed structures, such as streams. To enable such modelling, RV measurements in the outer halo need to include the full extent of the galaxy potential. This will allow us to gain a better understanding of the global potential of the host and its assembly history.

Acknowledgments

We thank the referee for the helpful comments and suggestions that have improved this manuscript. We thank Prashin Jethwa, Ryan Leaman and Laurane Freour for their fruitful discussions and feedback on the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 724857 (Consolidator Grant ArcheoDyn). T.V. acknowledges the studentship support from the European Southern Observatory. L.M.V. acknowledges support by the German Academic Scholarship Foundation (Studienstiftung des deutschen Volkes) and the Marianne-Plehn-Program of the Elite Network of Bavaria. The computational results presented have been achieved (in part) using the Vienna Scientific Cluster (VSC). This work made use of the following software: NUMPY (Harris et al. 2020b), SCIPY (Virtanen et al. 2020), MATPLOTLIB (Hunter 2007), ASTROPY (Astropy Collaboration 2013, 2018, 2022).

References

Appendix A: GC tracer density

For the binning, we first define bin edges based on the closest and outermost GC in the bright photometric population. We use those as limits to generate 16 equidistant bins in log space. We then apply these bin edges to the subsample C and in each bin compute the centre as the median position of GCs in a given bin. The uncertainties are given by the Poisson counting noise in each bin. Bins with fewer than 3 GCs are masked out.

The priors for the fitting of the analytic profile are uniform in the following regions: log R_e ∈ [1, 3] arcsec, α, n ∈ [0.2, 6] and logΣ_S,  logΣ_P ∈ [ − 10, 3] arcsec⁻². The fitted analytic profiles were then decomposed with Multi Gaussian expansion (MGE) into individual Gaussian components, These are listed in Table A.1 for the blue GCs and Table A.2 for the red GCs. The three columns show the parameters of the MGE: the central number density (N_GC pc⁻²), the width along the major axis (σ) in arcsec, and the last is the flattening of the Gaussian component.

The median profile was computed by randomly drawing 1000 samples from the posterior, evaluating the profile given by Eq. 1 and Eq. 2 and computing the median and the percentiles.

Appendix B: Mass – concentration relation

In the 2D grid model set-up there is a strong degeneracy between the total DM mass and the halo concentration. Because of it and in order to reduce the number of free parameters we introduced the literature relation of M₂₀₀ − c₂₀₀ at iteration step 5. The 2D χ² surface in these two parameters is shown in Fig. B.1 with the location of the minimum in the statistics highlighted by the blue cross. The literature M₂₀₀ − c₂₀₀ are shown as black dotted and dashed lines (Neto et al. 2007, for relaxed and non-relaxed haloes respectively), dash-dotted line from Prada et al. (2012) and black and red solid lines from Correa et al. (2015) using cosmological parameters from two different studies: Komatsu et al. (2009, WMAP5) and Planck Collaboration XIII (2016, Planck15). We adopt the latter relation highlighted in red.

Fig. B.1. Δχ2 surface for log M200 − c200 plane in iteration 4 (see Table 5 for a summary of the parameters). The black contour lines show the 68th, 90th and 95th percentile confidence levels, with minimal χ2 solution shown with the blue cross. Dotted, dashed, dot-dashed black lines and solid red and black lines show the different literature M200 − c200. The solid lines show the M200 − c200 computed with COMMAH code from Correa et al. (2015) using the cosmological parameters from Planck Collaboration XIII (2016) in red and the cosmological parameters from Komatsu et al. (2009) in black.

Appendix C: Convergence of the χ² solution

For completeness, we show the convergence of the best-fit solution for the blue GCs in Fig. C.1 and red GCs in Fig. C.2. Table 5 shows how the parameters were varied between the iteration steps.

Fig. C.1. Convergence of the parameters and minimal χ2 for the blue GCs. The symbol shapes and colours correspond to those in Fig. 8 (and are also summarised in the bottom left panel). The panels in the left column show the evolution of the minimal value of the statistics and the right panel the value of the parameter at the minimum χ2.

Fig. C.2. Convergence of the parameters and minimal χ2 for the red GCs. The symbols and colours correspond to those in Fig. 8 and columns are explained in Fig. C.1.

Appendix D: Comparing oblate and prolate models

The results of the comparison between our best fit prolate model with q_DM = 1.38 and an oblate model, where we enforced q_DM = 0.7 are shown in Fig. D.1 for blue and red GCs. To make these plots we evaluated CJAM on a grid of points linearly spaced in R and ϕ for the best fit parameters shown in Table 6. The left column shows the maps of the first moment projected along the line-of-sight and the right column shows the second moment. The first two rows show the model results from the blue GCs best-fit solutions for the prolate (top) and the oblate (bottom) cases same as for Fig. 13. The bottom two rows show the residual of the prolate map w.r.t. the oblate case for the blue GC in the top and red GCs in the bottom. The streaks along the major axis are an artefact of the plotting and not a real feature.

Fig. D.1. Comparison between the first and second moment velocity maps for axisymmetric Jeans modelling in the left and right column, respectively. The first two rows show the best-fit models for blue GCs, with the best fit qDM = 1.38 in the top qDM = 0.7 below. The bottom two rows show the residuals for the blue GCs on the top and red GCs on the bottom. The parameters used to evaluate the model were taken from Table 6 and for qDM = 0.7 we take the other parameters at the minimal χ2 based on the 2D χ2 grids as shown in Fig. 9.

The 2D maps reveal the full complexity of the velocity moments which can be probed by discrete dynamical modelling. The residuals both for blue (in the third row) and red GCs (in the bottom row), show that the best fit oblate model favours a stronger rotational signature. The residuals in the left column are higher for red than blue GCs and they consistently show higher residuals along the semimajor axis. This means that the gradient of the mean velocity from the minor to major axis depends on the flattening of the dark matter halo. Velocity dispersion maps in the right column show a similarly complex signature, with the opposite sign of the residuals. The prolate solution shows on average higher dispersion compared with the best-fit oblate solution. These maps highlight the importance and strength of discrete dynamical modelling, where the position of each GC is used to constrain the 2D velocity distribution, which shows distinct signatures for oblate and prolate solutions.

All Tables

All Figures

Fig. 1. Spatial distribution of the photometric GC candidates from T17 in small grey points and of the kinematic dataset from H21 shown with larger blue dots. The dashed and solid blue circles mark the radius of 30 arcmin and 130 arcmin respectively. The blue lines delimit the quadrants used in Sect. 3.1 to determine the magnitude-selected subsample of GC candidates from T17 that are consistent with a relaxed axisymmetric distribution. The Roman numerals indicate the number of the respective quadrant and the Arabic numbers below the number of GCs with kinematic measurements in the corresponding quadrant.

In the text

Fig. 2. Comparison between the observed GCLF from the T17 candidate GCs in grey and theoretical GCLF in red. The top panel shows the GCLF and the middle panel the cumulative GCLF, which allows a more precise determination of the magnitude subsamples. Both theoretical GCLFs are scaled, to match the total observed number of GCs in regions A–C. The bottom panel is based on the cumulative GCLF from the middle panel and it shows the difference between the theoretical and observed distributions. The coloured vertical lines delineate different magnitude subsamples labelled with A–E in corresponding colours. The light blue histogram on the top panel shows the magnitude distribution of the kinematic H21 dataset and is arbitrarily scaled for clarity.

In the text

Fig. 3. Fraction of candidate GCs in each of the eight angular bins as shown in Fig. 1. We focus on the subsamples A, B, and C in yellow, light blue, and purple, respectively as defined in Fig. 2 and Table 1. The top panel shows the azimuthal variation of the fraction of GCs in the inner halo (5 < R < 30 arcmin) and the bottom panel the outer halo (30 < R < 130 arcmin). The shaded regions show the associated 1σ uncertainties. The stellar position angle is along NE–SW. To guide the eye, the grey dashed lines show an azimuthal variation for a mock ellipsoidal distribution of points aligned along the Cen A stellar PA.

In the text

	Fig. 4. Mean velocity field of the PNe in the background obtained with harmonic expansion (Pulsoni et al. 2018) in comparison with the individual kinematic measurements of GCs. The photometric position angle of the stellar component is shown with a black solid line (Dufour et al. 1979).
In the text

Fig. 5. Tracer density distribution for subsample C of GCs. The top figure shows the blue and the bottom the red population. The blue and red points with associated uncertainties indicate the binned profile computed for the respective subpopulation GCs. A solid black line shows the best-fit model to the data, and dark and light grey shadings show 1 and 3σ uncertainties computed as the 16th and 99.7th percentiles from the posterior. The dotted blue and red lines show the best-fit double power-law from T17 for each population. H21 used different datasets to determine three different stages of the clean photometric datasets. They refer to the best sample as golden (Au) and the next best sample as the silver (Ag) sample. The dashed yellow line shows the best-fit model of the golden sample from H21 and the dashed grey line the silver and golden combined. The two literature results are arbitrarily scaled for comparison.

In the text

Fig. 6. Combined stellar surface brightness profile of Cen A from Dufour et al. (1979) and Rejkuba et al. (2022) and its MGE decomposition. The left panel shows the profile along the semimajor axis and the right along the semiminor axis. The blue points show the observed 2D profile we generated with Eq. (6) in the region highlighted by the grey solid line (D+79 fitting region) and Eq. (8) in the outer regions highlighted with the grey dashed line (R+22 fitting region). The red line shows the resulting MGE fit to the profile and the faint coloured lines are the individual Gaussian components. The blue vertical bars at the bottom of the figure show the positions of the kinematic tracers.

In the text

	Fig. 7. Sketch of a variation of χ2 statistics as a function of parameter θi for an iteration step (n) in blue and (n + 1) in orange. The points demonstrate the grid points in θi, with a grid step size of Δθi, where the χ2 is computed. The coloured arrows highlight the value of the parameter at the minimum of the χ2 statistics on the x axis and the absolute value of the statistics on the y axis.
In the text

Fig. 8. Evolution of the minimal χ2 and associated values of the parameter at the minimum, when different parameters were covarying for log M200. Each colour highlights the iteration step and the symbols correspond to the different parameters covarying in the 2D grid. In iteration 7, log M200 was kept constant as seen in Table 5. The top panel shows the evolution of the minimal χ2 value for the individual 2D covariation compared with the global minimum of χ2 for all grids and iterations. The bottom panel shows the difference between the value of the parameter at the minimum with 1σ uncertainties and the value of the parameter shown in Table 5 at a given iteration. The grey horizontal line is shown as a guide.

In the text

	Fig. 9. Smoothed Δχ2 surfaces for the different combinations of parameters from iteration 8 from the blue GCs. The black lines highlight the 68th, 90th, and 95th percentile corresponding to Δχ2 levels of 2.3, 4.6, and 6.3, respectively. The minimal χ2 solution is highlighted with the blue cross.
In the text

Fig. 10. Posterior distributions of the individual parameters for the final iteration. Red and blue coloured lines correspond to the constraints from the corresponding GC subpopulation and black shows the joint posterior. Thin lines are the marginalised posteriors from the different parameter pairs shown in Fig. 9 and the thick lines are the joint probabilities from all of the parameter pairs. The parameters in the bottom row are intrinsic to the two GC populations and we do not combine them.

In the text

Fig. 11. Comparison between the smooth velocity field of blue GCs and the model. The left column shows the mean velocity and the right column velocity dispersion. The first row shows the data and the second row the CJAM model evaluation for the best-fit values as given in Table 6. The third row shows the residuals between the smooth data and the model and the histogram of the residuals is shown in the bottom row.

In the text

Fig. 12. Enclosed mass measurements from this work and literature studies. Measurements from this work are shown with lines and literature studies with symbols. The grey line shows the stellar enclosed mass, with M⋆/LB determined in this work. The red line with the associated shaded 1σ uncertainty region shows the enclosed dark matter mass integrated in spherical shells with rs the scale radius of an NFW halo. The black line shows the total mass. Circles show literature measurements of the total mass with GCs as dynamical tracers, stars show PN-based measurements, and the square shows measurements from X-ray gas from Kraft et al. (2003). The green point shows the measurement with the whole sample of GCs from Woodley et al. (2007), the yellow point shows the most robust measurement from Woodley et al. (2010a), and the light blue point shows the measurement from Hughes et al. (2023). Blue and red circles show the enclosed mass estimates from blue and red GCs respectively from Hughes et al. (2023). The stars show constraints from PNe, with yellow showing estimate from Hui et al. (1995), light blue from Peng et al. (2004a), with correction from Woodley et al. (2007), and the measurement from the latter study is shown in purple. The location of the red and blue GCs used in this work is shown in the bottom panel.

In the text

Fig. 13. Comparing best-fit prolate (blue) and oblate (orange) CJAM model results. The top panel shows the mean LOS velocity along the major axis. The bottom panel shows the velocity dispersion along the major axis as solid lines and along the minor axis as dashed lines. The LOS velocity is shown symmetrically along the major and minor axes, to highlight the asymmetric rotational signature. However, the velocity dispersion on the bottom panel is shown for positive values of the axes due to the symmetry of σLOS.

In the text

Fig. B.1. Δχ2 surface for log M200 − c200 plane in iteration 4 (see Table 5 for a summary of the parameters). The black contour lines show the 68th, 90th and 95th percentile confidence levels, with minimal χ2 solution shown with the blue cross. Dotted, dashed, dot-dashed black lines and solid red and black lines show the different literature M200 − c200. The solid lines show the M200 − c200 computed with COMMAH code from Correa et al. (2015) using the cosmological parameters from Planck Collaboration XIII (2016) in red and the cosmological parameters from Komatsu et al. (2009) in black.

In the text

	Fig. C.1. Convergence of the parameters and minimal χ2 for the blue GCs. The symbol shapes and colours correspond to those in Fig. 8 (and are also summarised in the bottom left panel). The panels in the left column show the evolution of the minimal value of the statistics and the right panel the value of the parameter at the minimum χ2.
In the text

	Fig. C.2. Convergence of the parameters and minimal χ2 for the red GCs. The symbols and colours correspond to those in Fig. 8 and columns are explained in Fig. C.1.
In the text

Fig. D.1. Comparison between the first and second moment velocity maps for axisymmetric Jeans modelling in the left and right column, respectively. The first two rows show the best-fit models for blue GCs, with the best fit qDM = 1.38 in the top qDM = 0.7 below. The bottom two rows show the residuals for the blue GCs on the top and red GCs on the bottom. The parameters used to evaluate the model were taken from Table 6 and for qDM = 0.7 we take the other parameters at the minimal χ2 based on the 2D χ2 grids as shown in Fig. 9.

In the text