© ESO 2021

1. Introduction

Galaxies experience several star formation episodes and interaction events during their lifetimes, resulting in complex systems that might have significantly changed their morphological appearance, orbital structure, and stellar population properties. All the different components that a galaxy might contain (such as bulges, disks, bars, shells) are indeed the result of these multiple episodes. Therefore, understanding how a galaxy formed and built up its various components is challenging. One key tool for investigating the formation processes is to disentangle the properties of the stellar populations (kinematics, age, metallicity, chemical abundances, and luminosity) in each structural component (e.g. Corsini et al. 2016; Morelli et al. 2016; Zhu et al. 2018, 2020; Tabor et al. 2019; Oh et al. 2020).

In this respect, galaxies hosting two or more counter-rotating components offer a unique astrophysical laboratory to study how some galaxy components formed and evolved. A stellar and/or gaseous counter-rotating component is characterized by an opposite angular momentum with respect to the main stellar body of the host galaxy. This results in peculiar kinematic signatures that allow us to separate the relative contribution of the different components that are counter-rotating, and to independently study their spectral properties (see Galletta 1996; Bertola & Corsini 1999; Corsini 2014, for reviews).

Large-scale counter-rotating stellar disks are of particular interest among the several kinds of counter-rotation phenomena. The first detection of two large-scale counter-rotating stellar disks was done by Rubin et al. (1992) by studying the shape of the absorption lines in the long-slit spectra of the E7/S0 galaxy NGC 4550. Its counter-rotating disks are co-spatial and have similar sizes, luminosities, and masses (Rix et al. 1992; Johnston et al. 2013), but different thickness (Cappellari et al. 2007). Moreover, one of them is associated with a disk of ionized gas (Rubin et al. 1992; Sarzi et al. 2006; Coccato et al. 2013). In addition to NGC 4550, a number of galaxies hosting two counter-rotating stellar disks with different sizes, luminosities, and masses are also known (e.g., Bertola et al. 1996).

Counter-rotating stellar disks are more easily detected in the outermost regions of galaxies where the seeing blurring is negligible, the contamination of bulge and bar is not relevant, and the velocity separation between the counter-rotating components is maximum giving rise to double-peaked absorption-line spectra. The line-of-sight velocity distribution (LOSVD) turns out to be bimodal, and each of its two parts corresponds to one of the two counter-rotating stellar disks. Several techniques have been developed to perform a spectroscopic decomposition and derive the relative contribution of the counter-rotating components to the measured LOSVD as a function of position onto the sky-plane. In this way it is possible to recover their spatial distribution, kinematics, and stellar population properties (Coccato et al. 2011; Katkov et al. 2011; Johnston et al. 2013; Mitzkus et al. 2017; Méndez-Abreu et al. 2019; Falcón-Barroso & Martig 2021). All these methods suffer severe limitations when the velocity separation and/or light contribution of the two counter-rotating stellar components are small. However, when the spectroscopic decomposition is successful it is possible to constrain the formation timescale of the counter-rotating components and recover the assembly process of the host galaxy (e.g., Pizzella et al. 2014; Morelli et al. 2017; Nedelchev et al. 2019).

Different scenarios involving an external or internal origin were proposed to explain the formation of a counter-rotating stellar disk. It could be the end result of star formation in a counter-rotating gaseous disk, which was accreted through episodic or continuous acquisition of gas clouds from the environment (Thakar & Ryden 1996, 1998) or through the capture of a gas-rich satellite during a minor merger (Thakar et al. 1997; Bassett et al. 2017). Major mergers can form two counter-rotating stellar disks of similar size and mass, but only under some very specific initial conditions (Puerari & Pfenniger 2001; Crocker et al. 2009; Martel & Richard 2020). Another viable external formation process was revealed by cosmological simulations when the effects of the torque of two distinct filamentary structures triggered the star formation onto both prograde and retrograde orbits (Algorry et al. 2014). In addition, Kantharia (2016) suggested that the mutual gravity torques between nearby galaxies can play an important role in changing the disk angular momentum and stressed the importance of environment in the assembly of counter-rotating components. An alternative to the external origin is the internal dynamical instability of the system itself that can lead to the dissolution of a bar or triaxial stellar halo with the formation of a counter-rotating component (Evans & Collett 1994). However, the internal origin does not account for the observed properties of all the known counter-rotating galaxies, which are better explained by mergers and interactions, although the details of this building process are not fully understood yet (Corsini 2014). To this end, we first need to precisely address the demography of counter-rotating stellar disks to determine whether they are a widespread phenomenon in lenticulars and spirals or a class of rare, although fascinating, structures.

Early studies of the frequency of counter-rotating disks of gas and/or stars in lenticulars and spirals were mostly based on long-slit spectroscopic observations of limited samples. The frequency of counter-rotating gaseous disks in lenticulars is ∼30% (Pizzella et al. 2004). This reveals that gas accretion events are not rare in S0 galaxies. On the other hand, only ∼10% of them host a significant fraction (> 5%) of counter-rotating stars (Kuijken et al. 1996). Pizzella et al. (2004) found that < 12% and < 8% of spirals host a gaseous or a stellar counter-rotating disk, respectively. They explained the lower frequency of counter-rotating spirals with respect to counter-rotating lenticulars by the requirement of a larger amount of external gas to first form a counter-rotating gaseous disk that subsequently turns into a counter-rotating stellar disk. Spirals are indeed gas-rich systems and the acquired retrograde gas is swept away if it is less than the pre-exiting prograde one. On the contrary, gas-poor lenticulars can build a counter-rotating gaseous disk even if a small amount of external retrograde gas is acquired.

A step forward in disentangling the different kinematic components of galaxies and in revealing hidden structures like counter-rotating stellar disks has been made possible by integral-field unit (IFU) spectroscopic observations (e.g., Krajnović et al. 2011; Jin et al. 2016; Bryant et al. 2019). The analysis of large samples of galaxies gives important hints about their formation process (e.g., Li et al. 2021), and provides criteria for their classification (e.g., van de Sande et al. 2017). Krajnović et al. (2011) separated the 260 early-type galaxies of the ATLAS^3D survey (Cappellari et al. 2011) into five groups according to their stellar kinematic maps. The so-called double-sigma (2σ) galaxies are characterized by two offset and symmetric peaks of the stellar velocity dispersion along the galaxy major axis. This feature was already identified as the signature of the presence of two counter-rotating disks co-existing in the host galaxy from long-slit data (e.g., Rix et al. 1992; Bertola et al. 1996; Vergani et al. 2007) and used as a diagnostic in the search for counter-rotating galaxies. Krajnović et al. (2011) found that the frequency of 2σ galaxies is 4% when considering this as a lower limit because most of them are lenticular galaxies seen at high inclinations. Moreover, the 2σ signature in some objects could be undetected depending on the instrumental velocity dispersion and the properties of the counter-rotating stellar disks. Though the 2σ feature is the compelling evidence of counter-rotating stellar disks, a dedicated study is still necessary to investigate the observational limits for their detection in stellar kinematic maps.

The detection of a counter-rotating component depends on either instrumental (e.g., spectral and spatial resolution) and physical characteristics of both prograde and counter-rotating components. At present, the most suitable instrument to search for counter-rotating signatures is the Multi Unit Spectroscopic Explorer (MUSE). In particular, this instrument covers a wide wavelength range that includes the calcium triplet (CaT) absorption lines (λλ8498, 8542, 8662 Å) and its instrumental velocity resolution at λ = 8600 Å is σ_MUSE = 38 km s⁻¹. This resolution allows us to better resolve the absorption lines and obtain more precise kinematic measurements. Comparing the MUSE instrument to other IFU spectrographs used for galaxy surveys, such as ATLAS^3D (Cappellari et al. 2011), CALIFA (Sánchez et al. 2012), SAMI (Bryant et al. 2015), and MaNGA (Bundy et al. 2015), only the last is found to cover a wider spectral range, but its instrumental dispersion is larger than 50 km s⁻¹ and its spatial resolution is ∼2 arcsec. The kinematic measurements of MaNGA (and SAMI) are severely impacted by atmospheric seeing (van de Sande et al. 2021).

This work consists of a detailed study of the stellar kinematic maps of mock counter-rotating galaxies observed with MUSE to define new useful diagnostics for identifying the presence of two large-scale counter-rotating stellar disks in real galaxies, and hence addressing their statistics. We want to both investigate hidden and/or unclear kinematic signatures and unveil faint and/or unresolved counter-rotating components. In particular, we focus on those cases where the small velocity separation, similar velocity dispersions, similar stellar populations, or different light contributions make it difficult to detect and analyze the counter-rotating components.

The paper is organized in the following way. In Sect. 2 we explore how the presence of a counter-rotating component shapes the spectrum of a prograde and brighter component, and investigate the goodness of the recovery of the LOSVD. In Sect. 3 we extend our analysis to mock IFU data representing two counter-rotating stellar disks in order to highlight the counter-rotating signatures in the observed kinematic maps. In Sect. 4 we present the analysis of the stellar kinematics. In Sect. 5 we discuss some possible limitations to our approach, and we conclude in Sect. 6.

2. The LOSVD of two counter-rotating components

To search for counter-rotating signatures in the stellar kinematics we need to investigate the possible shapes of the LOSVDs that differ not only because of prograde and retrograde motions but also on account of the relative light contribution of the two decoupled components.

We assume that each component is represented by an E-MILES single stellar population (SSP) model (Vazdekis et al. 2010, 2015) to be convolved with different Gaussian LOSVDs. These models are characterized by spectral resolution FWHM_ssp  =  2.51 Å and spectral sampling (Δλ)_ssp  =  0.9 Å pixel⁻¹, and cover a wide spectral range (1680−50 000 Å) including the CaT triplet absorption lines. The SSP models are based on BaSTI isochrones (Pietrinferni et al. 2004) and bimodal initial mass function (IMF) with slope of 1.30. They span a wide range of ages (0.03 ≤ t ≤ 14 Gyr) and metallicities (−2.27 ≤ [M/H] ≤ 0.40 dex). These models do not explicitly consider any α-enhancement and assume that [M/H] = [Fe/H].

For our analysis we decided to adopt a model with age of t = 6 Gyr and metallicity [M/H] = −0.25 dex. These values are consistent with those measured at 1 r_e by Kuntschner et al. (2010) for the most widely studied counter-rotating galaxy NGC 4550. On the other hand, we do not expect our analysis to change significantly adopting different ages and metallicities. Stellar population properties mainly affect the continuum and equivalent width of the spectral lines and have a secondary effect on their positions. Using a single age and metallicity allows us to limit the number of free parameters and highlights the effects induced by the kinematics. For the two components we used the same SSP model, namely the same age and metallicity, to avoid biasing the fractional light contribution due to intrinsic differences (e.g., young and old stars have different spectral energy distributions) when combining the two templates. However, our aim was to model and measure the relative contribution of the counter-rotating component with respect to the prograde component (or equivalently to the total galaxy) without any attempt to recover the input properties of the stellar populations. Therefore, the choice of the stellar parameters is not relevant because they are equal for both components.

For the prograde component, we fix the relative light contribution to unity and let the other vary considering a multiplicative factor f ranging from 0.1 to 1.0 with a step of 0.1, such that the light contribution of the counter-rotating component with respect to the total galaxy f_cr varies from 0.09 to 0.50. The resulting mock spectrum of the counter-rotating galaxy is given by

$$\begin{array}{c}\hfill S_{\mathrm{m}}(\lambda )=T(\lambda )\ast G(\lambda |V_{\mathrm{pro}},{\sigma }_{\mathrm{pro}})+f\times T(\lambda )\ast G(\lambda |V_{\mathrm{cr}},{\sigma }_{\mathrm{cr}}),\end{array}$$(1)

where λ is the wavelength, T(λ) is the stellar template and G(λ|V, σ) is the Gaussian LOSVD characterized by mean velocity V and velocity dispersion σ. The symbol denotes convolution.

We built a set of mock galaxy spectra in the wavelength range between 8000 Å and 9000 Å that is centered on the CaT triplet. For the convolution, we define the mean velocity and velocity dispersion of the prograde (V_pro, σ_pro) and counter-rotating component (V_cr, σ_cr). We choose V_cr = −V_pro such that the velocity separation ΔV between the two components spans the range from 5 km s⁻¹ to 400 km s⁻¹ with a step of 5 km s⁻¹. This matches the maximum amplitude of stellar rotation curves, which typically does not exceed 200 km s⁻¹. For the velocity dispersion, we fix the value of the prograde component σ_pro = 50 km s⁻¹ and vary σ_cr from 40 km s⁻¹ to 100 km s⁻¹ with a step of 5 km s⁻¹. These values are consistent with those found in NGC 4550 (e.g., Coccato et al. 2013). With this choice the ratio σ_cr/σ_pro covers the interval from 0.8 to 2.0. The range of the parameters of the mock spectra is given in Table 1.

We assume that the velocity and velocity dispersion are the observed quantities along the LOS. Moreover, to mimic observations, we consider an instrumental broadening function and a spectral sampling matching those of MUSE (i.e., FWHM_MUSE = 2.54 Å and (Δλ)_MUSE = 1.25 Å pixel⁻¹). Whereas, we do not add detector, sky, and photon noise to the mock spectra, which is equivalent to having very high signal-to-noise ratio (S/N) data. We can account for noise a posteriori only if we require it.

The absorption lines of the mock spectra appear broadened, asymmetric, or even double-peaked depending on the relative light contribution and the velocity separation of the two counter-rotating components. This demonstrates the large variety of possible LOSVD shapes we have to deal with in the case of counter-rotating stellar systems. Since it is not straightforward to quantify resolved and unresolved stellar counter-rotation, we only define a criterion for the detection of counter-rotating components based on the recovered shape of the LOSVD by using the Gauss-Hermite (GH) parameterization (van der Marel & Franx 1993; Gerhard 1993).

2.1. The Gauss-Hermite parameterization

To understand the counter-rotating signatures we extract the stellar kinematics from our set of mock spectra by using the standard (single) GH parameterization. This analysis allows us to investigate the limitations of the recovery of the LOSVD of counter-rotating systems when adopting this commonly used method and, at the same time, to test a possible detection method for counter-rotating components based on the measurement of non-Gaussian absorption lines.

We use the Penalized Pixel-Fitting code (pPXF; Cappellari & Emsellem 2004; Cappellari 2017), which fits the galaxy spectra by convolving a stellar template with the LOSVD modeled as a truncated GH series. This implementation uses the Hermite polynomials assuming h₀, h₁, h₂ = 1, 0, 0. We adopt as stellar templates a full set of E-MILES SSP models degraded to instrumental resolution. The best-fit solution is a non-negative linear combination of these templates obtained from a χ² minimization in pixel space. In particular, we perform three sets of kinematic measurements: (1) considering v and σ only; (2) including h₃ and h₄; (3) truncating the series at the sixth-order term (V, σ, h₃, h₄, h₅, h₆). We do not include any additive polynomials, but use a multiplicative polynomial of order four to account for the shape of the continuum.

In Fig. 1 we show some examples of best-fit solutions for each kinematic set and different f_cr. We show our mock spectra with their best fit and corresponding residuals. The mock spectra are labeled with the value of the input velocity separation ΔV and are shifted, as are the residuals, to avoid overlapping. As expected, the residuals decrease with increasing order of the truncation of the GH series and with decreasing ΔV and f_cr. The fit performed with only the lower-order moments of LOSVD is obviously poor. In particular, for f_cr = 0.17 (Fig. 1, top panels) or less we do not detect the small contribution of the counter-rotating component, even when including in the fit the higher-order moments. On the contrary, for f_cr = 0.23 (Fig. 1, middle panels), using at least h₃ and h₄, we start detecting the secondary component producing asymmetries of the lines, although the residuals are higher than in the case with f_cr = 0.50 (Fig. 1, bottom panels). This problem could be solved, for example, by removing the constraints of h₁ = h₂ = 0, but this is beyond the scope of this work.

Fig. 1. Best fit (red lines) of some mock spectra (black lines) with their residuals (blue lines) obtained using the GH parameterization for three sets of kinematic measurements: (1) lower-order moments only (v, σ, left panels); (2) four moments (v, σ, h3, h4, middle panels); (3) six moments (v, σ, h3, h4, h5, h6, right panels). Each panel shows mock spectra with decreasing velocity separation from top to bottom as labeled, while the constant contribution of the counter-rotating component increases from top to bottom (fcr = 0.17, 0.23, and 0.50, respectively).

Despite the goodness of the fit for equal contributing components, in realistic cases we do not have very high S/N data. Thus, for comparison we include noise in the spectra of the last example to simulate data with S/N = 50 pixel⁻¹. It is not surprising that the noise reduces the ability of recovering small features or double-peaked lines due to the presence of the secondary component (Fig. 2). However, the broadening of the lines and their strong asymmetries are still recovered as can be seen from trend of the measured parameters. In Fig. 3 we show the measured GH parameters v, σ, h₃, and h₄ as a function of the input velocity separation ΔV, which are obtained from noise-free (top panels) and noisy (bottom panels) spectra with σ_cr/σ_pro = 1.0. These trends show how the different fraction of the counter-rotating component causes: (1) a shift of v from zero to the mean velocity of the primary (i.e., dominant) component V_pro; (2) an increase in σ starting from the input value σ_pro; (3) discernible Gaussian deviations that are given by the non-null values of h₃ and h₄. The most interesting feature is the sudden change that happens in all four parameters when f_cr ≥ 0.33 and ΔV is large (> 150 km s⁻¹). In particular, the parameter v jumps from a low velocity to a higher value, closer to V_pro, while the value of σ drops and then increases again. As expected, h₃ anti-correlates with v, and h₄ correlates with σ. Both exceed values of 0.3 that are usually not observed. Fabricius et al. (2012) studied the kinematics along the major axis of 45 S0-Sc galaxies and found 5 galaxies that exhibit large values of h₃ and h₄. The most extreme case is NGC 3521, which shows |h₃| and |h₄| as large as 0.24 and 0.35, respectively (see also Zeilinger et al. 2001). For this reason it is thought to host a counter-rotating component; however, by applying a spectral decomposition, Coccato et al. (2018) found that the two distinct components are actually a bulge and disk in co-rotation. They explained the counter-rotation as an artifact due to spurious signals that are generated in the Fourier space when the data have a poor spectral resolution. Zeilinger et al. (2001) and Fabricius et al. (2012) used the Fourier correlation quotient method (Bender 1990; Bender et al. 1994) for the LOSVD recovery of NGC 3521.

Fig. 2. As in Fig. 1, but shown are the best fit to mock spectra with fcr = 0.50 and S/N = 50 pixel−1 obtained using four GH moments.

Fig. 3. Comparison of the recovered stellar kinematic parameters in the case of mock spectra without noise (top four panels) and with S/N = 50 pixel−1 (bottom four panels). In both cases the mean velocity v (upper left panel), velocity dispersion σ (upper right panel), and higher-order moments h3 (bottom left panel) and h4 (bottom right panel) of the LOSVD of two counter-rotating components are known as a function of their input velocity separation ΔV. The points correspond to the measured values for mock spectra with fcr = 0.17 (red), 0.23 (orange), 0.33 (light-blue), and 0.50 (black), respectively. The points are connected to highlight the trend. The black solid lines indicate the input value of the mean velocity (v = Vpro), velocity dispersion (σ = σpro), and h3 = h4 = 0 as a function of ΔV. The black dashed lines indicate the default limits to the higher-order moments.

We also found some “off-trend points” (e.g., the orange ones at large ΔV) that are the results of a degenerate problem. Especially in the case where the two components have similar strength (blue line), the code returns different solutions depending on the initial guesses. For example, setting the initial guess of v closer to the velocity of the counter-rotating component the measured v is about zero. On the contrary, if a larger positive initial guess is used then v is closer to the expected value following the trend. In extreme cases, the solutions hit the positive boundary limit that is set to unity indicating that a better fit can be obtained by varying the initial guesses. For example, giving a v that is closer to zero and/or a σ that is larger than σ_pro. The jump from one solution to another is more evident for f_cr = 0.50 in the noisy spectra because here there is no dominant component that can bias the solution. The noise also contributes to smoothing the best-fit solutions, as can be seen in Fig. 2.

For completeness, we analyze these trends when including the variation of σ_cr/σ_pro. For a fixed value of f_cr, we store the mock spectra mimicking a datacube, where the abscissa indicates the input velocity separation ΔV, while the ordinate gives the input σ_cr/σ_pro, but we do not find any specific trend. However, the kinematic signatures of a counter-rotating component can be highlighted by extracting the stellar kinematics when adopting a single GH function: (1) measuring only v and σ makes clear the increase in the velocity dispersion, particularly evident for f_cr = 0.50; (2) including h₃ and h₄ gives a measurement of the strength of the counter-rotating component. The spectral fitting obtained with moments of order greater than four returns improved residuals, but it does not reflect any strong feature that could be helpful for the detection of counter-rotating components. In particular, we use the h₃ parameter to distinguish strong counter-rotation (corresponding to |h₃|> 0.2) from weak counter-rotation (|h₃|≤0.2).

3. Mock IFU observations

After having studied how the presence of two counter-rotating stellar components shapes the LOSVD for given values of ΔV, σ, and f_cr, it is necessary to investigate how this affects the kinematic map. In a real galaxy we expect that ΔV, σ, and f_cr vary with radius and position angle, producing typical patterns in the kinematic maps. To this end, we simulate IFU observations of galaxies hosting two counter-rotating thin stellar disks.

3.1. Modeling counter-rotating galaxies

To describe the physical properties of our galaxy models, we need to define their chemical composition (stellar population parameters), kinematics (LOS velocity and velocity dispersion), and photometry (surface brightness parameters).

The galaxy chemical properties are given by the adopted stellar templates, which are characterized by age and metallicity. The template spectrum is chosen from a stellar library and convolved with a Gaussian LOSVD to match the kinematics along the LOS. During the convolution process also the instrumental characteristics (e.g., the instrumental line spread function) are taken into account. At this early stage we decide to adopt the same stellar template for the two counter-rotating components, whereas the use of different stellar populations in our models and their effects on the kinematic signatures is left to a future analysis.

The observed velocity field is generated projecting into the sky-plane an input rotation curve V_rot(R) typical of disk galaxies, which is characterized by a linear rotation in the central region and flat rotation at larger radii. We define a circular disk of radius R to be projected into the sky-plane considering the inclination i and position angle PA of the galaxy. The observed radius r as function of the spatial coordinates (x, y) is expressed as

$$\begin{array}{cc}\hfill r(x,y)=& \{[-(x-x_0)sin\mathrm{PA}+(y-y_0)cos\mathrm{PA}{]}^2\hfill \\ \hfill & +{{\left[\frac{-(x-x_0)cos\mathrm{PA}-(y-y_0)sin\mathrm{PA}}{cosi}\right]}^2\}}^{\frac{1}{2}},\hfill \end{array}$$(2)

where (x₀, y₀) are the coordinates of the galaxy center. The projected velocity field is given by

$$\begin{array}{c}\hfill v_{\mathrm{los}}(x,y)=V_{\mathrm{rot}}(R)cos\theta sini,\end{array}$$(3)

where

$$\begin{array}{c}\hfill cos\theta =\frac{-(x-x_0)sin\mathrm{PA}+(y-y_0)cos\mathrm{PA}}{r}·\end{array}$$(4)

To keep our models as simple as possible we assume v_cr(x, y) = − v_pro(x, y) to describe the LOS velocity fields of the prograde and counter-rotating disks, while for their LOS velocity dispersion we assume σ_cr(x, y) = σ_pro(x, y) = constant.

The spectra are re-scaled according to the given surface brightness profile I(x, y). To describe the surface brightness of a disk component we adopt the exponential law (Freeman 1970)

$$\begin{array}{c}\hfill I(x,y)=I_{0}exp(-\frac{r}{h}),\end{array}$$(5)

where r is defined as in Eq. (2), while I₀ and h are the central surface brightness and scale length of the disk, respectively. Disk isophotes are ellipses centered on the galaxy center (x₀, y₀) with constant PA and constant ellipticity ϵ = 1 − q′, where q′ is the apparent axial ratio and it is related to galaxy inclination through the relation q′ = cos i. The counter-rotating model will have I_m(x, y) = I_pro(x, y)+I_cr(x, y). The final datacube (i.e., galaxy model) is obtained by adding the contribution of the counter-rotating disk to that of the prograde disk, which is kept fixed. We ignore the contribution of a bulge component, which is very small in most of counter-rotating galaxies, and blurry due to the seeing, since our analysis focuses on the effects in outer region of the galaxy. Furthermore, we do not add sky, detector, or Poisson noise, but they can be modeled and included a posteriori if required for the analysis.

3.2. Validation of the counter-rotating galaxy models

To test the reliability of our approach, we built a kinematic model of IC 719 for which the photometric, kinematic, and stellar population parameters needed to construct the mock data are available in Pizzella et al. (2018).

The decomposition of IC 719 provides very detailed information on the two counter-rotating components that allows us to drop some assumptions previously made to generate our model. For example, we do not use the same stellar population model and equal kinematics for the two counter-rotating components. We use instead two different E-MILES SSP models, one for the main (prograde) disk with mean ages t_pro = 3.5 Gyr and t_cr = 1.5 Gyr for the secondary (counter-rotating) one. We choose mean metallicities [M/H]_pro = +0.26 dex and [M/H]_cr = −0.25 dex, consistent with the available values for the SSP models and measured values at 1.7 kpc. The velocity fields are generated adopting two rotation curves that best reproduce those extracted along the major axis of IC 719, while the velocity dispersion fields are assumed to be constant with values of 60 km s⁻¹ and 40 km s⁻¹ for σ_pro and σ_cr, respectively. The surface brightness distribution of both disks follows an exponential law, although with different scale length h_pro = 1.5 kpc and h_cr = 1.0 kpc. The counter-rotating components equally contribute to the galaxy surface brightness at about 1.0 kpc. We use the relative flux contribution to constrain the central surface brightness of the main component I_0, pro, which is expressed in arbitrary linear units, while the central contribution of the counter-rotating disk is given as a fraction of the main one (i.e., I_0, cr = f₀ × I_0, pro). In the innermost region (< 0.8 kpc) the parameters are poorly constrained, thus we do not consider the central region and we do not even include the bulge contribution. Indeed, IC 719 has a very small bulge.

To directly compare data and model stellar kinematics we spatially re-sample the MUSE datacube of IC 719 with the same binning map. We use the map obtained by the Voronoi tesselletion method developed by Cappellari & Copin (2003). We choose a target signal-to-noise (S/N)_T = 90 and select only spaxels with S/N > 7. This allows us to avoid poor quality spaxels, which could introduce undesired systematic effects in the data at low surface-brightness regimes that are not straightforward to account for. We then extract the stellar kinematics with pPXF limiting the wavelength range to 8400−8800 Å. This short wavelength range permits us to decrease the computational time for the measurements and also allows us to assume that the instrumental line spread function does not vary as function of wavelength. We perform the kinematic measurements following the same scheme used in Sect. 2. The uncertainties are the 1σ formal errors returned by pPXF. The first measurements, obtained fitting only v and σ, were subsequently used as starting guesses for v and σ, when we also fit the higher-order moments h₃ and h₄ of the LOSVD. The best-fit solution of each bin is visually inspected to check the goodness of the fit and when a better solution is expected we re-perform the measurements varying the initial guesses.

We show in Fig. 4 the stellar kinematic maps of IC 719 and the good agreement between the radial profiles of the observed and model kinematics extracted along the galaxy major axis within 1 arcsec aperture. In particular, we nicely recover the reversal of v measured at |r|∼10 arcsec from the center, the remarkable 2σ feature with σ peaking at ∼130 km s⁻¹ for 10 < |r|< 25 arcsec, and the strongly asymmetric LOSVD with the higher-order moments steadily rising to |h₃|∼1.0 and h₄ ∼ 0.8 at |r|∼10 arcsec. We note that we do not fit the MUSE data, but plot the stellar kinematics resulting from the kinematic model built from the known parameters of IC 719.

Fig. 4. Stellar kinematic maps (top panels) and major-axis radial profiles (solid lines with gray error bars, bottom panels) of the mean velocity v (upper left panels), velocity dispersion σ (upper right panels), and higher-order moments h3 (lower left panels) and h4 (lower right panels) of the LOSVD for IC 719 derived from MUSE data. The green dashed lines show the major-axis radial profiles of the kinematic parameters obtained with our model.

3.3. Mock counter-rotating galaxies

We built our mock counter-rotating galaxies assuming the same distance for all models to make possible their direct comparison. We fix the spatial sampling Δx = 50 pc pixel⁻¹, which corresponds to a distance of 51.6 Mpc considering the MUSE spatial sampling of 0.2 arcsec pixel⁻¹. We decided to have the PA = 90° such that the major axis of the galaxy model is aligned with the direction of the x-axis. For these choices, a disk scale length of 1.0 kpc is located at the radius of 20 pixels along the x-axis. Within the field of view (FoV) of 60 arcsec, we cover a physical dimension of 15 kpc.

We model the two counter-rotating components with the same SSP model of age t = 6 Gyr and metallicity [M/H] = −0.25 dex. We consider a rotation curve with a maximum velocity of 100 km s⁻¹. Adopting three different inclinations i = 40°, 50° , and 70° for the two counter-rotating stellar disks their projected maximum velocity separations are ∼130 km s⁻¹, ∼150 km s⁻¹, and ∼190 km s⁻¹, respectively. These Δv address the cases in which the spectral lines are not strongly double-peaked. For both stellar disks we opt for a constant velocity dispersion σ_pro = σ_cr = 55 km s⁻¹ all over the FoV. This value is consistent with that observed at 10 arcsec for the main component of IC 719.

For the photometric parameters of the galaxy models we fix the central surface brightness I_0, pro and scale length h_pro = 1.5 kpc of the prograde disk. We adopt eight different values of the multiplicative factor f₀ such that the central surface brightness of the counter-rotating component can be expressed relative to that of the prograde component, as in Sect. 3.2. The ratio of the surface brightness of the counter-rotating disk to the surface brightness of the model (i.e., f_cr(x, y) = I_cr(x, y)/I_m(x, y)) spatially varies over the FoV, ranging from 0 to 1. We also consider three different values for the scale length of the counter-rotating disk h_cr = 1, 1.5, and 2 kpc. The ratio of the total luminosity of the counter-rotating disk to that of the galaxy model, L_cr/L_T, varies between 0.10 and 0.78 depending on h_cr.

We divide our models into three sets of 24 galaxies each depending on their inclination. Each set of galaxy models is then divided into three subsets according to the scale lengths of the two counter-rotating disks: subset 1 includes the models with more concentrated counter-rotating disks dominating the surface brightness distribution of the galaxy in the inner regions (h_cr < h_pro); subset 2 has the models with constant ratio of the surface brightness of the counter-rotating disk to that of the prograde disk at all radii (h_cr = h_pro); subset 3 counts models with a counter-rotating disk dominating the outer surface brightness distribution of the galaxy (h_cr > h_pro). The photometric parameters characterizing the galaxy models are given in Table 2. The naming convention (miinnn) provides the inclination of the galaxy in degrees (ii) and a sequential number (nnn) to identify the models. For instance, m40001 is the first model with inclination i = 40°.

4. Analysis of stellar kinematic maps

For all our models we perform a standard procedure that can be automatically run for the entire set. In this analysis we apply a spatial re-sampling of the datacubes. and then from each bin we extract the stellar kinematics, as done in Sect. 3.2.

Since our models are noise-free, in principle, our data already have enough S/N to perform the spectral fit, but we want to spatially re-sample the datacube as is usually done for real observations. For this reason, when applying the Voronoi binning, we assume that the noise associated with our data is Poissonian ($N\approx \sqrt{S}$). With this assumption it turns out that for all models (S/N)_T = 150 is a good compromise between our ability to measure the kinematics and the dimension and spatial distribution of the bins. This also avoids effects due to the averaging of spectra over a wide galaxy region. After extracting the stellar kinematics we first look at the qualitative effects due to the presence of the counter-rotating component on the kinematic maps, and then we perform a more quantitative analysis that consists in an automatic detection of the observed feature along the major axis radial profiles of all models. We present the two-moment analysis (Sect. 4.1) and the four-moment analysis (Sect. 4.2).

4.1. Counter-rotating signatures from two-moment analysis

The kinematic maps obtained fitting only v and σ are shown in Fig. 5. By inspecting these maps we confirm that the decrease in the mean velocity and double-peaked velocity dispersion are the main kinematic features due to the presence of two counter-rotating stellar disks.

Fig. 5. Stellar velocity (top panels) and velocity dispersion (bottom panels) maps of models belonging to subset 1 (top row), 2 (middle row), and 3 (bottom row) with an inclination of 50°. The white ellipse shows where fcr = 0.50; this level is not reached in all models. The black dashed lines indicate the velocity separation Δv = 100 km s−1. In the bottom left corner of each velocity map we give the minimum and maximum values of v. These kinematic measurements are obtained fitting only v and σ.

The decrease in the measured parameter v is observed because the mean velocity of the galaxy model is always biased towards the velocity of the dominant component being the one that contributes more than half of the total surface brightness. Thus, the slope and amplitude of the rotation curve decrease with the increasing contribution of the counter-rotating disk. The inversion of the velocity field occurs at the radius corresponding to f_cr = 0.50 (see, e.g., m50008). In subset 2, f_cr is constant over the entire FoV, thus we do not observe a change of sign within the same velocity map, but the direction of rotation changes for models m50013-16 compared to the velocity field of models m50009-11. For model m50012 we find a constant null velocity field because the prograde and counter-rotating disks have equal surface brightness at all radii and give the same contribution to the model (h_cr = h_pro and f_cr = 0.50). In subset 3, for models with f₀ < 1.00, the counter-rotating disk dominates the outer regions, while in the other cases the counter-rotating disk becomes the dominant component at all radii (m50020-24). This last subset can be considered the reverse situation of subset 1 with effects on larger scales.

The size, shape, and slope of the two peaks seen in the velocity dispersion map strongly depend on the velocity separation and contribution of the two counter-rotating stellar disks. It is worth clarifying that we consider as the 2σ feature also what we observe in the models of subset 2, although there is no decrease in σ after the maxima. This is because (1) both stellar disks equally contribute to the model at all radii and (2) at large distances from the center the rotation curves become flat (i.e., the Δv remains constant). Therefore, we observe a biconical shape instead of two symmetric, off-center, roundish regions of higher velocity dispersion. We consider the 2σ feature as the most important kinematic diagnostic for counter-rotating galaxies because it appears even for galaxy models with a regular rotation and no evidence of velocity lowering or reversal. This is the case of m50004, as can be seen from its velocity and velocity dispersion maps in Fig. 5. In subset 3 the 2σ peaks are more elongated than those in subset 1 because the counter-rotating contribution is higher in the outermost than in the innermost regions.

A comparison in Fig. 6 of three representative models, mii004, mii008, and mii012, with different inclinations shows that higher viewing angles emphasize the above features. Highly inclined disks cause an increase in the projected velocity separation between the two counter-rotating components. Thus, for example, the values of σ at the peak locations (r_peak) become larger and the slopes (from center to r_peak) are steeper for higher inclinations. For a quantitative analysis, we extract the radial profiles along the major axis of all models with different inclinations belonging to subset 1, and measure the maximum values of σ and their locations (σ_peak, r_peak). We find that the σ peaks for models of subset 1 are located at radii between ∼0.5 and 1.5 r_e. The value of σ_peak increases with increasing luminosity contribution of the counter-rotating disk to the model, as shown in Fig. 7. From these results we also confirm that the statistics of 2σ galaxies is biased by higher viewing angles.

Fig. 6. Comparison of stellar velocity and velocity dispersion maps (left and right panel of each pair of plots, respectively) of models mii004 (top row), mii008 (middle row), and mii012 (bottom row) with the three adopted inclinations of 40° (right column), 50° (middle column), and 70° (left column), respectively.

Fig. 7. Ratio of the peak velocity dispersion σpeak to the input velocity dispersion σpro as a function of the ratio of the total luminosity of the counter-rotating component Lcr to that of the model LT for subset 1 (filled circles). The yellow, blue, and red lines show the trend for inclination i = 40° ,50°, and 70°, respectively.

4.2. Counter-rotating signatures from four-moment analysis

Analyzing the second kinematic set (v, σ, h₃, h₄) we find that the mean velocity is even closer to the velocity of the brighter component compared to the previous two-moment analysis and the 2σ feature is not always present.

Instead of the 2σ peak we detect a new kinematic feature, which we name 2σ drop. It corresponds to the presence of two minima in the velocity dispersion profile along the major axis with σ_drop ≤ 0.9σ_pro at 0.5 r_e ≲ |r|≲1.5 r_e. This feature is seen only in galaxy models having f_cr > 0.30 and Δv > 100 km s⁻¹. Recovering a 2σ drop strongly depends on the data quality. The 2σ drop disappears in favor of a 2σ peak by adding noise to the kinematic model. In Fig. 8 we show that by adding noise to the model m50008 to have S/N = 5 pixel⁻¹, we do not measure σ lower than the input value, which results in an unphysical solution. This behavior allows us to clarify that we do not measure a real drop in the velocity dispersion at the center of the galaxy. It is an artifact due to adopted GH parameterization of the LOSVD rather than the signature of a peculiar orbital structure of the galaxy. On the contrary, the presence of young stars distributed in a thin and kinematically cold disk in the innermost regions of a galaxy may dominate the light distribution causing a real σ drop in the very center of the velocity dispersion map (Portaluri et al. 2017).

Fig. 8. Stellar kinematic maps of v, σ, h3, and h4 (from left to right) of the LOSVD for the model m50008 derived from data without noise (upper panels) and with S/N = 5 pixel−1 (lower panels). The inner and outer dotted ellipses show where fcr = 0.50 and 0.30, respectively. The solid vertical lines give the position of the 2σ drop detected along the major axis radial profile of the velocity dispersion.

Compared to the previous analysis, the change in v and σ values are related to the use of higher-order moments for the extraction of the LOVSD, which are worth analyzing as well. First of all, we recall that we assume both the prograde and counter-rotating disks to be thin stellar structures with a Gaussian LOSVD, and therefore with null h₃ and h₄ moments. Thus, we do not have a low-velocity tail producing asymmetric deviations traced by the h₃ parameter in any of the velocity distributions along the LOS of the two counter-rotating components, as usually observed in galaxy disks (Binney & Tremaine 1987; Krajnović et al. 2013; van de Sande et al. 2017).

We show in Sect. 2.1 that large h₃ and h₄ are an indication of strongly asymmetric or even double-peaked absorption lines, specifically strong counter-rotation. Hence, we focus on most critical cases of galaxy models with small total light contribution (i.e., the first four models of subset 1), which show weak counter-rotating signatures because of their small values of h₃ and h₄, except for model m70004, which is useful as a comparison.

Our models are characterized by a local anti-correlation of h₃ with v/σ where |v/σ|⪅1, as we see in Fig. 9. This anti-correlation becomes steeper for larger contribution of the counter-rotating component and stronger for larger velocity separations (or equivalently higher inclinations). Models m50004 and m70003 show some points at large h₃. The explanation relies on the possible degenerate solutions of the fitting method. Inspecting the kinematic maps of m50004, we see on the left side of the model that the measured σ is smaller, v is larger, and h₃ and h₄ are larger than on the other side. This means that for the same LOSVD we have two different parameterizations leading to different values of v/σ and h₃. If the solution with the 2σ drop is preferred for model m50004, then the local correlation would resemble that of model m70004 with a discontinuity at v/σ ≈ 1. Looking at the local values of h₃ as a function of x/r_e in Fig. 9, we see that h₃ shows a change in slope at around r_e, being negative for |r|> r_e and positive in the inner regions. This correlation was studied by Fisher (1997) from observations of 18 S0 galaxies. He found a wide variety of behaviors in the local h₃ − v/σ relation. In particular, the anti-correlation is always seen in the inner regions, while in the outer there are three different trends including the correlation between h₃ and v/σ. The high-speed tails can be explained as features that occur primarily within the bulge dominated region of the galaxies. In our case it is exclusively due to the presence of a secondary, non-dominant, counter-rotating component, and probably depends on the interplay of the two different exponential profiles. A systematic study of the higher-order parameters at large radii along the major axis for a large sample of disk galaxies, including those with bars, is necessary to better understand the h₃ behavior and its possible indication of counter-rotating stellar components.

Fig. 9. Local correlation of h3 − v/σ (top panels) and radial trend of h3 as function of x/re (bottom panels) for the first four models of subset 1 (from left to right) with inclination i = 40°, 50°, and 70° (from top to bottom, respectively). The filled circles are the data points that are connected through a solid line to better illustrate the trend.

5. Observational limits

The most critical situation in detecting a counter-rotating component is for the case of weak counter-rotation and when the 2σ feature is not well defined. To better understand whether the h₃ parameter is a good tracer of weak counter-rotation, we analyzed the h₃ radial profile of a simulated galaxy without counter-rotation evolved with an N-body code (Sect. 5.1). For the clear identification of the 2σ peak we decided to start from a very noisy model and to increase the S/N until we recovered this feature (Sect. 5.2).

5.1. Thin disk approximation

A possible limit of our approach is the assumption of two counter-rotating stellar disks being thin. Counter-rotating stellar disks may have a different thickness (Cappellari et al. 2007; Coccato et al. 2013; Pizzella et al. 2018). When we observe a thick stellar disk, the integration along the LOS may give rise to non-zero values for h₃ and h₄ that could affect our kinematic diagnostics preventing us from identifying the presence of a counter-rotating component. This is the reason why we exclude highly inclined disks in our analysis.

However, we quantify this effect by building and analyzing a mock datacube obtained with the MUSE standard setup for a simulated galaxy. In this simulation we model the evolution of a disk galaxy with the N-body+smoothed particle hydrodynamic code GASOLINE (Wadsley et al. 2004). We start with the pure gas corona+dark matter halo system described in Roškar et al. (2008). Over the first 4 Gyr, gas cools and forms stars. At 4 Gyr we stop star formation and let the system evolve to 6 Gyr. At this point we rotate the stars (only) by 180° about the x-axis, and restart star formation, evolving to 11 Gyr. While this is an ad hoc method for producing a counter-rotating galaxy, nonetheless it represents a better approximation to real counter-rotating galaxies than does our thin disk approximation.

We use the same numerical parameters as in Roškar et al. (2008, but see also Roškar et al. 2012). At the end of the simulation, the stellar component consists of 1.8 M particles corresponding to a total mass of 5.7 × 10¹⁰ M_⊙. The counter-rotating component contributes 34% to the total mass.

By considering only the star particles and neglecting the contribution of gas, we generate the luminosity-weighted spaxels of a datacube. We use the SYNTRA code (Portaluri et al. 2017) that translates the photometric, kinematic, and chemical properties of the simulated galaxy assigning a SSP model to each stellar particle. Here, the adopted SSP spectra are from the GALAXEV library computed using the isochrone synthesis code of Bruzual & Charlot (2003). For each simulation time step and given galaxy orientation (i.e., i and PA), SYNTRA allows us to build a datacube covering a FoV of 5.3 × 5.3 kpc².

We consider mock observations of the simulated galaxy with a PA = 0° and two viewing angles of 60° and 80°. For each of the two inclinations, we build two datacubes: (1) sim0-11i60/80 including all stars, with ages of 0−11 Gyr and (2) sim6-11i60/80 considering only stars with ages of 6−11 Gyr. The former includes both prograde and retrograde stars, while the latter only prograde stars. The stellar counter-rotating component is younger (t < 6 Gyr).

For the four datacubes, we measure the stellar kinematics with pPXF in the same way as we do for our kinematic models. In Fig. 10 we show the v, σ, h₃, and h₄ maps of the datacubes sim0-11i60/80 and sim6-11i60/80. Considering both prograde and retrograde star particles, the mock observations of the simulated galaxy further confirm the kinematic signatures of two counter-rotating stellar disks, in particular the 2σ signature, which is stronger at higher inclinations. In the case with only a prograde disk, the 2σ signature is no longer present, while the values of |h₃| only reach a maximum of 0.05. The local anti-correlation of h₃ with v is limited to the bulge region, which is identified by the higher velocity dispersion at the galaxy center. These results support the stellar counter-rotation diagnostics, which we investigated by our simple kinematic models relying on the thin disk approximation.

Fig. 10. Stellar kinematic maps of v, σ, h3, and h4 (from left to right) of the LOSVD for the simulated galaxy with (sim0-11i60/80, top panels) and without (sim6-11i60/80, bottom panels) counter-rotating stars.

5.2. The effect of the S/N

To investigate whether h₃ is a good diagnostic also in the case of weak counter-rotation when the 2σ feature is not well defined and S/N is poor, we examine the galaxy model m50001 in more detail because it shows a very faint 2σ signature. The counter-rotating disk of m50001 is characterized by a total luminosity ratio L_cr/L_T = 0.11 and a maximum surface brightness contribution f_cr = 0.20 in the center and decreasing outwards. By measuring only the first two moments of the LOSVD of the noise-free model, we find σ_peak/σ_pro ∼ 1.2. We measure a maximum value of |h₃| = 0.09 when we also take the higher-order moments of the LOSVD into account.

The stellar kinematics maps we measure for m50001, by increasing the noise to S/N = 5, 10, 15, and 20 pixel⁻¹, are shown in Fig. 11. We clearly identify the 2σ peak for the galaxy model with S/N = 20 pixel⁻¹ when the two offset peaks of velocity dispersion are nicely outlined by the contour level at σ = 60 km s⁻¹, even though the increase in σ is small relatively to σ_pro. The maximum value of |h₃| is ∼0.1 in all cases, and the h₃ signature extends as far as the counter-rotating component, out to surface brightness contribution f_cr = 0.05. Even at low S/N, the anti-correlation between h₃ and v is more clearly discernible than the 2σ peak. We point out that our analysis considers only pure disk structures. Thus, this large-scale structure in the h₃ map may be misleading when other structural components are present.

Fig. 11. Stellar kinematic maps of v, σ, h3, and h4 (from left to right) of the LOSVD for the model m50001 derived from data with S/N = 5, 10, 15, 20 pixel−1 (from top to bottom). The inner and outer dotted ellipses show where fcr = 0.10 and 0.05, respectively. The solid contours in the σ and h3 maps correspond to σ = 60 km s−1 and |h3| = 0.05, respectively.

In this respect, when analyzing real galaxies, where the presence of intrinsic irregular structures or inhomogeneities due to instrumental characteristics make the data analysis more challenging, an S/N treatment like the one suggested by González-Gaitán et al. (2019) will improve our results.

6. Conclusions

The signature of two large-scale counter-rotating stellar disks appears in a galaxy spectrum with a large variety of LOSVDs that are broadened, strongly asymmetric, or even bimodal. This pushes to the limit the most widely used method for extracting the stellar kinematics, namely the single GH parameterization, which nevertheless provides important clues to the presence of multiple kinematic components in a galaxy.

In this work we tested the LOSVD recovery for a set of mock spectra adopting the GH parametric fit. Though the residuals decrease with increasing number of GH moments and S/N, we find that the fit to strongly non-Gaussian or bimodal spectral lines is not appropriate. Even allowing |h₃| and |h₄| to be larger than 0.3, the best-fit solution is poor, especially when the two counter-rotating components have large velocity separation and one of them contributes less than 30% to the luminosity. However, the GH parameterization is useful for revealing counter-rotating signatures in spatially resolved kinematics.

The extraction of the stellar LOSVD using the GH series lead to the conclusion that by adding a counter-rotating component to a dominant prograde disk the amplitude of the observed rotation curve is lowered, while the observed velocity dispersion is increased. In particular, our analysis confirms that the stronger evidence of the presence of two counter-rotating stellar disks is a 2σ feature along the major axis of the galaxy. The size, shape, and slope of the 2σ peaks strongly depend on the velocity separation and contribution of the two counter-rotating disks. This signature is highlighted even by measuring only the lower-order moments v and σ of the GH truncated series.

Including the higher-order moments to the fit, we find a change in the observed mean velocity, which becomes closer to the velocity of the brighter component compared to that measured by fitting v and σ only. Concerning the velocity dispersion, we find that the 2σ feature is not always present. It disappears in favor of a new kinematic features, called the 2σ drop, which is an artifact of the adopted parameterization when fitting high S/N data. In other words, the code detects only one component, which is the dominant one. Thus, the best-fit solutions give velocities and velocity dispersions closer to those of the single component. However, the contribution of the secondary disk is still revealed by the higher-order moments. By adding noise to the spectra, the 2σ peak is recovered, while the higher-order signatures become less prominent.

At low S/N regimes and for faint counter-rotating stellar component, the 2σ feature may be unclear. Interestingly, we find that the local anti-correlation of h₃ with v is always a useful indicator of counter-rotation as long as we consider only the disk structure. In addition, the h₃ parameter provides a criterion to distinguish between strong and weak counter-rotation. When the asymmetries in the LOSVD become relevant, then the |h₃| assumes values larger than 0.2. Contrary to the 2σ feature, the large-scale h₃ signature weakly fades with the decrease in S/N. Independently of the S/N, all the signatures become fainter at lower viewing angles as an effect of the smaller projected velocity separation between the two counter-rotating stellar components. Hence, highly inclined galaxies are favored for the application of the kinematic diagnostics, limiting the statistics of counter-rotating stellar disks.

A fundamental step in improving the statistics of counter-rotating galaxies is to detect faint counter-rotating stellar disks while looking for the 2σ feature, but also analyzing the signatures in the higher-order moment maps. The counter-rotation diagnostics confirmed in this work can be used as a guideline to search for counter-rotating signatures in the kinematic maps of galaxies stored in the MUSE data archive. This research will increase the number of optimal candidates to test and fine-tune the state-of-the-art decomposition techniques, in particular those based on the dynamical modeling using kinematic constrains. The study of large-scale counter-rotating galaxies will help to determine the relative importance of the different formation scenarios probing the galaxy assembly models in the cosmological context.

Present-day research on extended counter-rotation is limited by the small number of objects known to host this phenomenon. The IFU data available in the archives still do not have the necessary spatial extent and spatial resolution to perform complete statistics. In the near future, the IFU mounted on next-generation telescopes will allow us to investigate counter-rotators or, more generally, decoupled components (Corsini et al. 2012), at large distances and therefore at an earlier stage of their evolution. For instance, the MCAO Assisted Visible Imager and Spectrograph (MAVIS) instrument at the Very Large Telescope (McDermid et al. 2020) with a 50 mas spaxels (FoV of 5 × 7.2 arcsec²) and the low-resolution spectrograph will allow us to investigate counter-rotation in galaxies at z = 0.4−0.6. The High Angular Resolution Monolithic Optical and Near-infrared Integral field spectrograph (HARMONI) at the European-Extremely Large Telescope will also provide a similar configuration.

Acknowledgments

We thank the anonymous referee for the constructive report that helped us to improve the paper. We thank A. Moiseev and W. W. Zeilinger for their valuable comments. MR is grateful to P. Jethwa, G. van de Ven, and J. Falcón-Barroso for inspiring discussions over the course of this project, and thanks the Department of Astrophysics, University of Vienna, for hospitality during the preparation of this paper. MR, AP, EMC, and EDB are supported by MIUR grant PRIN 2017 20173ML3WW_001 and Padua University grants DOR1885254/18, DOR1935272/19, and DOR2013080/20. VPD is supported by STFC Consolidated grant ST/R000786/1.

References

All Tables

All Figures

Fig. 1. Best fit (red lines) of some mock spectra (black lines) with their residuals (blue lines) obtained using the GH parameterization for three sets of kinematic measurements: (1) lower-order moments only (v, σ, left panels); (2) four moments (v, σ, h3, h4, middle panels); (3) six moments (v, σ, h3, h4, h5, h6, right panels). Each panel shows mock spectra with decreasing velocity separation from top to bottom as labeled, while the constant contribution of the counter-rotating component increases from top to bottom (fcr = 0.17, 0.23, and 0.50, respectively).

In the text

	Fig. 2. As in Fig. 1, but shown are the best fit to mock spectra with fcr = 0.50 and S/N = 50 pixel−1 obtained using four GH moments.
In the text

Fig. 3. Comparison of the recovered stellar kinematic parameters in the case of mock spectra without noise (top four panels) and with S/N = 50 pixel−1 (bottom four panels). In both cases the mean velocity v (upper left panel), velocity dispersion σ (upper right panel), and higher-order moments h3 (bottom left panel) and h4 (bottom right panel) of the LOSVD of two counter-rotating components are known as a function of their input velocity separation ΔV. The points correspond to the measured values for mock spectra with fcr = 0.17 (red), 0.23 (orange), 0.33 (light-blue), and 0.50 (black), respectively. The points are connected to highlight the trend. The black solid lines indicate the input value of the mean velocity (v = Vpro), velocity dispersion (σ = σpro), and h3 = h4 = 0 as a function of ΔV. The black dashed lines indicate the default limits to the higher-order moments.

In the text

Fig. 4. Stellar kinematic maps (top panels) and major-axis radial profiles (solid lines with gray error bars, bottom panels) of the mean velocity v (upper left panels), velocity dispersion σ (upper right panels), and higher-order moments h3 (lower left panels) and h4 (lower right panels) of the LOSVD for IC 719 derived from MUSE data. The green dashed lines show the major-axis radial profiles of the kinematic parameters obtained with our model.

In the text

Fig. 5. Stellar velocity (top panels) and velocity dispersion (bottom panels) maps of models belonging to subset 1 (top row), 2 (middle row), and 3 (bottom row) with an inclination of 50°. The white ellipse shows where fcr = 0.50; this level is not reached in all models. The black dashed lines indicate the velocity separation Δv = 100 km s−1. In the bottom left corner of each velocity map we give the minimum and maximum values of v. These kinematic measurements are obtained fitting only v and σ.

In the text

	Fig. 6. Comparison of stellar velocity and velocity dispersion maps (left and right panel of each pair of plots, respectively) of models mii004 (top row), mii008 (middle row), and mii012 (bottom row) with the three adopted inclinations of 40° (right column), 50° (middle column), and 70° (left column), respectively.
In the text

	Fig. 7. Ratio of the peak velocity dispersion σpeak to the input velocity dispersion σpro as a function of the ratio of the total luminosity of the counter-rotating component Lcr to that of the model LT for subset 1 (filled circles). The yellow, blue, and red lines show the trend for inclination i = 40° ,50°, and 70°, respectively.
In the text

Fig. 8. Stellar kinematic maps of v, σ, h3, and h4 (from left to right) of the LOSVD for the model m50008 derived from data without noise (upper panels) and with S/N = 5 pixel−1 (lower panels). The inner and outer dotted ellipses show where fcr = 0.50 and 0.30, respectively. The solid vertical lines give the position of the 2σ drop detected along the major axis radial profile of the velocity dispersion.

In the text

	Fig. 9. Local correlation of h3 − v/σ (top panels) and radial trend of h3 as function of x/re (bottom panels) for the first four models of subset 1 (from left to right) with inclination i = 40°, 50°, and 70° (from top to bottom, respectively). The filled circles are the data points that are connected through a solid line to better illustrate the trend.
In the text

	Fig. 10. Stellar kinematic maps of v, σ, h3, and h4 (from left to right) of the LOSVD for the simulated galaxy with (sim0-11i60/80, top panels) and without (sim6-11i60/80, bottom panels) counter-rotating stars.
In the text

	Fig. 11. Stellar kinematic maps of v, σ, h3, and h4 (from left to right) of the LOSVD for the model m50001 derived from data with S/N = 5, 10, 15, 20 pixel−1 (from top to bottom). The inner and outer dotted ellipses show where fcr = 0.10 and 0.05, respectively. The solid contours in the σ and h3 maps correspond to σ = 60 km s−1 and |h3| = 0.05, respectively.
In the text