Intrinsic Shapes of Brightest Cluster Galaxies

Massive elliptical galaxies, when observed projected on the plane of the sky, show smooth elliptical contours with mild boxy deviations (a₄ > 0, Bender & Möllenhoff 1987), small twists (typically ≲10°), and increasing ellipticity $\left(\varepsilon \right)$ profiles toward outer radii (Goullaud et al. 2018). The average projected flattening is $\langle q^{\prime} \rangle \sim 0.8$ (Tremblay & Merritt 1996; Weijmans et al. 2014; Chen et al. 2016; Ene et al. 2018), with $q^{\prime} \equiv 1-\varepsilon$. The presence of isophote twists, together with other factors such as the statistical distribution of ellipticity profiles, provides evidence for the triaxiality of these objects (Illingworth & King 1977; Bertola & Galletta 1978; Vincent & Ryden 2005).

Various works in the last 25 years have studied the average statistical distribution of intrinsic shapes for large galaxy samples at different redshifts z (Tremblay & Merritt 1996; Vincent & Ryden 2005; Chang et al. 2013; Weijmans et al. 2014; Chen et al. 2016; Li et al. 2018a; Ene et al. 2018). Deprojecting these distributions (e.g., Tremblay & Merritt 1996) allows us to recover the average or typical intrinsic shape of the galaxies. These studies find that most massive objects are indeed triaxial, with a mean triaxiality parameter $T=\left(1-{p}^{2}\right)/\left(1-{q}^{2}\right)$ (Franx et al. 1991) in the range $\left[0.4,0.8\right]$ (Vincent & Ryden 2005), where p ≡ b/a, q ≡ c/a, and a ≥ b ≥ c are the lengths of the three principal axes of the density ellipsoid. However, no study has yet attempted to directly measure radially resolved intrinsic shapes of individual galaxies in large samples.

In a recent paper, de Nicola et al. (2020) have presented a triaxial deprojection routine that fits the intrinsic shape of ellipsoidal galaxies and allows them to constrain the viewing angles under which an object is seen by photometric data alone. This can be refined further in combination with the dynamical modeling of appropriate stellar kinematics (S. de Nicola et al. 2022, in preparation), since the number of deprojections that need to be tested is drastically reduced, allowing the study of large samples of galaxies.

An interesting group of massive galaxies consists of the so-called Brightest Cluster Galaxies (BCGs). According to Kluge et al. (2020), a BCG is defined as the closest galaxy to the geometrical and kinematic center of a given galaxy cluster, although not necessarily the most luminous galaxy of the cluster itself. Lying deep in the potential well of the cluster, these giant ellipticals are able to increase their mass through processes such as galaxy mergers (Contini et al. 2018), cannibalism, or tidal stripping (Mo et al. 2008). In a recent paper by Kluge et al. (2020), a sample of 170 BCGs was analyzed in great detail using extremely deep photometric observations, revealing that BCGs follow different scaling relations to ordinary early-type galaxies (ETGs). BCGs are also interesting because their outer parts have probably grown predominantly by collisionless accretion and hence in a manner similar to the proposed growth of (collisionless) dark matter halos.

The first goal of this paper is to constrain the intrinsic shapes and viewing angles of a representative subsample (56 objects) of this BCG catalog with the deprojection method of de Nicola et al. (2020). Then, our second goal is to compare the recovered shapes to those of simulated massive galaxies and their dark matter halos. For this purpose we consider the IllustrisTNG (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018) and Magneticum pathfinder (Hirschmann et al. 2014; Teklu et al. 2015) ³ simulations. These cosmological (magneto)hydrodynamical simulations model the formation and evolution of galaxies in a ΛCDM universe including recipes for star formation and evolution, chemical enrichment of the interstellar medium, gas cooling and heating, and black hole and supernova feedback. These simulations produce galaxy populations with properties in reasonable agreement with observations (Remus et al. 2017; Teklu et al. 2017; Genel et al. 2018; Rodriguez-Gomez et al. 2019; van de Sande et al. 2019; Pulsoni et al. 2020; Remus & Forbes 2021).

The paper is structured as follows. Section 2 describes the galaxy sample used in this work. In Section 3 we explain the deprojection procedure. In Section 4 we present the results on the statistics of triaxial shapes and compare our findings with the TNG and Magneticum simulations. Finally, we draw our conclusions in Section 5. Throughout the paper we assume a flat cosmology with H₀ = 69.6 km s⁻¹ Mpc⁻¹ and Ω_m = 0.286.

The BCGs studied in this work come from a recently published sample (Table 1 of Kluge et al. 2020). Each BCG was observed in the $g^{\prime}$ band with the 2 m Fraunhofer telescope at the Wendelstein Observatory (see Kluge et al. 2020 for technical details). The photometry has exactly the requirements needed for the present work, being extremely deep (down to ${m}_{g^{\prime} }\sim 30$ mag) and reaching very large radii (typically well beyond 100 kpc). From the complete sample, we extract those objects for which supplementary F606W Hubble Space Telescope (HST) photometry (typical resolution ∼015) is available, excluding galaxies that are overall unrelaxed (see below). To combine these high-resolution data with those coming from Wendelstein observations, we first select the radii where we have data from both observation sets that are not affected by seeing (typically from 5'' to 15'' from the center), then we interpolate HST photometry at Wendelstein radii. Finally, we convert the HST data to the g' band by determining the sky level and the scaling factor that minimize the differences between the two photometric sets. We combine the two sets by taking the HST values in the inner 10''–15'' and the Wendelstein values at larger radii. In this way we have photometric data with very high resolution in the center and also extending out to ∼100 kpc for the majority of the objects. ⁴ We complement this list with eight further BCGs which we recently observed in the H and/or Ks bands at the 8.4 m Large Binocular Telescope (LBT) using adaptive optics, with typical resolution of ∼04. We combine the LBT photometry with the Wendelstein one using the approach described above. In Appendix A we show a comparison of the deprojection of two galaxies with and without high-resolution photometry to explore possible photometric effects, showing that they are small. Without high-resolution data, the deprojections cannot probe the central regions of the galaxies, but reliable profiles are derived at larger radii.

This allows us to add 16 more BCGs with only Wendelstein data, for a total of 56 galaxies. The average isophotal flattening $\langle q^{\prime} \rangle$ is ∼ 0.77, although almost every BCG becomes very flat $\left(q^{\prime} \lesssim 0.4\right)$ at large radii. In Appendix B we show the ε and position angle (PA) profiles for every BCG of the sample.

Since BCGs often show signs of interactions with other neighboring galaxies of the cluster or active galactic nuclei (AGN) activity in the central regions (Kluge et al. 2020), and given that our triaxial code works under the assumption of (smooth) "deformed ellipsoids" (see Equation (29) of de Nicola et al. 2020), we omit the innermost/outermost isophotes from the deprojection when we find signs of incomplete relaxation (for example, in the form of bumpy ε or PA profiles). This happens in the very center (e.g., AGNs, ongoing accretion) or in the very outer parts (where dynamical timescales are large). Notes on individual galaxies can be found in Appendix D. Since we are interested in comparing our findings with simulations at large radii, we try, when possible, to extend the deprojection up to 2–4 R_e , with the values for half-light radii R_e taken from Table 4 of Kluge et al. (2020).

2.1. Selection Effects

The full sample of Kluge et al. (2020) is drawn from the Abell–Corwin–Olowin (ACO) catalog (Abell et al. 1989) by adopting redshift- and volume-limiting constraints. It contains BCGs with redshift z ≲ 0.08 (with 15 outliers) and has a slight Malmquist bias (see their Figure 3). However, a comparison of this sample with other large samples, also drawn from ACO, shows that they have about 80%–90% of the objects in common.

In order to check the completeness of our subsample, we define R_N as the ratio of the number of galaxies in our subsample N_sub to the number of galaxies in the full sample N_full. We become progressively incomplete at larger redshifts: at z ≤ 0.04, we have 0.5 ≤ R_N ≤ 1.0, but at larger redshifts R_N ≤ 0.2. Moreover, Figure 1, left panel, shows that the mean redshift and redshift range covered by our subsample are smaller than those of the parent sample. This is expected, since both HST- and LBT-observed galaxies are at lower redshift than the average, and we selected them to perform dynamical modeling.

Zoom In Zoom Out Reset image size

We do not find significant selection effects when considering the size (Figure 1, middle panel) or the total absolute magnitude of the galaxies M_tot (Figure 1, right panel). The results are summarized in Table 1.

Table 1. The Ratio R_N of the Number of Galaxies in Our Sub-Sample N_sub For a Given Interval to the Number of Galaxies in the Full Sample N_full For Three Variables of Interest

Variable	Intervals	RN
	≤0.02	1.0
z	$\left[0.02,0.05\right]$	∼0.5
	≥0.05	∼0.15
Re	≥50	∼0.25
Mtot	$\left[-22,-26\right]$	∼0.25

Note. Our sub-sample is biased towards low-redshift galaxies, while no significant bias in magnitude and effective radii is found.

Download table as:  ASCIITypeset image

In this section we describe the deprojection parameters used for the BCGs. An extensive description of the deprojection routine itself is given by de Nicola et al. (2020). In short, given the observed surface luminosity I_obs = L/pc² on a polar elliptical grid, the code searches for the three-dimensional luminosity density ρ, placed onto an ellipsoidal grid, whose corresponding projected surface luminosity I_fit minimizes rms $=\sqrt{\langle {\left(\mathrm{ln}({I}_{\mathrm{obs}}/{I}_{\mathrm{fit}})\right)}^{2}\rangle }$. The algorithm works under the assumption that a galaxy can be described by what we call a deformed ellipsoid, namely an ellipsoid whose radius is given by

$$\begin{eqnarray}&&{m}^{2-\xi (x)}={x}^{2-\xi (x)}+{\left[\displaystyle \frac{y}{p(x)}\right]}^{2-\xi (x)}+{\left[\displaystyle \frac{z}{q(x)}\right]}^{2-\xi (x)}\end{eqnarray} \tag{ 1 }$$

where the exponent ξ can be used to generate disky (ξ > 0) or boxy (ξ < 0) bias. The three one-dimensional functions p(x), q(x), and ξ(x), along with the density on the x-axis ρ_x (x), specify ρ at each point of the grid. Finally, the code uses a one-dimensional radial smoothing on p(x), q(x), ξ(x), and ρ_x (x) to penalize against unsmooth solutions.

3.1. Choice of Parameters

Table 2 highlights for each BCG the parameters used for the deprojections. In the following we will give a detailed explanation of these parameters.

• 1.  
  Grid sampling. It is important to place both the observed I_obs and deprojected ρ on grids large enough to reproduce photometric information properly, but at the same time a very large grid would slow down the code significantly. Therefore, we start by interpolating the SB onto a finely sampled grid, and then gradually reduce both the number of radial points $\left({n}_{r^{\prime} }\right)$ and angular points $\left({n}_{\theta ^{\prime} }\right)$ as long as the comparison with the observations remains acceptable, that is, an error below 1% for every photometric variable (SB, ε, PA, a₄, a₆). Then, we set the number of radial points of the ρ-grid ${n}_{r}={n}_{r^{\prime} }+20$, while the numbers of angular points n_θ and n_ϕ are typically the same as $\left({n}_{\theta ^{\prime} }\right)$ or slightly larger.
• 2.  
  Grid extent in radius. The innermost radii are the same for both grids, namely the semiminor axis of the innermost isophote. These have to be estimated for every BCG by taking into account the spatial resolution of the observations and by checking whether the galaxy shows central activity. In this last case, central regions are omitted. The outermost radii for the SB grid are estimated by making use of our software for isophotal fitting (Bender & Möllenhoff 1987). We typically stop when the isophotal shape profiles become noisy and have to be set to a constant values (typically at SB ∼ 27–28 mag arcsec⁻²). The largest radii of the ρ-grid are then a few times those of the corresponding SB grid.
• 3.  
  Grid flattenings. The flattening of the SB grid can easily be estimated by considering the isophote PAs. For example, if a galaxy shows isophote structures with the major axes aligned along the vertical axis on the plane of the sky $y^{\prime}$, it makes sense to use an elliptical polar grid flattened in the $x^{\prime}$-direction, even with a 15°–20° twist. As far as ρ is concerned, we first assume a spherical grid, then we re-deproject the galaxies at the best-fit inclination(s) using the recovered 〈p(r)〉 and 〈q(r)〉 as flattenings.
• 4.  
  Smoothing. As shown in S. de Nicola et al. (2022, in preparation) we can recover the true intrinsic density of a triaxial N-body simulation with an rms of ∼10%. Since it is not entirely clear how to estimate the smoothing a priori, we take the four λ-values $\left[{\lambda }_{\rho },{\lambda }_{p},{\lambda }_{q},{\lambda }_{\xi }\right]$ (see Equation (30) of de Nicola et al. 2020) used with the simulation divided by a factor of 2, to take into account that our data are less noisy than the N-body simulation (the smoothing scales as λ⁻²). Since the smoothing value affects the rms one gets at the end (the higher the smoothing, the higher the rms), we verified that for the best galaxies the rms was comparable to the one we got for the simulation. The values we chose are [0.6, 0.03, 0.03, 0.3]. A more rigorous implementation would be the minimization of the Akaike information criterion (AIC, Akaike 1974), as shown by Lipka & Thomas (2021) (see also Thomas & Lipka 2022). We defer this to a forthcoming paper.
• 5.  
  Constraints on p, q. Our code allows the possibility of deprojecting by imposing constraints on p, q. Since the code has shown excellent results in terms of recovering the right profiles, we only impose p, q ≥ 0.2, to prevent solutions that are too flat and may cause problems for the fit.

Table 2. Parameters Used for Deprojections

Galaxy	Conversion (arcsec kpc−1)	${r}_{\min }$	${r}_{\max }$	Iobs grid	ρ-grid	Photometry
2MASX J0753	1.17	2.13	128.1	40 × 15	60 × 16 × 16	L + W
2MASX J0900	1.426	2.35	104.2	50 × 15	70 × 15 × 15	W
2MASX J1358	1.225	1.98	74.04	40 × 10	60 × 11 × 11	W
IC 613	0.653	0.0940	116.1	50 × 10	70 × 11 × 11	H + W
IC 664	0.679	0.0984	76.0	40 × 10	60 × 11 × 11	H + W
IC 1101	1.48	1.06	70.2	40 × 12	60 × 13 × 13	H + W
IC 1565	0.765	0.101	243.2	40 × 10	60 × 11 × 11	H + W
IC 1634	1.336	1.58	156.8	50 × 10	70 × 11 × 11	W
IC 1695	0.987	0.150	149.1	40 × 10	60 × 11 × 11	H + W
IC 1733	0.714	0.528	33.8	30 × 12	50 × 13 × 13	W
IC 2378	0.990	0.762	58.5	40 × 10	60 × 11 × 11	H + W
IC 5338	1.10	4.91	178.3	40 × 10	60 × 11 × 11	H + W
LEDA 1518	1.248	4.16	117.1	40 × 10	60 × 11 × 11	W
LEDA 2098	1.467	2.82	127.3	40 × 10	60 × 11 × 11	W
MCG+01-60	1.16	0.622	38.4	30 × 12	50 × 13 × 13	W
MCG-02-02	1.083	3.44	173.7	40 × 10	60 × 11 × 11	W
MCG+02-04	0.869	0.213	184.3	50 × 12	70 × 13 × 13	H + W
MCG+02-27	0.653	0.114	103.4	40 × 10	60 × 11 × 11	H + W
MCG+02-58	1.52	1.69	100.1	40 × 10	60 × 11 × 11	H + W
MCG+03-04	1.375	2.76	66.34	40 × 10	60 × 11 × 11	W
MCG+03-38	0.886	0.165	87.6	40 × 10	60 × 11 × 11	H + W
MCG+04-28	2.53	1.15	299.7	40 × 12	60 × 13 × 13	H + W
MCG+05-32	1.44	1.19	312.9	40 × 10	60 × 11 × 11	H + W
MCG+05-33	1.23	2.64	67.1	40 × 10	60 × 11 × 11	H + W
MCG+09-13	1.362	3.78	103.1	40 × 10	60 × 11 × 11	W
MCG+09-20	1.296	2.49	133.4	40 × 10	60 × 11 × 11	W
NGC 708	0.332	0.0551	61.2	40 × 15	60 × 15 × 15	H + W
NGC 910	0.354	0.611	70.4	40 × 10	60 × 11 × 11	H + W
NGC 1128	0.486	0.0895	31.7	40 × 10	60 × 11 × 11	H + W
NGC 1129	0.361	0.160	98.0	60 × 12	80 × 12 × 12	H + W
NGC 1275	0.359	0.739	81.9	40 × 10	60 × 11 × 11	W
NGC 2329	0.396	0.534	43.3	50 × 12	70 × 13 × 13	H + W
NGC 2804	0.559	1.58	35.6	40 × 10	60 × 11 × 11	W
NGC 3550	0.703	0.505	99.9	40 × 12	60 × 13 × 13	L + W
NCG 3551	0.640	0.382	33.4	30 × 12	50 × 13 × 13	W
NGC 4104	0.577	0.970	60.2	40 × 12	60 × 13 × 13	L + W
NGC 4874	0.469	0.0905	98.6	50 × 10	70 × 11 × 11	H + W
NGC 6166	0.622	0.718	94.9	40 × 10	60 × 11 × 11	H + W
NGC 6173	0.592	0.0613	109.2	40 × 10	60 × 11 × 11	H + W
NGC 6338	0.552	0.334	132.4	50 × 10	70 × 11 × 11	H + W
NGC 7647	0.818	0.323	97.5	50 × 10	70 × 11 × 11	H + W
NGC 7649	0.835	0.659	94.1	40 × 10	60 × 11 × 11	H + W
NGC 7720	0.611	0.720	153.1	40 × 10	60 × 11 × 11	H + W
NGC 7768	0.545	0.307	102.3	40 × 10	60 × 11 × 11	H + W
SDSS J0837	2.67	1.66	118.0	40 × 20	60 × 20 × 20	L + W
UGC 716	1.19	0.832	175.7	40 × 10	60 × 11 × 11	L + W
UGC 727	1.135	1.89	114.4	40 × 10	60 × 11 × 11	W
UGC 1191	1.21	0.638	124.5	40 × 20	60 × 20 × 20	L + W
UGC 2232	0.958	0.152	95.3	40 × 10	60 × 11 × 11	H + W
UGC 2413	0.690	0.379	113.4	50 × 10	70 × 11 × 11	H + W
UGC 4289	0.587	0.587	50.1	40 × 10	60 × 11 × 11	H + W
UGC 6394	0.847	1.99	97.2	40 × 10	60 × 11 × 11	W
UGC 9799	0.691	4.13	97.7	40 × 10	60 × 11 × 11	H + W
UGC 10143	0.708	0.584	110.6	50 × 10	70 × 11 × 11	H + W
UGC 10726	1.15	0.704	189.8	60 × 12	80 × 12 × 12	L + W
VV 16IC	0.354	0.401	181.8	40 × 10	60 × 11 × 11	L + W

Note. Column 1: galaxy name. Column 2: arcsec kpc⁻¹ conversion factor. Columns 3 and 4: smallest and largest isophotal radii, in kpc. Columns 5 and 6: I_obs and ρ-grid dimensions. Column 7: available photometry (W: Wendelstein, H: HST, L: LBT).

Download table as:  ASCIITypeset image

3.2. Viewing Angles

Using the cutoff described in Section 3.2 we are able to reduce the number of inclinations from the initial value of 1800 by at least a factor of 3. The typical rms values for the best-fit solutions are 0.01–0.03. The results are summarized in Table 3.

4.1. Reliability of the Deprojections

As a first step we verify that our deprojections do reproduce the average photometry of the sample. First, we calculate for every galaxy the mean ε and the twist, defined as ${\rm{\Delta }}\mathrm{PA}=\max (\mathrm{PA})\ -\ \min (\mathrm{PA})$, for both the observed and the recovered photometry. Moreover, we reproject the best-fit densities ρ_g for every galaxy g at three different random viewing angles, computing the same averages as above. This is a good test to statistically verify that the recovered intrinsic shapes are compatible with the observed shape distribution. In Figure 2 we show the histograms for ε (top row) and the twist (bottom row). A Kolmogorov–Smirnov (KS, Kolmogorov 1933; Smirnov 1939) test returns p-values above the canonical 5% threshold ⁶ for both the ε and PA distributions, with this being valid for both the best-fit angles and the reprojections at random viewing angles. This confirms that the recovered photometric variables are statistically representative of the BCG sample.

Zoom In Zoom Out Reset image size

In a second step, we check the distribution of the two best-fit angles 〈(θ, ϕ)〉, which specify the LOS position on the plane of the sky. In the upper panels of Figure 3 we plot the two octants with the best-fit 〈(θ, ϕ)〉 marked on them. We see that there is a lack of solutions near the principal axes, but this is only because we are plotting only the best-fit solution for each galaxy. This clearly disfavours such viewing angles, because isophotal twists cannot occur along the principal axes of an ellipsoidal body. Hence, fits along an assumed LOS that coincides with one of the principal axes will deliver larger values of the rms. In the bottom panel of Figure 3 we show the entire distribution of deprojections over the two octants for the galaxy NGC 708. Solutions on the principal axes are not excluded but lead on average to less good fits. Other examples are provided in the notes in Appendix D.

Zoom In Zoom Out Reset image size

4.2. Distribution of Intrinsic Axis Ratios

We now present the measured shapes of the BCGs. Our deprojection algorithm directly yields the intrinsic axis ratios p(r) and q(r) as a function of the distance from the galaxy center. From these profiles, we compute the triaxiality parameter as a function of the radius T(r). The profiles that we obtain by averaging over all good deprojections are shown in Figure 9.

The left and central panels of Figure 4 show the histograms of the average over all radial bins and over all acceptable deprojections p(r) and q(r) for every galaxy of the sample, i.e., averaged on every "good" deprojection. We get 〈p(r)〉 = 0.84 and 〈q(r)〉 = 0.68, with scatters of ∼0.1 (see Table 3 for the values we get for each BCG). For comparison, Ene et al. (2018) used an ellipticity distribution and found 〈p(r)〉 = 0.88 and 〈q(r)〉 = 0.65 for a sample of slow rotators. The histogram of the mean triaxiality parameter, presented in the right panel of Figure 4, shows that although BCGs follow different scaling relations from ordinary ETGs, they have 0.39 ≤ 〈T〉 ≤ 0.72, in agreement with the findings of Vincent & Ryden (2005) for a sample consisting only of ordinary ETGs. The conclusion here is that the triaxiality is extremely high for every object of the sample, with no object showing a mean triaxiality outside the interval $\left[0.39-0.72\right]$.

Zoom In Zoom Out Reset image size

We do not detect correlations between 〈p(r)〉, 〈q(r)〉, or 〈T(r)〉 and the size. A weak correlation is seen with absolute magnitudes: bright BCGs appear rounder than fainter ones that have approximately the same triaxiality. The trend is more clearly seen when considering the radial profiles (see Figures 5 and 6).

Zoom In Zoom Out Reset image size

4.3. Comparison with the TNG and Magneticum Simulations

In order to compare the recovered shape profiles with the shape profiles of simulated galaxies, we use the IllustrisTNG and Magneticum simulations. We consider the 110.7³ Mpc³ and 68³ Mpc³ cosmological volumes, respectively, as a good compromise between resolution and number of massive galaxies simulated. In TNG100, the mean mass of the stellar particles is 1.4 × 10⁶ M_⊙ while the dark matter (DM) particles have masses of 7.5 × 10⁶ M_⊙. The Plummer equivalent gravitational softening length for both stars and dark matter at redshift z = 0 is r_soft = 0.74 kpc. Magneticum, instead, has stellar particles with masses of 2.6 × 10⁶ M_⊙ and DM particles with masses of 5.1 × 10⁷ M_⊙, while r_soft = 2 kpc for DM and 1 kpc for stars.

We select simulated galaxies with total mass larger than 10¹³ M_⊙ that are the most massive members of their group (so-called 'central'). We divide these galaxies into two mass bins, with the number of objects in each mass bin summarized in Table 4. From these, we derive p(r) and q(r) profiles for the dark matter and the stellar component separately. This is done by diagonalizing the inertia tensor

$$\begin{eqnarray}&&{I}_{{ij}}=\displaystyle \frac{{\sum }_{n}{m}_{n}{x}_{n,i}{x}_{n,j}}{{\sum }_{n}{m}_{n}},\end{eqnarray} \tag{ 2 }$$

where x_n,i are the coordinates of the stellar particles and m_n their mass, calculated in ellipsoidal shells (Zemp et al. 2011). We choose 10 radial bins logarithmically spaced along the intrinsic major axis of the galaxies from 3 to 100 kpc. In each step, the iterative procedure adjusts the flattening of the ellipsoidal shell and the direction of principal axes to the isodensity contours, until it converges within 1% in both p and q. We verified that the variation in the direction of the principal axes is generally within 5° between 3 and 100 kpc and that fixing their position to a mean direction (for example, measured within one effective radius) slightly overestimates the axis ratios by a few percent, up to a median ∼3% in p and 7% in q at 100 kpc. This allows a comparison with the shape profiles derived for our BCGs with our deprojection code, which keeps the direction of principal axes in the 3D deprojected model fixed.

Table 4. The Number of Galaxies for Every Total Mass Bin from the TNG100 and Magneticum Simulations.

${\rm{\Delta }}\mathrm{log}\left({M}_{\mathrm{tot}}/{M}_{\odot }\right)$	TNG	Magneticum
13–13.35	116	31
≥13.35	86	20

Download table as:  ASCIITypeset image

For each radial bin we compare the average profiles with those derived for our BCGs with our deprojection code, doing the same for the triaxiality parameter T(r). We split the BCG sample into a bright one with M_tot < −23.7 and a faint one with M_tot > −23.7, and each of which has 21 galaxies. Assuming an M/L ratio of 6, this corresponds to a stellar mass of 2.2 × 10¹² M_⊙. Similarly, we split the simulated galaxies into two samples considering a total mass cut of 2.2 × 10¹³ M_⊙.

In Figures 5 and 6 we show the comparison between the average profiles for stars and for DM respectively, plotting the BCGs using lines, the simulated TNG100 galaxies using squares, and the simulated Magneticum galaxies using triangles. For the BCGs we also show the rms in each radial bin as an error bar, with a typical value of ∼0.1. This implies a typical error on the mean value of ∼0.02. The rms for the simulated galaxies is of the same order. The left panels show the faint BCG sample together with the less massive simulated galaxies; the right panels show the bright BCG sample together with the more massive simulated galaxies (see Table 4).

Zoom In Zoom Out Reset image size

Bright BCGs appear slightly rounder than faint BCGs by Δp ∼ 0.04 and Δq ∼ 0.08, but with the same triaxiality. This trend is not obvious when looking at simulated galaxies.

The comparison of the profiles of BCGs and simulated galaxies shows that there is a strong disagreement in the inner regions, especially when the simulated stellar component is considered. The disagreement is less pronounced for the simulated DM halos. In particular, p and q of TNG100 galaxies have values down to 0.2–0.3, which would imply the presence of squashed structures in almost all galaxies. The Magneticum galaxies are generally rounder, but still flatter than the observed ones. This shortcoming of simulations in reproducing the correct distribution of the ellipticity of massive (slowly rotating) systems is well documented: different sets of simulations predict a population of flat slow rotators with ellipticities as high as 0.55–0.6 (Li et al. 2018b using Illustris, ⁷ Schulze et al. 2018 using Magneticum, Naab & Burkert 2003 using collisionless N-body simulations). In contrast, the observed ellipticity profiles on the sky (see top panels of Figure 8) demonstrate that most BCGs are round near the center, and it is statistically impossible that all of them are axisymmetric systems viewed close to face-on or pole-on. Turning to the outermost regions, we find a much better agreement between the profiles of BCGs and of simulated galaxies, for both stars and DM. In particular, the p(r) and q(r) profiles of BCGs and DM from TNG (and Magneticum at the high-mass end) follow a similar decreasing trend with a slight offset.

The average profile T(r) of our BCGs is almost flat with r at a value of ≈0.55 with an rms scatter of about 0.08, showing that these objects are overall triaxial. The TNG100 simulation generates objects that tend to be prolate in the center and as triaxial as our BCGs in the outer parts. The Magneticum simulation almost matches the observed average profile at intermediate masses when looking at the dark component, but produces more oblate/prolate profiles in the lower/higher mass bins.

Assuming that the flattening of the stellar component in the simulations compared to the observations is due to implementation in the hydrodynamics scheme and that the dark matter is unaffected by this, the similarity between the observed and simulated DM properties and the fact that the observations show similar triaxiality in the outskirts make it plausible that the light distribution of the outer regions of BCGs is tracing the underlying DM halos and may even allow us to probe the nature of dark matter. In particular, recent N-body simulations (Robertson et al. 2019, Fischer et al. 2022) that study mergers of galaxy clusters show that the shape of dark matter subhalos depends on their physical properties: self-interacting dark matter produces rounder halos than classical ΛCDM.

We have deprojected the photometry of a representative sample of 56 BCGs covering a large radial range with good resolution, from the innermost to the outermost regions and probing into the intracluster light. The deprojection algorithm is able to generate SB profiles that are representative of the observed photometry. Moreover, the results show that the BCGs are consistent with random orientations in space. For the first time, we have measured radial profiles p(r), q(r), and T(r). The recovered shapes point to strongly triaxial galaxies, rounder at the center and flatter at large radii. A comparison with the results of the TNG100 and Magneticum simulations shows that BCGs at large radii are a tracer of the DM halo they are embedded in, possibly probing the nature of dark matter. Extending this analysis to galaxies at higher redshifts can probe the formation history of such objects, although getting SB profiles with high enough signal-to-noise ratio at large values of z certainly represents a challenge for present-day facilities.

The extremely strong triaxiality of these objects stresses the need for triaxial dynamical modeling of the stellar kinematics (e.g., Neureiter et al. 2020) in order to recover unbiased black hole mass and M/L estimates, reconstruct the anisotropy profiles of these galaxies, and evaluate the effects of the different numbers of degrees of freedom. We will address these issues in two forthcoming papers (S. de Nicola et. al 2022, in preparation; B. Neureiter et. al 2022 in preparation).

We thank the anonymous referee for carefully reading the manuscript and providing us with useful comments that helped us to improve the paper. We thank Moritz Fischer for showing us the results of his PhD thesis previous to submission and for interesting discussions.

The Magneticum Pathfinder simulations were partially performed at the Leibniz-Rechenzentrum with CPU time assigned to the Project "pr86re", supported by the DFG Cluster of Excellence "Origin and Structure of the Universe". We are especially grateful for the support by M. Petkova through the Computational Center for Particle and Astrophysics (C2PAP).

The LBT is an international collaboration among institutions in the United States, Italy, and Germany. LBT Corporation partners are: LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, the Astrophysical Institute Potsdam, and Heidelberg University; The University of Arizona on behalf of the Arizona university system; Istituto Nazionale di Astrofisica, Italy; The Ohio State University, and The Research Corporation, on behalf of The University of Notre Dame, University of Minnesota and University of Virginia.

We analyze here the systematic effects stemming from the lack of high-resolution HST or LBT data with Wendelstein observations as well as from the residual degeneracy between the p and q profiles and the viewing angles.

The first point can easily be investigated by deprojecting galaxies with and without high-resolution data, verifying how much the deprojection differs between the two cases. We choose two galaxies, NGC 7647 (for the HST case) and UGC 10726 (for the LBT case), and re-perform the deprojection using Wendelstein data only for the best-fit viewing angles. These two galaxies represent stringent tests given the very low scatter for both p and q (see Table 3) among different solutions at different viewing angles and the code yielding a low best-fit rms value. Moreover, the central regions of these galaxies are relaxed, meaning that we can exploit HST and LBT data up to the innermost radii.

In the two top panels of Figure 7 we show the p(r) and q(r) profiles from the HST(LBT)+Wendelstein case as solid lines, while the Wendelstein-only profiles are shown as dotted lines for p and q. The Wendelstein-only deprojection cannot probe the inner region of the galaxy, but remains within the region delimited by the rms (shown as a colored area in the figures) at larger radii.

Zoom In Zoom Out Reset image size

As a second test, we take the galaxy NGC 7647 and deproject it at the best-fit viewing angles (with HST photometry), stopping the deprojection as soon as the rms reaches $1.5\times {\mathrm{RMS}}_{\min }$. In the bottom panel of Figure 7 we compare this solution to the best-fit one. Also in this case the deviations from the best-fit solution are smaller than the scatter due to the different viewing angles for which an acceptable deprojection is found.

Thus, we conclude the lack of the high-resolution photometry does not change the conclusions reported in this paper. Moreover, considering only the best-fit solution (in terms of the rms) for a given set of viewing angles probes the range of acceptable p and q profiles.

In Figure 8 we show the ε and PA profiles for the BCGs of our sample. Omitted points (see also notes in Appendix D) are not shown.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 9 we show intrinsic axis ratio profiles, along with the corresponding triaxiality profiles, for every BCG of the sample. The profiles are computed by averaging over all acceptable deprojections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

• 1.  
  2MASX J0753: This is one of the galaxies observed at LBT. The point-spread function (PSF) effects are clearly visible way beyond the 04 seeing value estimated during the observations, and therefore all the affected points are not taken into account. The galaxy shows bumpy/noisy ε and PA profiles that cannot be described accurately, despite the fact that the viewing angles are not close to the principal axes.
• 2.  
  2MASX J0900: This well-fitted galaxy has a sudden 25° twist in the outermost regions, which is well reproduced. This may be due to the intracluster light, given that in the innermost the regions there is no significant twist.
• 3.  
  2MASX J1358: This is an example of a galaxy whose ε does not change much as a function of the radius. Our best-fit slightly overestimates it, while the twist is underestimated. We stop the deprojection at 75 kpc since the isophotal parameters cannot be adequately fitted anymore beyond this radius.
• 4.  
  IC 613: The galaxy is very round (ε ≲ 0.1 until beyond 10 kpc). The outermost radii hint at something not in equilibrium, which the code does fit well. The PA at the center is not well reproduced, but since ε is small in the central region, this is not a serious issue.
• 5.  
  IC 664: This is a rare example of a galaxy that is flat in the center and round in the outskirts. This low ε generates an unrealistic twist at large radii, but the galaxy still shows a nice constant PA at lower radii, which is well reproduced by the code.
• 6.  
  IC 1101: The best-fit angles $\left(\theta ,\phi ,\psi \right)$ = (50°, 80°, 160°) almost lie on the (y, z) plane, which does not allow for a precise recovery of the twist. Possibly these angles are not the correct ones, because one of the random reprojections does produce a better fit to the observed twist. The somewhat bumpy ε and PA profiles point to the presence of not fully relaxed structures.
• 7.  
  IC 1565: The galaxy is well fitted with the exception of the outermost ε points, which likely belong to the intracluster light. The huge twist shows some oscillations, hinting at a not completely relaxed galaxy.
• 8.  
  IC 1634: The galaxy has very noisy profiles, which hint at a not yet relaxed galaxy. However, the code reproduces the trends very well.
• 9.  
  IC 1695: Although both the ε and PA profiles are complex, the code reproduces these profiles well, with the exception of the major bump in the PA, which is underestimated by a factor of 2.
• 10.  
  IC 1733: We stop the deprojection at 35 kpc because of possible unrelaxed structures at larger radii. It is the best-fitted galaxy of the sample, with an rms of 0.0073.
• 11.  
  IC 2378: The best-fit angles $\left(\theta ,\phi ,\psi \right)$ = (70°, 10°, 150°) are close to the x-axis. The somewhat noisy but almost constant PA profile allows a good fit with low rms and indicates that the galaxy might be oriented along one of the principal axes.
• 12.  
  IC 5338: We start the deprojection at 4.9 kpc because of possible AGN contamination. The last points do give a suspicious increase in twist, probably linked to the intracluster light; however, this is nicely reproduced.
• 13.  
  LEDA 1518: The galaxy has a 20° twist and a smooth ellipticity profile, although we note that ε > 0.2 in the central regions already. The deprojection reproduces the profiles very well except for the outermost points.
• 14.  
  LEDA 2098: The galaxy has a huge 80° twist, but this is mostly given by the very low ellipticity in the central regions. The code reproduces these profiles well, slightly overestimating ε at the outer radii.
• 15.  
  MCG+01-60: We start the deprojection at 1.16 kpc, since for this galaxy only Wendelstein data are available. This is another galaxy with RMS_best ≤ 0.01.
• 16.  
  MCG-02-02: This well relaxed galaxy has a typical ellipticity profile rising steadily and a small 10° twist. No solutions compatible with viewing angles along the principal axes are found, but we get acceptable deprojections at θ = 10°.
• 17.  
  MCG+02-04: The twist is well reproduced with the exception of a bump around 10 kpc. In the central regions we suspect the galaxy to be not fully relaxed, because of an unrealistic ∼150° twist at small radii. Therefore, we start the deprojection at 1.3 kpc.
• 18.  
  MCG+02-27: The 30° twist present in the central regions is not reproduced very well; however, since ε is very low this does not significantly affect the goodness of the fit. The bumpy ε profile suggests that relaxation is not complete.
• 19.  
  MCG+02-58, MCG+05-33: The best-fit angles $\left(\theta ,\phi ,\psi \right)$ = (80°, 0°, 165°) and $\left(\theta ,\phi ,\psi \right)$ = (60°, 10°, 130°) almost lie along the x-axis and on the (x, z) plane, respectively. This does not allow for a good twist recovery; however, the true PA profiles oscillate around 15° (which would indeed give ψ = 165° if the galaxy were along x) for MCG+05-33 and around 50° for MCG+02-58 (for which ψ = 130° would be the right value). We measure twists oscillating around 20° and 9°, respectively. This indicates that the two galaxies are oriented along one of the principal axes, but possibly not fully relaxed yet.
• 20.  
  MCG+03-04: This interesting galaxy is flat in the central regions, where our code slightly underestimates the ellipticity, and gets rounder in the outskirts. The PA profile is tricky, since the twist is small in the central regions before getting significantly bigger at large radii, where ε is small. The fact that the twist is small in the central regions enables the code to obtain good fits close to the principal axes, as for IC 2378.
• 21.  
  MCG+03-38: Another flat galaxy in the central regions. Both ε and PA jump wildly at ∼10–15 kpc, as if a decoupled structure were there in the center. Nevertheless, the code returns a good fit to the observed photometry.
• 22.  
  MCG+04-28: The high rms (0.122) indicates that the fit is not satisfactory beyond ∼10 kpc, where the galaxy becomes extremely flat. The somewhat bumpy twist is overall small, therefore the systematic offset between model and data is not worrying.
• 23.  
  MCG+05-32: The galaxy is well fitted. The photometry shows an unrealistic twist in the first ∼10 isophotes, which is probably the result of the low ellipticity.
• 24.  
  MCG+09-13, MCG+09-20: These are both galaxies showing typical isophotes of massive ellipticals, although with some bumps. The twist in the central regions for MCG+09-20 is due to ε almost going to 0.
• 25.  
  NGC 708: Although the best-fit viewing angles are not exactly on the principal axes, there are several good solutions compatible with such inclinations, as shown in Figure 3. The photometry has not been deprojected within the first 1.2 kpc because of a dust lane. The scale of the plot in Figure 8 might give the wrong impression of a poorly recovered twist, which is not the case.
• 26.  
  NGC 910: Like NGC 1129, ε goes down and then up again. The 30° twist is well fitted.
• 27.  
  NGC 1128: We do not include the galaxy in the twist histogram, as ε is almost always below 0.1, except for the outermost radii.
• 28.  
  NGC 1129: The best-fit angles $\left(\theta ,\phi ,\psi \right)$ = (60°, 10°, 0°) almost lie on the (x, z) plane. ε goes down to 0 and then increases again. Given that the twist is roughly 90°, this could be a galaxy compatible with viewing angles along the principal axes (as suggested by the best-fit viewing angles) and with intersecting p, q profiles.
• 29.  
  NGC 1275: This peculiar galaxy shows a high ellipticity in both the innermost and the outermost regions, with a dip in between. The 80° twist is very well recovered.
• 30.  
  NGC 2329: We omit the innermost points because of an unrealistic bump in ε.
• 31.  
  NGC 2804: ε decreases toward the outermost regions. In the first 10 kpc, the twist is completely absent. We stop the deprojection at 65 kpc because the isophotal parameters cannot be measured anymore beyond this radius.
• 32.  
  NGC 3550: This galaxy was observed at LBT under poor seeing conditions. Moreover, both the ε and PA profiles hint at a not fully relaxed galaxy.
• 33.  
  NCG 3551: We omit the innermost three isophotes because of resolution problems when deriving the isophotal parameters (only Wendelstein images are available for this galaxy). The deprojection beyond 35 kpc also becomes unfeasible since the galaxy shows signs of nonequilibrium; however, the deprojection yields rms ≤ 0.01.
• 34.  
  NGC 4104: The first 15 arcseconds must be discarded because of poor seeing. Also the outermost points (from 60 kpc) are omitted due to contamination from a neighboring galaxy. It is one of the flattest galaxies of the sample, with ε always between 0.4 and 0.6.
• 35.  
  NGC 4874: Very round galaxy with noisy profiles. We include it in the histogram, although the only region where ε stabilizes above 0.1 is beyond 10 kpc.
• 36.  
  NGC 6166: The best-fit angles $\left(\theta ,\phi ,\psi \right)$ = (40°, 0°, 60°) lie on the (x, z) plane. The low rms might be explained by the fact that with the exception of the innermost radii (where most of the twist occurs but ε is low) the true PA oscillates around the constant PA recovered by the code.
• 37.  
  NGC 6173: We dropped the poorly fitted central region of the galaxy.
• 38.  
  NGC 6338: The deprojection starts at ∼06 because of possible AGN activity.
• 39.  
  NGC 7647: Nothing to remark on here.
• 40.  
  NGC 7649, NGC 7720, NGC 7768: It is not entirely clear how relaxed the galaxies are. The PA profiles are somewhat noisy with very small twists, while the ε profiles increase smoothly with radius (with the exception of the central regions of NGC 7768) with minor dips.
• 41.  
  SDSS J0837: The same considerations about the observations made for 2MASX J0837 also apply for this galaxy. However, this galaxy does not show signs of nonequilibrium.
• 42.  
  UGC 716: The best-fit angles $\left(\theta ,\phi ,\psi \right)$ = (40°, 0°, 60°) almost lie on the (y, z) plane. This is another galaxy that might indeed be close to the principal axes despite the small ∼10° twist near the round center.
• 43.  
  UGC 727: See comments on MCG+09-13 and MCG+09-20.
• 44.  
  UGC 1191: We start the deprojection at ∼0.6–0.7 kpc because of PSF effects. The galaxy has a large twist (∼40°), which is well reproduced.
• 45.  
  UGC 2232: Nothing to remark on here.
• 46.  
  UGC 2413: We note a slight offset in the central regions between the true photometry and the recovered one, probably because of resolution effects given by the spherical ρ-grid.
• 47.  
  UGC 4289: The galaxy shows a somewhat noisy ε profile along with a PA profile with an abrupt ∼30° twist starting from ∼20 kpc. We stop the deprojection at 50 kpc because of possible contamination from neighboring galaxies.
• 48.  
  UGC 6394: The ellipticity decreases as a function of radius and shows some bumps, while the PA profile is much smoother and very well recovered. For this galaxy we also obtain prolate deprojections compatible with the observed photometry.
• 49.  
  UGC 9799: We start the deprojection at 4 kpc to avoid the center, which is affected by probable AGN contamination.
• 50.  
  UGC 10143: The same considerations made for UGC 2413 also apply to this galaxy. We start the deprojection at 0.6 kpc because of a sudden 100° twist in the innermost regions.
• 51.  
  UGC 10726: The same considerations made for UGC 1191 also apply to this galaxy. We omit the outermost isophotes.
• 52.  
  VV 16IC: See comments on UGC 9799.