© The Authors 2024

1. Introduction

Galaxies are complex self-gravitating stellar systems consisting mainly of stars, dark matter, gas, and dust. From early catalogs (Messier 1781; Herschel 1864; Dreyer 1888) through Hubble (Hubble 1926, 1936) and to the present day, the images of the galaxies show a morphological variety that reveals a great diversity of stellar components (Sandage & Tammann 1981).

Analysis of galaxy light distribution allows one to make the most basic morphological distinction. By determining whether a galaxy is dominated by the presence of a rotation-supported disk or a velocity dispersion-supported spheroid, one can categorize a galaxy as being late type or early type. As instrumentation progress grew, this simple classification became quite insufficient because substructures gradually became revealed in galaxies, including the thin and thick disk, the nucleus, the stellar halo, and the bar (Dalcanton & Bernstein 2002; Kormendy & Kennicutt 2004; Trujillo & Bakos 2013). These substructures, typically referred to as stellar components, interact with each other and also follow their own individual temporal evolution. Consequently, the description of the formation and evolution of galaxies is extremely related to the formation and evolution of their components and their coupling.

Cosmological N-body numerical simulations, such as the EAGLE project (Schaye et al. 2015; Crain et al. 2015) or the Next Generation Illustris Simulations (IllustrisTNG) (Pillepich et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Marinacci et al. 2018; Springel et al. 2018), have become one of the fundamental tools for studying galaxies and their evolution over time. The progress in the development and optimization of the implementation of galaxy formation and evolution models has allowed for an increase in particle resolution and the box size of simulations (Vogelsberger et al. 2020). The possibility of having well-resolved galaxies with multiple substructures allows for the study of their physical properties, assembly, and temporal evolution. Thus, to be able to select stellar particles in different components (e.g., thin disk, thick disk, bulge, stellar halo bar), the community has implemented different strategies to identify substructures based on dynamic quantities. The pioneering work in carrying out this task was that of Abadi et al. (2003), who presented a method for selecting which stellar component each stellar particle belonged. Due to its simplicity and efficiency, this dynamical decomposition technique has been widely used with some simplifications or modifications (e.g., Okamoto et al. 2008; Scannapieco et al. 2009; Tissera et al. 2012; Doménech-Moral et al. 2012; Vogelsberger et al. 2014; Marinacci et al. 2014; Obreja et al. 2018; Park et al. 2019; Du et al. 2019; Xu et al. 2019; Gargiulo et al. 2019; Kumar et al. 2021; Zana et al. 2022; Yu et al. 2023; Liang et al. 2024). However, there is no study or agreement on which implementation recovers the components more accurately for specific metrics, such as the completeness and purity of stellar particle classification or the ability to recover physical features such as scale lengths, mass density profiles, velocity dispersion, and rotational velocity profiles. Authors usually choose among some variations of the method available in the literature or modify one of them to carry out their studies. The approach may depend on factors such as the number of stellar components to be identified or the computational cost required to perform the dynamical decomposition of the galaxies of interest. Therefore, it is useful to have the information provided by these different implementations in order to compare the advantages and disadvantages of each, in addition to the differences that arise from using one or the other. With all of this in mind, in this work, we present the dynamical decomposition method named JEHistogram, which is an improvement of the method presented by Abadi et al. (2003) that takes into account the energy distribution of particles as well as the angular momentum distribution of the particles. Its implementation as well as the other five dynamic decomposition methods, namely JThreshold, JHistogram, KMeans, GaussianMixture, and AutoGaussianMixture (see Sects. 2.1 and 2.2 for a brief description of each one), with which the analysis is carried out are present in a software package called GalaxyChop, which enables the dynamical decomposition of galaxies in N-body simulations. Moreover, GalaxyChop can calculate dynamical quantities and easily obtain the results.

This paper is organized as follows. In Sect. 2, we briefly describe the theory behind the different methods of dynamical decomposition of simulated galaxies as well as the methods implemented in the GalaxyChop package. The performance of the dynamical decomposition methods applied to model galaxies is in Sect. 3. In Sect. 4, we dynamically decompose a galaxy belonging to the IllustrisTNG simulations and analyze the results obtained with the different methods. Finally, in Sect. 5, we summarize our results and conclusions.

2. Dynamic decomposition

A simple way to describe the dynamic properties of the particles inside an isolated system in equilibrium is to use the integrals of motion since they do not vary as a function of time. An example of this is the Lindblad diagram (Lindblad 1933), which is constructed based on the energy (Eq. (1)) and angular momentum (Eq. (2)) integrals. Assuming a galaxy as an isolated system in equilibrium with axial symmetry, the orbit of the i-th particle is a fixed point in the Lindblad diagram and can therefore be characterized by its specific total energy,

$$\begin{array}{c}\hfill E_i=\frac{v_i^2}{2}-{\displaystyle \sum _{j\ne i}}\frac{G{m}_j}{r_{\mathit{ij}}},\end{array}$$(1)

where v_i is their total velocity, m_j is the mass of the rest of the particles in the galaxy, and r_ij is the distance between the i-th and j-th particle. The other quantities that characterize the orbit of the particle in the Lindblad Diagram are the z-component of the specific angular momentum of the i-th particle,

$$\begin{array}{c}\hfill J_{z,i}={(x{v}_y-y{v}_x)}_i,\end{array}$$(2)

where x, y correspond to the x and y position and v_x, v_y are the velocities in the x and y directions. We note that it is assumed that the z-axis coincides with the spin axis of the galaxy. Figure 1 shows the Lindblad diagram for the star particles of a simulated galaxy model. The stellar particles in red belong to the spheroid, while the blue ones belong to the disk.

Fig. 1. Lindblad diagram of a model galaxy in equilibrium (see galaxy model 04 in Table 1). The z-component of the angular momentum is shown as a function of the specific binding energy for all stellar particles. In red are represented the particles corresponding to the spheroid, while in blue are the stellar particles of the disk. The black dashed line is the Jcirc. The top inset shows the edge-on view of the galaxy to which the diagram belongs.

In the case of galaxies, an important aspect is that the positions of the particles in the Lindblad diagram are energetically bounded by E < 0 since it is a gravitationally bound system. Also, the inequality −J_circ ≤ J_z ≤ J_circ, where J_circ corresponds to the specific angular momentum of a circular orbit for a given energy value (van den Bosch et al. 1999), provides an upper and lower bound on the region that orbits can populate at in angular momentum values. This is because given an energy value E, the maximum angular momentum value J_circ corresponds to the angular momentum of a star particle in a co-rotating circular orbit, while the minimum angular momentum value −J_circ corresponds to the angular momentum of a star particle in a counter-rotating circular orbit. Therefore, the Lindblad diagram is a tool that allows for detailed analysis of the components of the galaxy through a statistical study of the dynamical properties of the particles that constitute it; since the particles belonging to different components populate different areas of the diagram. Thus, the particles that belong to the disk, which have circular orbits, have values of J_z ∼ J_circ, while the particles that belong to the spheroid have J_z ∼ 0, as can be seen in the Fig. 1. The ratio between the z-component of the angular momentum J_z and the circular angular momentum J_circ of the stellar particles for a given value of energy (ϵ = J_z(E)/J_circ(E)) is called “orbital circularity” or the “circularity parameter” (van den Bosch et al. 1999). This parameter makes it possible to distinguish in the first instance two stellar substructures in a galaxy since stellar particles belonging to the disk component take values of ϵ ∼ 1, while those that form the spheroid are around ϵ ∼ 0.

The analysis of particle clustering in the dynamical space is a powerful tool for classifying stars as belonging to different stellar components. Therefore, most decomposition models use methods that combine some or all of the properties that comprise the dynamic space. With the aim of untangling the stellar components of galaxies, in the following, we present a new decomposition method, JEHistogram, along with five other decomposition models developed by different authors. All six methods have been implemented in the Python package GalaxyChop¹, which is used throughout the rest of the paper for comparison.

2.1. JEHistogram method

Our method improves upon the original implementation by Abadi et al. (2003). We summarize the main assumptions here and refer readers to that paper for further details. The fundamental assumption is that the spheroidal component has zero net angular momentum. All counter-rotating particles – those with a negative z-component of angular momentum (where the positive z-axis is defined by the total angular momentum of all stellar particles) – belong to the spheroidal component. To counterbalance these counter-rotating stars, we add an equal number of co-rotating particles, randomly selected from those with a positive z-component of angular momentum. These co-rotating particles have a circularity parameter ϵ = J_z/J_circ following the same distribution but with a positive rather than negative sign. This approach essentially builds up a spheroidal component with a circularity histogram that is symmetric around ϵ = 0. All remaining galaxy particles are co-rotating and assigned to the disk component, ensuring that no counter-rotating particles are assigned to the disk (see Abadi et al. 2003, Fig. 2, inset lower panel). We note that in the original Abadi et al. (2003) method, there are no restrictions on the binding energy of the particles, so co-rotating particles in the spheroid can have different energies compared to the counter-rotating ones. To correct this asymmetry, we randomly select spheroidal co-rotating particles from all co-rotating stars, ensuring that their binding energy distribution is similar to that of the counter-rotating particles. As a result, the stellar particles belonging to the spheroidal component have an (almost) symmetric circularity distribution, peaked around ϵ = 0, but co-rotating and counter-rotating stars have very similar binding energies. Technically, we built a grid in the ϵ = J_z/J_circ versus E/|E|_max space, where |E|_max is the binding energy of the most bound particle. Thus, instead of working with the original Lindblad diagram (see Fig. 1), we worked with the normalized Lindblad diagram (see as an example Fig. 4) to make it more straightforward to apply to galaxies of a wide range of masses. To decide which co-rotating stellar particles are assigned to the spheroid, we select from each two-dimensional bin with ϵ = J_z/J_circ > 0 the same number of particles that have the bin symmetric to it, with respect to ϵ = J_z/J_circ = 0. Finally, all stellar particles not assigned to the spheroid are classified as part of the disk. In this case we have used 100 bins in ϵ = J_z/J_circ and 20 bins in E/|E|_max.

2.2. Other methods of dynamical decomposition

We briefly describe below the other five different ways of assigning particles to one or another component, implemented in the GalaxyChop package.

• JThreshold: This method is the simplest model used to decompose galaxies. It is an implementation that assigns stellar particles to a spheroid or disk component using only the circularity parameter ϵ. A threshold (ϵ_cut) is defined by the user and the particles with ϵ > ϵ_cut are assigned to the disk and the particles with ϵ ≤ ϵ_cut are assigned to the spheroid component. Several authors use different values of ϵ_cut, such as Vogelsberger et al. (2014), Marinacci et al. (2014), Park et al. (2019). This model is usually implemented to classify galaxies by morphology.

• JHistogram: This method is the implementation of the model described by Abadi et al. (2003). As in the two previous cases, this method decomposes the galaxy into two components: disk and non-rotating spheroid. It is assumed that the distribution of the ϵ for the spheroid component is symmetric around ϵ = 0. For this purpose, the distribution of ϵ < 0 is taken and a distribution of ϵ > 0 reflected around ϵ = 0 is constructed, randomly selecting stellar particles with ϵ > 0. The rest of the particles are assigned to the disk component.

• KMeans: This method is an implementation of the K-means package of scikit-learn (Pedregosa et al. 2011) library. The clustering model of K-means is an unsupervised machine learning technique. In this implementation, the algorithm looks for a specific number of clusters (k) inside a three-dimensional space formed by e/|e|_max, ϵ and ϵ_p, where e is the specific binding energy, |e|_max is the specific binding energy of the most bound particle and ϵ_p = J_proj(E)/J_circ(E) with J_proj = J − J_z is the projected angular momentum. These k clusters are interpreted as different stellar components of the galaxy.

• GaussianMixture: This case is an implementation of the model developed by Obreja et al. (2018). The implementation of this model applies the Gaussian Mixture Model (Pedregosa et al. 2011) algorithm to identify the different clusters in the three-dimensional space formed by e, ϵ, and ϵ_p. In the same way as Kmeans, the number of Gaussians used for the dynamic decomposition is associated with the number of stellar components of the galaxy to be identified. This method associates with each stellar particle a probability of belonging to each Gaussian to be identified.

• AutoGaussianMixture: This is an implementation of the model developed by Du et al. (2019). This model, as the previous one, uses the Gaussian Mixture Model in the parameter space (e, ϵ, ϵ_p). The number of Gaussian is automatically selected using a variation of the Bayesian Information Criterion (BIC) (Schwarz 1978). Each component is associated with a physical structure taking into account the mean values of e, ϵ, and ϵ_p of each Gaussian. In this case, as in the methods described above, the three-dimensional parameter space over which the dynamic decomposition is carried out is (e, ϵ, and ϵ_p). To separate the disk component from the spheroid one, it uses ϵ = 0.5 as the threshold. The disk is subdivided into thick and thin disk using the threshold ϵ = 0.85. The spheroidal component is subdivided into halo and bulge, using, in this case, a threshold on specific normalized energy e = −0.75. We note that the specific energy is normalized to the value of the specific energy of the most bound particle. Similar to the GaussianMixture method, it returns for each stellar particle the probability of belonging to each of the components of the galaxy.

We note that the specific computation time of each method (except AutoGaussianMixture) is practically negligible (∼10⁻² − 10⁻¹) compared to the CPU time invested in the potential energy calculations. Since the implementation of the potential energy calculation is parallelized, the dependence on particle number is much better than theoretical N² and scales approximately as t_CPU ∝ N^0.75 on average.

3. Performance of dynamical decomposition models

To understand the differences obtained in the identification of components by the six variations of the dynamical decomposition, described in Sects. 2.1 and 2.2, we performed an analysis of the assignment of stellar particles to each component. For this purpose, we built a set of nine models of galaxies in equilibrium with a discoidal stellar component and a spheroidal one using the AGAMA package (Vasiliev 2019). Then we assign the stellar particles to the disk and spheroidal components through dynamic decomposition methods and compare some of the physical properties of the resulting components with the original ones.

3.1. Sample of model galaxies in equilibrium

Each of the nine models of galaxies in equilibrium in the sample has a spheroidal stellar component and a discoidal stellar component. They were generated using the AGAMA code (Vasiliev 2019), which allowed the simulation of isolated galaxies in equilibrium from ad hoc initial conditions. The range of stellar masses spanned by the sample is 10 < log₁₀(M_gal/M_⊙) < 12, and the range of spheroidal mass fractions is f_sph = 0.5 − 0.9. Table 1 shows the stellar masses of galaxies, their disk and spheroid components, the mass fraction of the spheroid, and the number of particles in each component. Furthermore, in all cases, the stellar mass density profile of the discoidal component corresponds to an exponential profile. In contrast, the mass density profile of the spheroid corresponds to a power law with an exponential decay (see Table 1 of Vasiliev 2018).

3.2. Stellar component mass fraction

The zeroth-order test compares the stellar component mass fractions that each method reports, concerning the original stellar component mass fractions for each galaxy. In Fig. 2 we present the spheroid mass fraction f_sph (upper panel) and disk mass fraction f_dsk (lower panel) as a function of the M_gal for the nine equilibrium galaxy models and the six implementations of dynamic decomposition methods. As shown in the upper panel of Fig. 2, most of the six methods can reproduce the spheroidal mass with differences relative to the original mass fraction between 0.1% and 17.4%. The JHistogram and JEHistogram methods achieve the best performance for this stellar component since the differences become even smaller than 3.8%. The differences to f_sph for the method JThreshold reach up to 10%, while AutoGaussianMixture achieves up to 17%. The KMeans method underestimates the f_sph regardless of the stellar mass of the galaxy, reaching the largest differences, above 27%, for the three most massive galaxies. On the contrary, it is observed that the GaussianMixture method overestimates the mass fractions in all cases, exceeding differences of 21% for the two less massive galaxies.

Fig. 2. Mass fraction of stellar components. Upper panel: Mass fractions of the spheroidal component fsph as a function of galaxy stellar mass Mgal. Lower panel: Mass fractions of the discoidal component fdsk as a function of galaxy stellar mass Mgal. In both cases, the crosses with black solid line represent the original fsph or fdsk values for the nine equilibrium galaxies built using the AGAMA code (Vasiliev 2019) and the circles with a solid line show the fsph or fdsk values recovered with each of the six dynamical decomposition methods used.

It is interesting to know how the methods recover the disk and spheroid mass fractions, and how these fractions depend on the numerical resolution at which the galaxy is simulated for a given method. We used the fiducial galaxy model (AGAMA-04, see Table 1 and Fig. 1) represented by N ∼ 35 000 stellar particles, and we tested whether the recovered mass fraction f_sph depends on N. We decreased and increased N to N ∼ 3500 and N ∼ 350 000, respectively, and found that f_sph changes from 59% to 65%, compared to 64% for N ∼ 35 000. We note that for the highest resolution tested, N ∼ 350 000, f_sph = 65% has already converged to the original mass fraction with which the model was built, see black curve labeled AGAMA in the upper panel of the Fig. 2. The JHistogram method shows similar convergence behavior, though it is slightly more accurate, while the JThreshold, KMeans, and GaussianMixture methods do not show any convergence, only small fluctuations of the order of 1%.

As the mass fraction of the disk is 1 − f_sph, the trends are reversed. So those methods that overestimate the spheroid underestimate the disk and vice versa (see bottom panel of the Fig. 2).

Although the methods can retrieve the stellar mass fractions of the components in most cases, this does not guarantee that the particles have been assigned to the correct component.

3.3. Purity and completeness

One way to quantify how successful dynamic decomposition methods are is through parameter completeness and purity. We define completeness C as

$$\begin{array}{c}\hfill C=\frac{N_{\mathrm{dyn}}}{N_{\mathrm{AGAMA}}},\end{array}$$(3)

where N_dyn is the number of particles assigned to the component by the method and N_AGAMA is the actual number of particles in the original component. Since the mass of the stellar particles is the same in both components, in the nine equilibrium galaxy models, we can interpret this quantity as the ratio between the mass fraction of the component identified with the dynamical decomposition method (f_sph o f_dsk) and the mass fraction of the original component. In Fig. 2 this corresponds to, in the upper panel, the ratio of colored dots to black crosses, and in the lower panel, the ratio of colored squares to black crosses. However, recovering the correct number of particles per component does not guarantee that these have been correctly assigned to the corresponding stellar component. Additionally, we define the purity P as

$$\begin{array}{c}\hfill P=\frac{N_{\mathrm{dyn},\mathrm{true}}}{N_{\mathrm{dyn}}},\end{array}$$(4)

that is the fraction between the number of particles correctly assigned to the component N_{dyn,  true} and the number of particles assigned to that component N_dyn by the dynamic decomposition methods.

In Fig. 3 we show the results of calculating C and P for the disks and spheroids identified by the six dynamical decomposition methods for the nine equilibrium galaxy models as a function of the stellar mass of the galaxy. The upper-left panel shows C of the disk, while the upper-right panel shows C of the spheroid. The lower-left panel shows the P of the disk, while the lower-right panel shows the P of the spheroid.

Fig. 3. Completeness (C) and purity (P) of the recovered components in the sample of nine equilibrium galaxy models with each of the six methods of dynamical decomposition. The upper panels correspond to the completeness of the disks (left) and the spheroids (right), while the lower panels correspond to the purity of the disks (left) and the spheroids (right). The dashed black lines in the upper panels represent a completeness of one, indicating that the number of particles assigned to the component by the method matches the actual number of particles in the component. The dashed black lines in the lower panels represent a purity of 1, indicating that the number of particles assigned to each component was correctly assigned.

In all galaxies, the spheroids identified have C roughly between 0.66 and 1.25 (Fig. 3, upper-right panel). JHistogram and JEHistogram methods present the smallest variation of C around 1, with differences lower than ∼0.04. GaussianMixture systematically identifies spheroids with 1.01 ≲ C ≲ 1.25 (i.e., more massive), while KMeans systematically identifies spheroids with 0.91 ≲ C ≲ 0.66 (i.e., less massive). Additionally, the values of P are higher than 0.79 for spheroids in the case of the two less massive galaxies, and higher than 0.92 in the rest of the cases, irrespective of the method used (Fig. 3, lower-right panel).

For disks, the differences of C are reversed and exacerbated (Fig. 3, upper-left panel). This is due to the difference in mass between the two stellar components, which can be seen in Table 1. Mis-assigning a particle in the less massive stellar component has a greater impact on C than on the more massive component. GaussianMixture systematically identifies disks with 0.72 ≲ C ≲ 0.89, which is less massive than the originals, although it achieves P values closest to one (Fig. 3, lower-left panel). The C for the KMeans disks is always greater than one and is the method that reaches the highest values of C. This is because the f_dsk for this method, in each galaxy model, is between one to four times the f_dsk original. Consequently, P for these disks goes from 0.86 to 0.24, making this method the one with the worst performance. In addition, as with the spheroids, the JHistogram and JEHistogram methods have the smallest variations in C for the disks with differences of less than 0.06. Of these last two, the JEHistogram method is the one with the highest P values, since the differences of the values of P and one are always smaller than 0.17, for all the simulated galaxy models analyzed in this work (see the red curve in the lower-left panel, Fig. 3).

Therefore, taking into account all the above, we can conclude that the JEHistogram method achieves the best compromise between completeness and purity.

3.4. Stellar components

As the disks and spheroids exhibit different dynamics, the stellar particles of each one are located in different regions in the Lindblad diagram (see as an example Fig. 1). This normalized diagram corresponds to the parameter space where the dynamical decomposition takes place. The z-component of the specific angular momentum is normalized with the specific angular momentum of a circular orbit for a given energy value. Additionally, the specific total energy is normalized by the specific total energy of the most bounded particle. Then, as the dynamical decomposition methods used in this work separate the stellar components in different ways, this results in differences in the distribution of their particles in parameter space of the normalized Lindblad diagram. Therefore, this could be seen in the resulting spatial distribution of the components, as discussed below. We selected the model galaxy AGAMA-04 as an example of the differences that can arise when identifying disk and spheroid in a galaxy using different implementations of dynamical decomposition. Fig. 4 shows the particle distributions of the original disk and spheroid (first column) of this galaxy, and those obtained using the six implementations to decompose them dynamically (six remaining columns). The top row corresponds to the disk component, and the bottom corresponds to the spheroid component. Variations in this diagram are a consequence of the variations in completeness and purity achieved in each case. These differences translate into differences in the properties of the recovered components. Therefore, we analyzed the differences that appear and how they vary for this galaxy model.

Fig. 4. Distribution of stellar particles of the disk (first row) and spheroid (second row) of the AGAMA-04 galaxy in the normalized Lindblad diagram. The first column corresponds to the stellar particle distribution of the original disk and spheroid, while the remaining columns correspond to that obtained by each of the six methods. The color code represents the number of stellar particles in each region.

In the case of the JThreshold method, the assignment of stellar particles depends only on a cut-off circularity parameter ϵ_cut. Therefore, adopting different values of this parameter may result in variations in the mass of the disk and spheroid, in their mass distribution and characteristic velocities. We have considered one of the most commonly used values of ϵ_cut = 0.7 (e.g., Marinacci et al. 2014; Vogelsberger et al. 2014). As can be seen in the second column in Fig. 4, we find that the bulk distribution of the particles in the normalized Lindblad diagram for the stellar components has a certain similarity to the originals. The main difference is the lack of the particles with J_z/J_circ < 0.7 for the disk component and complementary the lack of the particles with J_z/J_circ ≥ 0.7 for the spheroid component.

The JHistogram method presents an improvement over the JThreshold method since it ensures an equal number of counter-rotating and co-rotating stellar particles. This is reflected in its ability to recover the total mass of the spheroid correctly and, consequently, the total mass of the disk (see Fig. 2), regardless of the stellar mass of the model galaxy used. Although the co-rotating particles are selected so that the J_z/J_circ distribution of the spheroid is symmetric about J_z/J_circ = 0, the selection is carried out randomly (see the almost symmetric histogram labeled spheroid in the bottom inset panel of Fig. 2 of Abadi et al. 2003). This impacts the energy distribution of the stellar particles in the resulting components. The Fig. 4 (in the third column) shows how, in the example galaxy, this results in the assignment to the spheroid of a group of stars of high circularity but with energies less bounded to those corresponding to the original spheroid.

Compared to the JHistogram method, the JEHistogram method considers the energy distribution of the counter-rotating particles to improve the assignment of the co-rotating particles to the two stellar components. This is mainly reflected in a better similarity in the original distribution of the particles in the Lindblad diagram especially that corresponding to the spheroidal component as shown in the case of the model galaxy in equilibrium, in the fourth column of the Fig. 4. Thus, taking into account the energy distribution improves the particle assignment. This, together with the fact that this method reproduces the mass fractions better (see Fig. 2), indicates the ability of the method to recover the components correctly and, in consequence, to replicate the surface mass distribution and mass density profiles of the stellar components of galaxies, and their characteristic velocity profiles, as is discussed later.

In the case of the KMeans, its performance depends on the number of particles in each component and their form in the space of parameters. This is because the method performs the separation by minimizing a criterion that measures how coherent the clusters are. Also assumes that the clusters are convex and isotropic, which is not always the case (as an example see Obreja et al. 2016). Therefore, the lack of the ability to distinguish the stellar components of the galaxy reasonably increases with mass, due to the rising difference in the number of particles that form them. This results in less massive (more massive) spheroids (disks) than the original ones (see Fig. 2) and also is reflected in the particle distribution of the resulting stellar components, as a lack (excess) of particles with high circularity in the spheroid (disk) (as we can see, for example, in the fifth column of Fig. 4).

The GaussianMixture method makes a probabilistic assignment of the particles because it assumes that it is a mixture of Gaussians in the parameter space that produces the point distribution. It therefore assigns to each stellar particle the probability of belonging to one or the other stellar component, rather than assigning membership in a taxative way. Despite this, if the distribution of the stellar components in the Lindblad diagram does not have a Gaussian origin, the assignment of the particles to the disk or spheroid may not be optimal. An example of this can be seen in the sixth column of Fig. 4, where the method seems to mis-assign the most bounded particles with values of circularity ∼1. This impacts the mass of the component, resulting in spheroids (disks) that are more (less) massive than the originals (see Fig. 2 as an example), and also in his mass distribution.

Finally, we analyzed the resulting distributions for the components identified by the AutoGaussianMixture method. For this, we took into account that this method identifies by default the halo, bulge, cold disk, and warm disk. Therefore, we decided to unify the halo and bulge in the spheroidal component. Also, we unified the cold disk and warm disk in the discoidal component. This method, developed by Du et al. (2019), presents an improvement in recovering the stellar components concerning the previous method. This is because the particles are distributed in such a way that the use of a single Gaussian to characterize each component seems insufficient. An example of this behavior can be observed in the last column of Fig. 4, which is reflected in the fact that the particle distributions in the normalized Lindblad diagram of the recovered components are more similar to those of the original ones.

3.5. Spatial mass distribution

As an example of the result of the stellar components identified through the different methods of dynamic decomposition, we show in Fig. 5 the spatial mass distribution of the disks and spheroids of one of the nine equilibrium galaxy models included in Table 1. The top two rows show the spatial distribution of the particles identified as belonging to the disk component, while the bottom two rows show those of the spheroid. The bins are colored according to the value of the amount of mass divided by the area of the bin, as indicated in the colorbar on the right-hand side of the plot. From a visual inspection, we can see that the identified components look similar to each other and similar to the original ones. The differences are present where the components overlap spatially, which is in the center of the galaxy because it is where distinguishing which component each stellar particle belongs to becomes more difficult. In the case of the disks, the edge-on views appear indistinguishable from each other, while in the face-on views, the disks appear to have an excess of mass density in the center, except for the case of the GaussianMixture method, which has a hole in the center (see first and second rows of Fig. 5). In the case of the spheroids, also there do not seem to be many differences. The JThreshold, JHistogram, and GaussianMixture methods seem to show misclassified particles in the outer regions (see third and fourth row of Fig. 5).

Fig. 5. Face-on and edge-on view of the surface mass distribution of the disk (first row) and spheroid (second row) corresponding to the AGAMA-04 galaxy. The first column corresponds to the original stellar component distributions, while the rest corresponds to the components obtained by the six dynamical decomposition methods implemented. The color scale represents the value of the surface mass density.

To extend this analysis to a more quantitative one we decided to compare the stellar half-mass radius of the components obtained by all methods for each of the nine model galaxies. In Fig. 6 we present the ratios between the half-mass radius of the component identified by dynamical decomposition and the half-mass radius of the corresponding original component. In the upper panel, the radii were measured in cylindrical coordinates with the disk viewed face-on, while in the lower panel, the radii were measured in spherical coordinates. First of all, regardless of the method used, the spheroids reproduce the radius at half mass with differences of less than 7.4%. The JEHistogram method is the one that presents the smallest variations since, in any case for all simulated galaxy models, they do not exceed 1.2%. If we now analyze the results for the disks, more notable variations appear. However, the JEHistogram method is still the one with the smallest differences in radius at half mass with respect to that of the original components, with values of less than 7.0%.

Fig. 6. Comparison of the half mass radius obtained by the six dynamic decomposition methods and those of the original components of the nine equilibrium galaxy models. Top panel: Comparison of the half mass radius of the disks in cylindrical coordinates. Right panel: Comparison of the half mass radius of the spheroids in spherical coordinates.

3.6. Component mass density profiles

The mass density profiles of the disks are typically exponential (Freeman 1970), while those of the spheroids are de Vaucouleurs profiles (de Vaucouleurs 1959). Thus, they serve to test the ability of different implementations of the dynamical decomposition to reproduce the mass density profiles of the stellar components of the galaxy. Therefore, we performed a comparison between the original mass density profiles and those six recovered by dynamic decomposition. For the disks, we calculated the radial surface mass density profiles in rings Σ = m/π(R_out² − R_in²) in a face-on projection (where R_out is the outer cylindrical radius of the ring, R_in is the inner cylindrical radius of the ring, and m is the mass inside the ring) for the original stellar component Σ_AGAMA and the components identified with the dynamic decomposition Σ_dyn. Also, we estimated the differences between the original disk surface mass density profile. For spheroids, we calculated the radial volumetric mass density profiles in spherical shells ρ = 3m/4π(r_out³ − r_in³), where r_out is the outer spherical radius of the shell, r_in is the inner spherical radius of the shell, and m is the mass inside the shell. These profiles were calculated for the original stellar component ρ_AGAMA and the components identified with the dynamic decomposition ρ_dyn. All profiles were calculated with bins of an equal number of particles. An example of this is shown in Fig. 7, which corresponds to the stellar components of the AGAMA-04 galaxy. The left panel shows the original profile with those of the obtained components (top) and the comparison (bottom) in the case of the disk. The profiles corresponding to the original spheroid and those obtained by the dynamic decomposition are in the top right panel, while the comparison is in the bottom right. The part of the galaxy where identifying which component each stellar particle belongs to becomes difficult corresponds to where the disk and the spheroid spatially overlap, as can be seen in Fig. 7. This refers mainly to the more central part, where it is more difficult to distinguish which component each stellar particle belongs, and can be translated into an excess or deficit in the mass density profiles. For this galaxy, the mass profiles of spheroids identified dynamically show differences lower than 17.2% in comparison with the original profile, except for the last bin where the JHistogram and GaussianMixture methods reach differences of 36.8% and 62.5% respectively. The JEHistogram method achieves the best performance, whereas differences are always lower than 6.5%. This behavior is in agreement with the visual inspection in the two bottom rows in Fig. 5. The discs show differences between the original profile and the profiles of the components recovered by the methods are similar for R > 4.28 kpc. In this range, the differences are less than 19.2% with the JEHistogram method achieving the best performances since differences are less than 3.5%. The exception is the case of the GaussianMixture method, which achieves a difference of 41.4%. In the inner part of the disk, where the disk and spheroid coexist spatially, the impact of mis-assigning a particle is greater. Therefore, in the first bin, the profiles obtained with the dynamic decomposition methods differ by approximately one order of magnitude from the original. We note that the imbalance between the mass of the spheroid and the disk must be taken into account to interpret how the miss-assign impacts the performance of the dynamical decomposition methods. This is because if the mass or number of particles of one of them is much greater than that of the other, the erroneous assignment of one of the stellar particles has a greater impact on the less massive one.

Fig. 7. Surface mass density profiles of the original disk of the AGAMA-04 galaxy and those recovered by the six dynamical decomposition methods (upper-left panel) and comparison between them (lower-left panel) as a function of the cylindrical radius R. Mass density profiles of the original spheroid and those recovered by the different dynamical decomposition methods (upper-right panel) and comparison between them (lower-right panel) as a function of the 3D radius r. In all cases, we observed that the largest differences occur where both components coexist spatially.

3.7. Profiles of characteristic velocity components

Since the disks are supported by rotation, while the spheroids are supported by velocity dispersion, the comparison of the recovered velocity profiles can provide an additional measure of the performance of the dynamic decomposition methods. To this end, we compute the velocity profiles in cylindrical coordinates v_R, v_θ, v_z for the discoidal components and the integrated velocity dispersion in spherical coordinates σ_r(< r), σ_ϕ(< r), σ_θ(< r) for the spheroidal components, of one of the nine models of galaxies in equilibrium as an example. The components of the velocities in cylindrical coordinates are given by

$$\begin{array}{cc}& v_R=(x{v}_x+y{v}_y)/R\hfill \end{array}$$(5)

$$\begin{array}{cc}& v_{\varphi }=(x{v}_y-y{v}_x)/R,\hfill \end{array}$$(6)

where $R=\sqrt{x^2+y^2}$, x, y, z are the position of the particles in Cartesian coordinates and v_x, v_y, v_z are the velocities of the particles in Cartesian coordinates. The integrated profiles of the components of the dispersion velocity in spherical coordinates are given by

$$\begin{array}{c}\hfill {\sigma }_j(<r)=\frac{1}{N(<r)}\sqrt{{\displaystyle \sum _{i=1}^{N(<r)}}{(v_j(r_i)-\overline{v_j}(<r))}^2},\end{array}$$(7)

where j corresponds to one of the spherical coordinates; N(< r) is the number of particles whose radius is less than r of the spheroidal component; v_r = (xv_x + yv_y + zv_z)/r; $v_{\varphi }=(x{v}_y-y{v}_x)/\sqrt{x^2+y^2}$; $v_{\theta }=(v_{r}z-v_{z}r)/\sqrt{x^2+y^2}$; v_j(r_i) is the velocity in direction j of the i-th particle at radius r_i; and $\overline{v_j}(<r)$ is the average velocity in direction j of the particles inside the radius r. In Fig. 8, we present the original and identified velocity profiles for the discoidal component of the galaxy model AGAMA-04. These profiles were calculated in bins with equal particle numbers. In addition, in Fig. 9 we present the integrated velocity profiles as a function of the radius r of the spheroidal components in the same galaxy model, also for the original and recovered components.

Fig. 8. Comparison of velocity profiles as a function of cylindrical radius R of the AGAMA-04 galaxy disk with those identified using the six different methods of dynamical decomposition. Top panel: Comparison of rotational velocity vϕ. Middle panel: Comparison of the radial velocity profile vR. Bottom panel: Comparison of the vertical velocity profile vz. In all cases, the solid black line corresponds to the velocity profile of the original stellar disk. These profiles are calculated with bins with an equal number of particles.

We find that the discoidal components obtained using any of the employed methods exhibit v_R profiles and v_z profiles with values around zero, similar to the case of the original stellar disks, as can be observed in the middle and bottom panels of the Fig. 8. These profiles generally assume values between −7 km s⁻¹ and 8 km s⁻¹. Therefore, the reproduction of these profiles is robust regardless of the method employed for dynamic decomposition. In the case of rotational velocity profiles of the disks identified for this galaxy example (see top panel of Fig. 8), they exhibit differences compared to the original profiles whose relative errors are below 23.2%, especially in the central region, where the disk and spheroid overlap in the galaxy. In the inner bin of the disk profile, which corresponds to the GaussianMixture method, the value of the relative error rises to 54.5%. For all dynamical decomposition methods, these values fall to differences of 5%. This shows the general robustness of the methods in obtaining discoidal components that manage to reproduce their characteristic velocity profiles. Additionally, it confirms that it is more difficult to recover the correct velocity profiles at the center.

Regarding spheroids obtained from dynamical decomposition for this galaxy model, we find that the integrated velocity dispersion profiles in the three directions are reproduced with differences of less than 11.2% concerning the original profiles (see Fig. 9). In all cases, the different extensions of these velocity profiles reflect the position of the stellar particle most distant from the center of the stellar component. We note that the JEHistogram method is the one that most accurately reproduces these profiles because, in the worst-case scenario, the differences are 1.5%.

Fig. 9. Comparison of the integrated velocity dispersion profiles as a function of the spherical radius r of the AGAMA-04 galaxy spheroid with those identified by the six different dynamical decomposition methods. Top panel: Comparison of the radial integrated velocity dispersion profiles σr(< r). Middle panel: Comparison of the azimuthal integrated velocity dispersion profiles σϕ(< r). Bottom panel: Comparison of the polar integrated velocity dispersion profiles σθ(< r). In all cases, the profiles corresponding to the original spheroid are shown with a black solid line.

4. Analysis of stellar component properties in a galaxy from hydrodynamic simulation

In this section, we show an illustrative example of how the properties of the stellar disk and spheroid components can present variations in their properties when different dynamical decomposition methods are used in a galaxy belonging to a hydrodynamic simulation. For this purpose, we select a galaxy belonging to the IllustrisTNG hydrodynamic cosmological simulations (Pillepich et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Marinacci et al. 2018; Springel et al. 2018), which are state-of-the-art simulations. We have decided to choose a galaxy from the TNG100 run. The selected galaxy has a total stellar mass of log₁₀(M_gal/M_⊙) = 10.74, measured within 0.1 × r_vir = 22.7 kpc, with r_vir being the radius of a sphere enclosing an average density 200 times the critical density of the universe, and a rotational velocity v_rot = 219 km s⁻¹ at half mass radius. Both properties were chosen to be similar to those of the Milky Way. Figure 10 displays the distribution of the surface mass density of the galaxy face-on (left) and edge-on (right). It can be observed that the chosen galaxy, used for the subsequent comparison of the identified component properties using different methods of dynamical decomposition, exhibits a discoidal morphology. In addition to the disk component, a spheroidal component is observed in the form of an increased density in the central region of the galaxy. In this case, we rotate the galaxy by aligning the z-axis with the spin of the galaxy calculated with the stellar particles that lie within 0.1 × r_vir.

Fig. 10. Surface mass density of the stellar component for a galaxy from the IllustrisTNG simulations with a stellar mass and rotational velocity similar to those of the Milky Way (ID 486966). Left panel: Face-on view of the galaxy. Right panel: Edge-on view of the galaxy. In both cases, the presence of a disk and spheroid component is evident at first glance. The color scale corresponds to surface stellar mass density.

We dynamically decomposed the selected galaxy using the six methods implemented in the GalaxyChop package. For the case of the JThreshold method, we use ϵ_cut = 0.7. In the KMeans and GaussianMixture methods, we use n = 2 because we want to identify two stellar components: disk and spheroid. Finally, since the AutoGaussianMixture method identifies apriori four components, we call the union of the cold disk and warm disk components a disk, while we call the union of the bulge and halo a spheroid. We compare properties such as stellar mass, surface mass distribution, mass density profiles, and characteristic velocity profiles of the resulting components. We also compare the mass density profiles with those obtained from classical photometric decomposition. Finally, we analyze some astrophysical properties, such as their age and color.

4.1. Mass fraction of the spheroid and mass distribution of the components

In Table 2 we present the stellar masses and the mass fraction of the spheroids and disks identified, using the six dynamical decomposition methods implemented in the package. As we already saw in Sect. 3.2, both disks and spheroids recovered by dynamical decomposition present differences in the amount of stellar mass assigned to them in each case. This results in the f_sph varying between 0.17 − 0.42, with JHistogram and JEHistogram being the methods that recover the less massive spheroidal component.

Figure 11 shows the face-on and edge-on view of the surface mass distribution of the disks (top two rows) and spheroids (bottom two rows) obtained by dynamical decomposition to the IllustrisTNG galaxy. Each column corresponds to each of the six dynamical decomposition methods and the color scale represents the value of the surface mass density. Additionally, Fig. 12 presents the mass density profiles of the stellar components of the IllustrisTNG galaxy recovered by the six dynamical decomposition methods. In the left panel are the surface mass density profiles of the discoidal component as a function of the cylindrical radius R, whereas in the right panel are the volumetric mass density profiles of the spheroidal component as a function of the 3D radius r.

Fig. 11. Face-on and edge-on view of the surface mass distribution of the disk (first row) and spheroid (second row) corresponding to the IllustrisTNG galaxy. Each column corresponds to the stellar component distributions of the components obtained by the six dynamical decomposition methods implemented. The color scale represents the value of the surface mass density.

From a visual inspection of the surface mass density of the disks, we can see that the edge-on distribution seems to be very similar regardless of the dynamic decomposition method used. In the case of the face-on view, some differences appear as gaps in the central region, except for the JHistogram and JEHistogram methods. The presence of this drop in density in the center of the disks is also observed in the surface mass density profile (see left panel of Fig. 12). This behavior, which we already saw in the face-on view of the disk identified by the GaussianMixture method (Fig. 5) and in his corresponding surface mass density profile (Fig. 7), is due to the misclassification that the methods have when assigning a particle from the galaxy central region to the disk or spheroid. If we analyze visually the surface mass densities of the spheroids, we find a roughly similar distribution between the six, as can be seen in the last two rows of Fig. 5. This similarity is reflected in the volumetric density profiles of these stellar components since the differences observed between the profiles do not exceed 0.63 in the logarithm in any case.

Fig. 12. Surface mass density profiles of the discoidal component of the IllustrisTNG galaxy recovered by the six dynamical decomposition methods as a function of the cylindrical radius R (left panel). Mass density profiles of the spheroidal component of the IllustrisTNG galaxy, recovered by the six dynamical decomposition methods as a function of the 3D radius r (right panel). In all cases, we observed that the largest differences occur in the inner part of the disk.

Finally, to complete the analysis quantitatively, we measured the radii at half stellar mass in cylindrical coordinates and the height at half mass of the disks, and the radii at half stellar mass in spherical coordinates of the spheroids. The values obtained are presented in Table 3. We find that the radii for the disks differ by less than 1.41 kpc. These differences are due to different density profile decreases in the central region of the disks, as we can see in the left panel in Fig. 7, since the deeper the profile drop, the larger the radius at half mass of the component. In addition, the height at half mass of the disks differs by less than 0.04 kpc when we compare the values obtained by the six dynamic decomposition methods. For spheroidal components, the half-mass radius differs by at most 0.66 kpc. This further reinforces the idea that the results obtained are robust regardless of the dynamical decomposition method used to identify the stellar components, except for the half-mass radius of the discoidal component.

4.2. Photometric versus dynamical decomposition

Traditionally, to distinguish between the different stellar components of a galaxy, the distribution of their light is usually taken into account. This is done by fitting the radial profile of the total surface brightness through the sum of the radial profiles of the stellar components. The radial surface brightness profile of the spheroidal component follows a de Vaucouleurs profile (de Vaucouleurs 1959), while on the other hand, the disk light distribution is given by an exponential type profile (Freeman 1970). Both profiles are generalized in the Sérsic profile (Sérsic 1963). However, only taking into account the spatial distribution of the stars and ignoring their dynamics makes it much more difficult to elucidate to which component a star in the central region belongs due to the overlapping of the stars. Additional factors that contribute to the complexity of photometric decomposition include the alignment of galaxies to the plane of the sky. This is because a single galaxy, when oriented differently, exhibits a distinct surface brightness distribution. In addition, we have the phenomenon of seeing, causing the blurring of surface brightness distribution as well as extinction, which can notably fade the surface brightness within specific galaxy regions, further amplifying this complexity. It may also happen that the light coming from different regions of the galaxy does not reflect the amount of mass found in that region.

To explore the differences obtained from analyzing the light and mass distribution of the stellar components, we compared the mass and light density profiles for the six dynamical decomposition methods implemented in the GalaxyChop package. For this purpose, we fit a Sérsic profile to each luminosity-weighted and mass-weighted surface density profile, obtaining the values of the n index in each case. All profiles were calculated as a function of the cylindrical radius R from the face-on projection of the components and using bins of an equal number of particles. Moreover, we calculated the cylindrical radii enclosing half the mass and half the light of disks and spheroids. The results are shown in Fig. 13.

Fig. 13. Sérsic index versus cylindrical half mass radius Rhalf for disks and spheroidal half mass radius rhalf for spheroids. The squares represent values obtained for disks and circles for spheroids. The filled symbols represent quantities obtained from the mass-weighted density profiles, and the empty symbols correspond to quantities obtained from the light-weighted density profiles. Each color represents the stellar components obtained with the six dynamic decomposition methods implemented. We note that the error bars for the filled squares are smaller than the size of the symbols.

The first thing we observed in the case of spheroids is that, given one of the dynamical decomposition methods, the values of n and R_half are very similar regardless of whether light or mass is considered. This is because the mass-luminosity relationship of the spheroids is constant throughout, as can be seen in Fig. 14, so the light distribution reflects the mass distribution. The values of R_half are very similar and independent of which dynamical decomposition method was used to identify the component. They range from 1.96 kpc to 2.63 kpc. The n values are around n ∼ 3.

Fig. 14. Ratio between the radial surface density profile of the mass and the luminosity of the galaxy, calculated with bins with an equal number of particles, as a function of cylindrical radius R of the IllustrisTNG galaxy. The black solid line is the mass-luminosity relation considering all the stellar particles. The solid lines correspond to the mass-luminosity relation for spheroids identified with the six methods of dynamic decomposition, while the dashed line corresponds to the discs. The dotted line represents the one-to-one ratio of these quantities.

In contrast, when we perform this analysis for the disks, we find that the values of n and R_half present wide differences whether the light or mass distribution is considered. We find that the range of R_half for the light profiles is between 8.39 kpc and 9.93 kpc, while for the mass profiles, the range of values goes between 5.83 kpc and 7.24 kpc, depending on the dynamical decomposition method used (Fig. 13). This indicates that the light profiles are less concentrated than the mass profiles. The difference is because the mass-to-light ratio of the discoidal components varies with radius (see Fig. 14). In addition, we find a large scatter in the values of the Sérsic indices obtained using the light distribution (n ∼ 1.62 − 3.56), compared to those obtained from the mass distribution (n ∼ 0.95 − 1.66). As already mentioned, these discrepancies are because the light and mass distributions of the galaxy vary along the projected radius (see Fig. 14). In the outer regions where the discoidal component predominates, the light contribution increases compared to the mass, causing the M/L ratio of the galaxy to drop from ∼11 kpc. Another thing we can notice when analyzing the mass-luminosity profile of the whole galaxy is that the largest contribution in light, such as the mass, is due to the disk component. This is due to the similarity between the mass-luminosity profiles of the entire galaxy and the disk, regardless of the method used to identify the stellar components of the galaxy.

These results indicate that the components obtained from the typical photometric decomposition are very different from those obtained from the dynamic decomposition, mainly because the light distribution of the stars does not always reflect the mass distribution of the stars. In the simulations, this behavior seems to be independent of the dynamical decomposition process by which the identification of the stellar components was carried out.

4.3. Characteristic velocity profiles

We compare the velocity profiles in cylindrical (v_R, v_ϕ, v_z) and the integrated velocity dispersion in spherical coordinates (σ_r(< r),σ_θ(< r),σ_ϕ(< r)) for the disk and the spheroid, respectively. Both profiles were obtained similarly to those presented in Sect. 3.7. Figure 15 shows v_ϕ profile at the top panel, v_R profile at the center panel, and v_z profile at the bottom panel. There it can be seen that regardless of the method, v_R and v_z assume values around zero with variations that do not generally exceed 10 km s⁻¹ in absolute value. In the case of the v_ϕ (top panel of Fig. 15) the profiles are similar, regardless of the method used. For R ≲ 5 kpc the differences are less than 20% (which represents ∼30 km s⁻¹), while for R ≳ 5 kpc the differences are reduced to 9% (which represents ∼20 km s⁻¹). This is the same behavior that we had observed in the model galaxy (Fig. 8, in Sect. 3.7), and it is a consequence of the fact that in the central part of the galaxy, it is difficult to assign the stellar particles. The integrated velocity dispersion profiles are shown in Fig. 16. The σ_r(< r) are shown in the upper panel, the σ_θ(< r) are shown in the middle panel, and the σ_ϕ(< r) are shown in the bottom panel. As in the case of the spheroid of the equilibrium galaxy model, all methods identify spheroids with integrated velocity dispersion profiles similar to each other. Regardless of direction, profiles have variations of up to 12%. For r ≲ 2 kpc, a small peak appears, which corresponds to the difficulty of correctly assigning the particles of the innermost region of the components.

Fig. 15. Comparison of velocity profiles identified using the six different methods of dynamical decomposition as a function of cylindrical radius R of the IllustrisTNG galaxy disk. Top panel: Comparison of rotational velocity vϕ. The differences are more evident in the internal part of the disk. Middle panel: Comparison of the radial velocity profile vR. Bottom panel: Comparison of the vertical velocity profile vz. These two profiles adopt values around zero regardless of the method used.

Fig. 16. Comparison of integrated velocity dispersion profiles for spheroids identified using the six different methods of dynamical decomposition as a function of the 3D radius r of the IllustrisTNG galaxy. Top panel: Comparison of the radial integrated velocity dispersion profiles σr(< r). Middle panel: Comparison of the azimuthal integrated velocity dispersion profiles σϕ(< r). Bottom panel: Comparison of the polar integrated velocity dispersion profiles σθ(< r). In all cases, the profiles corresponding to the original spheroid are indicated with a black solid line.

It follows that the differences in the identified components with different implementations of dynamical decomposition methods would have little impact on their characteristic velocity profiles. This impact depends on how the methods assign the particles to the inner part of the galaxy, where both stellar components coexist.

4.4. Age and color profiles of stellar components

Age and color are two fundamental intrinsic properties that provide information about galaxy evolutionary processes. Here we analyze the relationship between the identified stellar components’ dynamics and their astrophysical properties. For this purpose, we calculated the age and g − r color profiles of the disks and spheroids identified by the six dynamic decomposition methods implemented. These profiles were calculated in bins of cylindrical radius R of equal number of particles, in the face-on projection. Figure 17 shows the comparison of the age profiles corresponding to the disk in dashed lines and that of the spheroids in solid lines. The g − r color profiles for the discs in dashed lines and for the spheroids in continuous lines are presented in Fig. 18.

Fig. 17. Comparison of the radial age profiles of the stellar components of the galaxy from IllustrisTNG simulations for the six dynamical decomposition methods as a function of cylindrical radius R. The dashed lines are the profiles corresponding to the disk, while the solid lines are the profiles corresponding to the spheroid.

Fig. 18. Comparison of the radial color g − r profiles of the stellar components of the galaxy from IllustrisTNG simulations for the six dynamical decomposition methods as a function of cylindrical radius R. The dashed lines are the profiles corresponding to the disk, while the solid lines are the profiles corresponding to the spheroid.

We find that all the profiles show a similar trend regardless of the method used. The scatter in the profiles between the methods is < 0.5 Gyr, except for the GaussianMixture method in the spheroid profile. This could be because the method is assigning an excess of particles to the spheroid that the rest of the methods assign to the disk (as can be seen from the drop in the disk mass density profile in the left panel of Fig. 11). The result is a bluer spheroid and a bluer disk compared to the rest. For R ≳ 5 kpc, the components become indistinguishable, indicating that the methods give robust results to each other. Whereas, for R ≲ 5 kpc, the differences are a consequence of the possible error in the assignment of the particles to one stellar component or the other because, in the central part of the galaxy, both components coexist spatially.

In the case of the disks, they all show a negative radial gradient, from the center to the outside of the disk. This gradient changes its slope at R ∼ 11 kpc and becomes more pronounced for R ≳ 11 kpc. This gradient in age is responsible for the change in the M/L ratio, which has a change in slope on approximately the same scale (see Fig. 14). We also analyzed the vertical age profiles of the disks and found that they decrease as we approach the plane of the galaxy, suggesting the presence of a thin disk and a thick disk.

The spheroids have a positive radial gradient (see Fig. 17). These profiles are less pronounced than disk profiles. These profiles also change his slopes, becoming flatter at R ≳ 5 kpc, as a result of the difficulty that the methods have for assigning the stellar particles and that it affects the properties of the less massive stellar component. Another possible explanation for this phenomenon is that the spheroids are the result of the superposition of the bulge and the halo. This results in a slightly younger population in the center and an older population extending to larger radii. As for the vertical age profiles of the spheroids, there is only a variation from younger to older stars as we move away from the plane of the galaxy within the first 5 kpc of radius.

In the case of g − r color profiles (see Fig. 18), as is the case for all the properties discussed so far, we find that regardless of the dynamic decomposition method applied, the profiles are similar to each other for both the disk and the spheroid. The spheroid profiles are redder than the disk profiles, independently of the radius R. Also, the spheroid profiles have an approximately constant value of g − r ∼ 0.74, with a scatter of ∼0.3 in the inner part and ∼0.03 in the outer part except for the GaussianMixture, JHistogram, and JThreshold methods.

As in the case of the age profile, the g − r color profiles recovered with all methods show two different behaviors for R ≲ 11 kpc and R ≳ 11 kpc. On the one hand, for R ≲ 11 kpc, the profile is approximately constant value of g − r = 0.55 and a scatter of 0.05. On the other hand, for R ≳ 11 kpc, it becomes bluer as R increases, reaching a value of g − r = 0.15. This change in the slope of the profile coincides with the radius at which both the galaxy and the disk show a drop in the luminosity-mass ratio (see dashed line in Fig. 14) and a drop in the age profile (see dashed line in Fig. 17), which would suggest that the outer region of the disk is dominated by a young, luminous, bluer stellar population. This could suggest an inside-out disk formation scenario. However, it should be noted that a more exhaustive analysis of the galaxy formation would be necessary to confirm these hypotheses. The differences in the color profiles between the different methods are negligible.

5. Conclusions

We have presented a variation of the dynamic decomposition method of Abadi et al. (2003) that takes into account the energy distribution in addition to the circularity parameter. We implemented this method, together with some other variations in the literature, in a python package (GalaxyChop) that allows one to perform the task of dynamical decomposition of galaxies simply and independently of whether it belongs to a cosmological simulation or an isolated equilibrium galaxy model. Additionally, we analyzed the ability of the six different dynamical decomposition methods implemented to recover the stellar components of nine models of isolated galaxies in equilibrium built with the AGAMA code (Vasiliev 2019). We find that all the methods tested here produce dynamical components with similar global properties, leading to the conclusion that their choice does not dramatically influence the recovered stellar components. However, it is important to note that the choice of the method for performing the dynamical decomposition is decisive for some fundamental properties, such as stellar mass fractions (f_sph and f_dsk) or scale radii, which are under- or overestimated by up to a factor of two (see upper panel Fig. 6). It should also be noted that the largest differences are found in the inner regions of galaxies, where the spatial intersection of the disk and spheroid makes disentangling the two components more challenging. A detailed analysis of purity and completeness showed that the JEHistogram method presented here offers substantial improvements in recovering more detailed galaxy features, such as density profiles and velocity profiles, where the task is more difficult and challenging.

As an example of how our code can be applied not only to galaxies in equilibrium but also to galaxies obtained from a cosmological simulation, we applied the six decomposition methods to a galaxy from the TNG100 simulation with a stellar mass similar to that of the Milky Way. We identified the disk and spheroid using the GalaxyChop package, and we carried out a comparison of some properties. Regarding the mass fractions of the spheroids as well as the stellar mass spatial distribution of the spheroids and disks, we found similar results regardless of the method. This similarity in results was also observed when we compared the half-mass radius and the half-mass height, obtained regardless of the method used. We again note that the innermost part of the galaxy (R ≲ 1.5 kpc) is where it becomes more difficult to carry out the assignment of stellar particles to one component or another. In addition, we compared the profiles obtained by the dynamic decomposition and the classical photometric decomposition, resulting in components that are very different as a consequence of the fact that the mass-luminosity relation varies throughout the galaxy. We find that one should be cautious when comparing dynamic and photometric decomposition. Finally, we find that the dynamically identified disk and spheroid components have different astrophysical properties. The disk has a younger and bluer stellar population than the spheroid. Furthermore, the g − r color profile for the spheroid is approximately constant, while the disk shows a change in both the age profile and the color profile, suggesting that the two components underwent different formation processes. This reinforces the idea that the dynamical processes are intrinsically related to the astrophysical processes of the stellar components.

Acknowledgments

The authors are thankful to the anonymous referee for his comments and suggestions, which helped greatly improve the first version of this paper. VAC and MA thank the financial support from Agencia Nacional de Promoción Científica y Tecnólogica, the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina) and the Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (SeCyT-UNC, Argentina).

References

Appendix A: Application example of the package

As an example of the application of the package, we show how to reach the dynamical decomposition for a simulated galaxy using the JHistogram method. For this example, we used a galaxy selected from TNG100 simulation (Pillepich et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Marinacci et al. 2018; Springel et al. 2018). This dataset is available in GalaxyChop repository for the user to test the functionalities. The models return the star particles and their label indicating the component to which they belong. Depending on the model, it is the number of components in which the galaxy is decomposed.

# Import necessary packages
>>> import numpy as np
>>> import astropy.units as u
# Import GalaxyChop package
>>> import galaxychop as gchop
# Loading the data
>>> gal = gchop.read_hdf5("galaxy.h5")
# Galaxy centring and alignment
>>> galc = gchop.center(gal)
>>> gal = gchop.star_align(galc , r_cut = 30)
# Dynamical decomposition models
# JHistogram Model (Abadi+03)
>>> decomposer = gchop.models.JHistogram()
>>> components = decomposer.decompose(gal)
>>> components.describe()
Particles Deterministic mass
Size Fraction Size
Fraction
Spheroid 12785 0.342918 1.223594e+10
0.32792
Disk 24498 0.657082 2.507786e+10
0.67208
# Index of the component each stellar
# particles belongs to.
>>> labels=
... components.labels[components.ptypes == '
stars']
>>> print(label)
array([0. 1. 0. ... 0. 0. 0.])
# Index = 0: correspond to spheroid particles.
# Index = 1: correspond to disc particles.

A more detailed tutorial of GalaxyChop applications can be found in the repository².

All Tables

All Figures

Fig. 1. Lindblad diagram of a model galaxy in equilibrium (see galaxy model 04 in Table 1). The z-component of the angular momentum is shown as a function of the specific binding energy for all stellar particles. In red are represented the particles corresponding to the spheroid, while in blue are the stellar particles of the disk. The black dashed line is the Jcirc. The top inset shows the edge-on view of the galaxy to which the diagram belongs.

In the text

Fig. 2. Mass fraction of stellar components. Upper panel: Mass fractions of the spheroidal component fsph as a function of galaxy stellar mass Mgal. Lower panel: Mass fractions of the discoidal component fdsk as a function of galaxy stellar mass Mgal. In both cases, the crosses with black solid line represent the original fsph or fdsk values for the nine equilibrium galaxies built using the AGAMA code (Vasiliev 2019) and the circles with a solid line show the fsph or fdsk values recovered with each of the six dynamical decomposition methods used.

In the text

Fig. 3. Completeness (C) and purity (P) of the recovered components in the sample of nine equilibrium galaxy models with each of the six methods of dynamical decomposition. The upper panels correspond to the completeness of the disks (left) and the spheroids (right), while the lower panels correspond to the purity of the disks (left) and the spheroids (right). The dashed black lines in the upper panels represent a completeness of one, indicating that the number of particles assigned to the component by the method matches the actual number of particles in the component. The dashed black lines in the lower panels represent a purity of 1, indicating that the number of particles assigned to each component was correctly assigned.

In the text

Fig. 4. Distribution of stellar particles of the disk (first row) and spheroid (second row) of the AGAMA-04 galaxy in the normalized Lindblad diagram. The first column corresponds to the stellar particle distribution of the original disk and spheroid, while the remaining columns correspond to that obtained by each of the six methods. The color code represents the number of stellar particles in each region.

In the text

Fig. 5. Face-on and edge-on view of the surface mass distribution of the disk (first row) and spheroid (second row) corresponding to the AGAMA-04 galaxy. The first column corresponds to the original stellar component distributions, while the rest corresponds to the components obtained by the six dynamical decomposition methods implemented. The color scale represents the value of the surface mass density.

In the text

	Fig. 6. Comparison of the half mass radius obtained by the six dynamic decomposition methods and those of the original components of the nine equilibrium galaxy models. Top panel: Comparison of the half mass radius of the disks in cylindrical coordinates. Right panel: Comparison of the half mass radius of the spheroids in spherical coordinates.
In the text

Fig. 7. Surface mass density profiles of the original disk of the AGAMA-04 galaxy and those recovered by the six dynamical decomposition methods (upper-left panel) and comparison between them (lower-left panel) as a function of the cylindrical radius R. Mass density profiles of the original spheroid and those recovered by the different dynamical decomposition methods (upper-right panel) and comparison between them (lower-right panel) as a function of the 3D radius r. In all cases, we observed that the largest differences occur where both components coexist spatially.

In the text

Fig. 8. Comparison of velocity profiles as a function of cylindrical radius R of the AGAMA-04 galaxy disk with those identified using the six different methods of dynamical decomposition. Top panel: Comparison of rotational velocity vϕ. Middle panel: Comparison of the radial velocity profile vR. Bottom panel: Comparison of the vertical velocity profile vz. In all cases, the solid black line corresponds to the velocity profile of the original stellar disk. These profiles are calculated with bins with an equal number of particles.

In the text

Fig. 9. Comparison of the integrated velocity dispersion profiles as a function of the spherical radius r of the AGAMA-04 galaxy spheroid with those identified by the six different dynamical decomposition methods. Top panel: Comparison of the radial integrated velocity dispersion profiles σr(< r). Middle panel: Comparison of the azimuthal integrated velocity dispersion profiles σϕ(< r). Bottom panel: Comparison of the polar integrated velocity dispersion profiles σθ(< r). In all cases, the profiles corresponding to the original spheroid are shown with a black solid line.

In the text

Fig. 10. Surface mass density of the stellar component for a galaxy from the IllustrisTNG simulations with a stellar mass and rotational velocity similar to those of the Milky Way (ID 486966). Left panel: Face-on view of the galaxy. Right panel: Edge-on view of the galaxy. In both cases, the presence of a disk and spheroid component is evident at first glance. The color scale corresponds to surface stellar mass density.

In the text

	Fig. 11. Face-on and edge-on view of the surface mass distribution of the disk (first row) and spheroid (second row) corresponding to the IllustrisTNG galaxy. Each column corresponds to the stellar component distributions of the components obtained by the six dynamical decomposition methods implemented. The color scale represents the value of the surface mass density.
In the text

Fig. 12. Surface mass density profiles of the discoidal component of the IllustrisTNG galaxy recovered by the six dynamical decomposition methods as a function of the cylindrical radius R (left panel). Mass density profiles of the spheroidal component of the IllustrisTNG galaxy, recovered by the six dynamical decomposition methods as a function of the 3D radius r (right panel). In all cases, we observed that the largest differences occur in the inner part of the disk.

In the text

Fig. 13. Sérsic index versus cylindrical half mass radius Rhalf for disks and spheroidal half mass radius rhalf for spheroids. The squares represent values obtained for disks and circles for spheroids. The filled symbols represent quantities obtained from the mass-weighted density profiles, and the empty symbols correspond to quantities obtained from the light-weighted density profiles. Each color represents the stellar components obtained with the six dynamic decomposition methods implemented. We note that the error bars for the filled squares are smaller than the size of the symbols.

In the text

Fig. 14. Ratio between the radial surface density profile of the mass and the luminosity of the galaxy, calculated with bins with an equal number of particles, as a function of cylindrical radius R of the IllustrisTNG galaxy. The black solid line is the mass-luminosity relation considering all the stellar particles. The solid lines correspond to the mass-luminosity relation for spheroids identified with the six methods of dynamic decomposition, while the dashed line corresponds to the discs. The dotted line represents the one-to-one ratio of these quantities.

In the text

Fig. 15. Comparison of velocity profiles identified using the six different methods of dynamical decomposition as a function of cylindrical radius R of the IllustrisTNG galaxy disk. Top panel: Comparison of rotational velocity vϕ. The differences are more evident in the internal part of the disk. Middle panel: Comparison of the radial velocity profile vR. Bottom panel: Comparison of the vertical velocity profile vz. These two profiles adopt values around zero regardless of the method used.

In the text

Fig. 16. Comparison of integrated velocity dispersion profiles for spheroids identified using the six different methods of dynamical decomposition as a function of the 3D radius r of the IllustrisTNG galaxy. Top panel: Comparison of the radial integrated velocity dispersion profiles σr(< r). Middle panel: Comparison of the azimuthal integrated velocity dispersion profiles σϕ(< r). Bottom panel: Comparison of the polar integrated velocity dispersion profiles σθ(< r). In all cases, the profiles corresponding to the original spheroid are indicated with a black solid line.

In the text

	Fig. 17. Comparison of the radial age profiles of the stellar components of the galaxy from IllustrisTNG simulations for the six dynamical decomposition methods as a function of cylindrical radius R. The dashed lines are the profiles corresponding to the disk, while the solid lines are the profiles corresponding to the spheroid.
In the text

	Fig. 18. Comparison of the radial color g − r profiles of the stellar components of the galaxy from IllustrisTNG simulations for the six dynamical decomposition methods as a function of cylindrical radius R. The dashed lines are the profiles corresponding to the disk, while the solid lines are the profiles corresponding to the spheroid.
In the text