Uncovering the formation of the counter-rotating stellar disks in SDSS J074834.64+444117.8

Using the integral field spectroscopic data from Mapping Nearby Galaxies at Apache Point Observatory survey, we study the kinematics and stellar population properties of the two counter-rotating stellar disks in a nearby galaxy SDSS J074834.64+444117.8. We disentangle the two stellar disks by three methods, including CaII $\lambda$8542 double Gaussian fit, pPXF spectral decomposition, and orbit-based dynamical model. These three different methods give consistent stellar kinematics. The pPXF spectral decomposition provides the spectra of two stellar disks, with one being more luminous across the whole galaxy named primary disk, and the other named secondary disk. The primary disk is counter-rotating with ionized gas, while the secondary disk is co-rotating with ionized gas. The secondary disk has younger stellar population and poorer stellar metallicity than the primary disk. We estimate the stellar mass ratio between the primary and secondary disks to be $\sim$5.2. The DESI $g$, $r$, $z$ color image doesn't show any merger remnant feature in this galaxy. These findings support a scenario that the counter-rotating stellar disks in SDSS J074834.64+444117.8 formed through gas accretion from the cosmic web or a gas-rich companion.


INTRODUCTION
Counter-rotating stellar disks are characterized by the presence of two stellar disks that are co-spatial in one galaxy but rotating in opposite directions.Such a feature was first discovered in NGC 4550 by absorption line analyses (Rubin, Graham, & Kenney 1992;Rix et al. 1992).Moreover, 2σ feature was proposed by Krajnović et al. (2011) to describe the observation of off-center but symmetric peaks along the major axis in the stellar velocity dispersion field of a galaxy with counter-rotating stellar disks, where the co-existing blueshifted and redshifted absorption components contribute to the broadening of absorption lines.
Individual galaxies hosting counter-rotating stellar disks have been studied, including NGC 5719 (Coccato et al. 2011), NGC 4550 (Johnston et al. 2013), NGC 3593 (Coccato et al. 2013), NGC 4138 (Pizzella et al. 2014), NGC 448 (Katkov et al. 2016), NGC 5102 (Mitzkus, Cappellari, & Walcher 2017), and IC 719 ⋆ E-mail: chenym@nju.edu.cn† E-mail: myang@shao.ac.cn (Pizzella et al. 2018).Using the spectral decomposition method, it is possible to disentangle the properties of the co-spatial stellar disks.Johnston et al. (2013) decomposed the spectra along the major axis of S0 galaxy NGC 4550, finding that the stellar disk that co-rotates with the gas disk has younger stellar population than the counter-rotating one, which indicates the co-rotating stellar disk lately forming from the accreted gas.Coccato et al. (2011) applied the decomposition method to the spectra of spiral galaxy NGC 5719, and revealed that the co-rotating stellar disk not only has younger stellar population, but also has poorer stellar metallicity.Considering that NGC 5719 is interacting with a companion NGC 5713, Coccato et al. (2011) suggested these findings as the accreted gas fueling the in-situ formation of a new stellar disk.
With the development of integral field spectroscopic (IFS) surveys, studies on samples of counter-rotating stellar disks became feasible, which revealed the dependence of formation scenarios on numerous physical factors.Bevacqua, Cappellari, & Pellegrini (2022) selected 64 galaxies with counter-rotating stellar disks from ∼4,000 galaxies in the MaNGA survey Data Release 16, and found that in most cases the younger stellar disk co-rotates with the arXiv:2401.11179v1[astro-ph.GA] 20 Jan 2024 gas disk, which supported the formation scenario of gas accretion.However, the gas co-rotates with the older stellar disk in two galaxies, which implies a disk galaxy merging with a gas-poor galaxy with younger stellar population in a retrograde orbit.Moreover, Bao et al. (2022) collected a sample of 101 galaxies with counterrotating stellar disks from 9,456 galaxies in the MaNGA survey Product Launch 10, divided the sample into four types based on the stellar and gas kinematics, and proposed different formation scenarios for different types.They suggested that the key factors in the formation of counter-rotating stellar disks are the abundance of pre-existing gas in the progenitor and the efficiency of angular momentum consumption between the pre-existing and external gas.Given the complexity of formation mechanisms, it is necessary to adopt reliable methods to disentangle two stellar disks, and compare properties between them as completely as possible to uncover the formation of counter-rotating stellar disks.
Orbit-based dynamical model provides a physical method to quantify the structure and dynamics of different galaxy components based on their orbit structure distribution.This method has been applied to galaxies to explore the relation between galaxy components and galaxy assembly history in several IFS survey, such as CALIFA (Zhu et al. 2018), MaNGA (Jin et al. 2020), SAMI (Santucci et al. 2022) and MUSE (Ding et al. 2023).And the method has recently been updated to include bar properly (Tahmasebzadeh et al. 2022).However, it has not been applied to the galaxies with counter-rotating stellar disks, since the line-of-sight velocity distribution (LOSVD) of these galaxies cannot be well described by the widely-used Gaussian-Hermite expansion of single dynamical component (Rubino et al. 2021).The solution to this problem relies on the non-parametric LOSVDs which capture the multiple dynamical components (e.g.Katkov et al. 2011;Falcón-Barroso & Martig 2021), providing a new prospect of modelling counterrotating stellar disks with orbit-based dynamical models and public DYNAMITE package (van den Bosch et al. 2008;Jethwa et al. 2020;Thater et al. 2022).
In this paper, we focus on a nearby (z ∼ 0.02) galaxy SDSS J074834.64+444117.8 (hereafter SDSS J0748+4441).We disentangle the two stellar disks by three methods, including Ca II λ8542 double Gaussian fit, pPXF spectral decomposition, and orbit-based dynamical model.The involved data are presented in Section 2. The spatially resolved galaxy parameters, the disentangle results, as well as the properties of the two stellar disks are presented in Section 3. Finally, we discuss the formation scenario for counter-rotating stellar disks in Section 4.

THE DATA
MaNGA is an integral field spectroscopic (IFS) survey conducted as a part of the fourth-generation Sloan Digital Sky Survey (SDSS-IV; Bundy et al. 2015;Law et al. 2016).MaNGA survey provides a representative sample of 10,010 unique galaxies with a flat stellar mass distribution in the range of 10 9 − 10 11 M⊙ and redshift in the range of 0.01 < z < 0.15 (Blanton et al. 2017).This survey utilizes the Baryon Oscillation Spectroscopic Survey (BOSS) spectrographs (Smee et al. 2013) on the 2.5-m Sloan Foundation Telescope (Gunn et al. 2006).The dual-channel BOSS spectrographs (Smee et al. 2013) provide simultaneous wavelength coverage from 3,600 to 10,000 Å.
For each target, the MaNGA Data Reduction Pipeline (DRP; Law et al. 2016) produced sky-subtracted spectrophotometrically calibrated spectra, and generated three-dimensional datacube con-taining spatially resolved spectra.The wavelength calibration of the MaNGA data is accurate to 5 km s −1 (rms), with a median spectral resolution of 72 km s −1 (R ∼2,000).Moreover, the MaNGA Data Analysis Pipeline (DAP; Westfall et al. 2019) provided measurements with pixel size of 0.5 ′′ /spaxel (∼211 pc/spaxel), including stellar kinematics, ionized gas kinematics, emission line flux and equivalent width, as well as spectral indices such as 4000 Å break indicating the light-weighted stellar population age.The gas kinematics in this paper is traced by Hα emission, while all the emission line centers are tied together in the velocity space in MaNGA DAP.Besides, we extract the stellar mass in each spaxel from Pipe3D (Sánchez et al. 2016), which used a spectral fitting tool FIT3D to analyse the physical properties of stellar populations of a galaxy.
SDSS J0748+4441 was chosen from a sample of 101 galaxies hosting counter-rotating stellar disks (CRDs; Bao et al. 2022) since two distinguishable absorption components are clearly shown in one of the Ca II triplet lines at λ8542 which is free from skylines.The MaNGA spectra have high enough signal-to-noise ratio (SNR) for disentangling the absorptions of two stellar disks even on the outskirts (∼1.5 Re).The SDSS g, r, i color image of this galaxy is displayed in Figure 1(a).The global stellar mass (M⋆) and attenuation corrected star formation rate (SFR) from the MaNGA DRP and DAP indicate that SDSS J0748+4441 locates at the green-valley region.

Spatially resolved properties
Figure 1 displays the kinematics of stars and ionized gas in SDSS J0748+4441.Figures 1(b) and 1(c) display the stellar velocity and velocity dispersion fields for spaxels with median spectral SNR per spaxel higher than 3, where the black dashed line shows the photometric major axis (hereafter major axis).The two black stars in Figure 1(c) mark the 2σ peaks along the major axis, where the maximum stellar velocity dispersion locates.The presence of two regions with enhanced stellar velocity dispersion along the major axis originates from the existence of two CRDs.Meanwhile, the stellar velocity field in Figure 1(b) following a regular pattern indicates that it is dominated by the more luminous stellar disk across the whole galaxy.
Figures 1(e) and 1(f) display the ionized gas velocity and velocity dispersion fields for spaxels with Hα emission line SNR higher than 3.The gas is regularly rotating as displayed in Figure 1(e), with low gas velocity dispersion around the 2σ peaks (black stars) in Figure 1(f).We fit the gas velocity field using the KINEMETRY package (Krajnović et al. 2006) to derive the kinematic major axis, which is shown by the grey dashed lines in Figures 1(e) and 1(f).Combining the stellar and gas velocity fields in Figures 1(b) and 1(e), we conclude that the more luminous stellar disk and the gas disk in SDSS J0748+4441 are counter-rotating.Referring to the kinematic classification in Bao et al. (2022), this galaxy belongs to Type 2b, where the old stellar disk outshines the newly formed one and counter-rotates with gas disk.In Section 3.2, we will disentangle and compare the properties of two stellar disks.
Figure 2(a) displays a map of stellar mass surface density (Σ⋆), for spaxels with median spectral SNR per spaxel higher than 3.The Σ⋆ is defined as the stellar mass of each spaxel divided by the physical size of the spaxel.Figure 2(b) displays Σ⋆ as a function of radius along major axis.Σ⋆ monotonically decreases with increasing radius, but the gradients in the central region and on the outskirts are obviously different.The grey area with |R| ⩽ 2.8 ′′ separates the central region with high surface density from the outskirts with low surface density.We fit two Gaussian functions to the data points inside and outside the grey area.The blue solid profile shows the Gaussian model for the Σ⋆ gradient with |R| ⩽ 2.8 ′′ , while the blue dashed profile shows the Gaussian model for the data points with |R| > 2.8 ′′ .It is clear that the central region and outskirts follow different Σ⋆ distributions, dominated by bulge and disk, respectively.We mark a red dashed ellipse with a half major axis of |R| = 2.8 ′′ on the Σ⋆ map, which basically covers the densest bulge region.
Mapping diagnostic emission line ratios across a galaxy gives the information on the ionization state distribution of gas. Figure 3(a) displays the [N II] BPT diagram (Baldwin, Phillips, & Terlevich 1981) for spaxels with Hβ, [O III]λ5007, Hα, and [N II]λ6583 emission line SNRs higher than 3.The black solid curve separates the star forming and composite regions (Kauffmann et al. 2003), while the black dashed curve is the demarcation between the composite and AGN regions (Kewley et al. 2001).The dots show the line ratios measured in different spaxels, and they are color-coded by the distances of corresponding spaxels to the black solid curve.The star formation rate describes the ongoing activity of star formation that occurred within the last 10 6−7 years.We derive the global SFR of SDSS J0748+4441, which is ∼6.7×10 −2 M⊙ yr −1 , based on attenuation corrected Hα flux from MaNGA DAP. Figure 3(c) displays a map of Hα equivalent width (EQW), which is a good indicator of the specific star formation rate in the star forming region.The black contours in Figure 3(c) outline three regions with enhanced Hα EQW, which are also marked in Figure 3(b).It is clear that the Hα EQW enhanced regions are dominated by star formation, indicating a higher sSFR in these areas.It is interesting that the two black stars marking the 2σ peaks locate within the contours with the enhanced Hα EQW.
Figure 3(d) displays a map of 4000 Å break (Dn4000), which is influenced by the intensity of star formation in Gyr-timescale.Along the major axis, Dn4000 presents the highest value of ∼1.8 in the bulge region (inside the red dashed ellipse), while it has a value of ∼1.5 in the disk region, including the 2σ peaks (black stars) and the enhanced Hα EQW regions (black contours).The homologous Dn4000 distribution but clumpy Hα EQW enhancement in the disk region can be explained by the different star formation timescales represented by Hα emission and Dn4000.

Double Gaussian fit on Ca II λ8542
SDSS J0748+4441 shows obvious asymmetric absorption structure in one of the Ca II triplet lines at λ8542, which provides us an opportunity to disentangle the two stellar disks without model dependence.The other two of Ca II triplet lines at λλ8498, 8662 are contaminated by the skylines.
Figure 4(a) displays the stellar velocity dispersion field given by MaNGA DAP, in which the two black rings outline the regions around the 2σ peaks (black stars), with the top one named 'Region A' and the bottom one named 'Region B'.To intuitively display the two absorption components in the Ca II λ8542 absorption line, we extract and stack the spectra within Regions A and B, respectively.The Ca II λ8542 absorption lines in the stacked spectra of Regions A and B are shown in the Figures 4(b) and 4(c), where we use the python-based tool curve fit to conduct double Gaussian fit on them.The Gaussian models with stronger absorption are in red, the Gaussian models with weaker absorption are in blue, and the bestfit models are in green which is the combination of blue and red components.

pPXF spectral decomposition
Taking advantage of the different velocities and velocity dispersions of two CRDs, we decompose their contributions to the observed spectrum in each spaxel along the major axis, using the penalized pixel fitting (pPXF; Cappellari 2017) code.The code builds two synthetic templates (one for each stellar disk) as linear combination of stellar spectra from the MILES library, convolves the templates with two Gaussian LOSVDs with different velocities and velocity dispersions, and ensures the combination of two convolved templates fitted to the observed spectrum by χ 2 minimization.
Figure 4(d) displays an example of pPXF two kinematic component fit.The observed spectrum shown in black is extracted from the spaxel marked by the black filled star in Figure 4(a).The red and blue spectra show the optimal models for the two stellar disks.The red spectrum corresponds to the stellar disk with higher flux (hereafter primary disk), and the blue spectrum corresponds to the stellar disk with lower flux (hereafter secondary disk).The green spectrum is the combination of the primary and secondary components, it is our best-fit model.The insert drawing of Figure 4(d) highlights the spectra in [4800, 5300] Å wavelength range.The four grey dashed lines mark the rest-frame wavelength of Hβ and Mg I triplet lines, where all the absorption profiles are well fitted.
Figure 5 displays the line-of-sight velocities of stars and gas as functions of radii.The line-of-sight velocities of the two stellar disks are obtained by the pPXF two kinematic component fit on the continuum and absorption lines for spaxels along the major axis.The red and blue circles in Figure 5 represent the line-ofsight velocities of the primary and secondary disks, respectively.The line-of-sight velocity of the primary disk is overall lower than that of the secondary disk.We also obtain stellar velocity dispersion by the pPXF two kinematic component fit.The stellar velocity dispersions of both primary and secondary disks follow flat distribution along the major axis, with comparable values of 84 km s −1 and 79 km s −1 , respectively.The grey area marks the bulge region, where the stellar kinematics of primary disk is dominated by velocity dispersion, resulting in larger decomposition error.
We collect the velocity of the ionized gas traced by Hα from the MaNGA DAP file (Westfall et al. 2019), which are represented by the grey squares in Figure 5.The line-of-sight velocities of the secondary disk (blue circles) and the gas disk (grey squares) are totally consistent, indicating the secondary disk is co-rotating with the gas disk.Meanwhile, the primary disk (red circles) is counterrotating with the gas disk.The fit on Ca II λ8542 lines in Figures 4(b) and 4(c) also provides velocities of two stellar disks in Regions A and B, which are shown as the crosses in Figure 5, with the primary and secondary disks in red and blue, respectively.The line-ofsight velocities of the two stellar disks measured from Ca II λ8542 absorption lines match well with the results from the pPXF spectral decomposition.

Orbit-based dynamical model
The Schwarzschild orbit-superposition method (Schwarzschild 1979) is a powerful orbit-based dynamical modelling technique to reveal orbit structures in a galaxy.This method computes all possible stellar orbits in a given potential and assigns weight to each orbit to recover galaxy kinematics, including widely used parametric LOSVDs (e.g., Gaussian-Hermite expansion with velocity and velocity dispersion) in ordinary galaxies and non-parametric LOSVDs which are superior in describing the counter-rotating spectral features generated by two stellar disks.Compared with the spectral decomposition method, it makes use of spatial information contained in all IFS spectra of a galaxy and sidesteps the high SNR requirement when decomposing two stellar components directly from the spectra.In this section, we confirm the kinematics of the CRDs in SDSS J0748+4441 with the orbit-based dynamical models using the non-parametric LOSVDs of this galaxy.
We first extract the non-parametric LOSVDs of this galaxy from the MaNGA datacube with the BAYES-LOSVD package 1 , which employs Bayesian inference to obtain LOSVDs and the corresponding uncertainties.This package improves its efficiency by adopting a principal component analysis to reduce the requirement of stellar templates, and provides several prior options to regularize the output LOSVDs.We refer readers to Falcón-Barroso & Martig (2021) for more details.For SDSS J0748+4441, the MaNGA spectra are first Voronoi binned to SNR > 50 using the median SNR within the fitting spectra range [4750, 5500] Å, and then fitted with the MILES stellar template library (Sánchez-Blázquez et al. 2006;Falcón-Barroso et al. 2011).The output LOSVDs have a velocity range [-600, 600] km/s with a step of 30 km/s (∼1/2 MaNGA instrumental resolution).As an example, we show the LOSVD extraction of bin 130 in Figure 6, which is the closest to the blackfilled star in Figure 4(a).
We then construct orbit-based dynamical models from the non-parametric LOSVDs using the DYNAMITE package (van den Bosch et al. 2008;Jethwa et al. 2020;Thater et al. 2022) and obtain the orbit distribution of this galaxy.We build the stellar mass distribution of this galaxy from its r-band image (Kaiser et al. 2002) with the MGE method (Emsellem, Monnet, & Bacon 1994;Cappellari 2002).Moreover, we correct for its mass-to-light ratio with the stellar mass obtained from the datacube of Pipe3D (Sánchez et al. 2018), and introduce a variable to indicate the scale of the total stellar mass.We adopt the NFW dark matter profile described with virial mass M200 and concentration c.M200 is defined as the enclosed mass within the virial radius r200, where the average density is 200 times the critical density (ρcrit = 1.37 × 10 −7 M⊙ pc −3 ), and c is defined as the ratio of r200 to the scale radius of the NFW dark matter profile.We also include a central black hole with a fixed mass of 10 6 M⊙.More details will be described in Yang et al. (in preparation).We show the fit along the major axis for SDSS J0748+4441 in Figure 7, which includes bin 130 mentioned above.
We further characterize different orbit types with the circularity λz = Lz/rVc. (1) For each orbit, Lz is the mean angular momentum along the short axis, r is the average radius in equal time step for each orbit, and Vc is the mean circular velocity (see Equation 9in Zhu et al. 2018 1 https://github.com/jfalconbarroso/BAYES-LOSVDfor more detailed mathematical definitions).We divide the orbits into three stellar components by clear circularity distinctions at λz = ±0.16,shown as two blue dashed lines in the first panel of Figure 8.The secondary disk is comprised of the co-rotating component classified as λz > 0.16, while the primary disk is comprised of the non-rotating and counter-rotating components classified as −0.16 < λz < 0.16 and λz < −0.16.In the second panel of Figure 8, we display the line-of-sight orbital velocities as functions of radii along the top half of the major axis for the two stellar disks, with the primary and secondary disks represented by the orange and purple squares.The line-of-sight velocities obtained by the orbit-based model are consistent with those obtained by the spectral decomposition, with the primary and secondary disks represented by red and blue circles in the second panel of Figure 8, respectively.

Properties of two stellar disks
Based on the disentangled spectra of the two stellar disks in SDSS J0748+4441, we compare the properties between them.We measure the Lick indices Hβ, Mgb, Fe5270, and Fe5335 (Worthey et al. 1994), for the primary and secondary components.Hβ can be used as an indicator of stellar population age.
We also calculate the combined magnesium-iron index [MgFe] ′ as Mgb • (0.72 • Fe5270 + 0.28 • Fe5335), which is independent of α-enhancement and can be an effective indicator of stellar metallicity (Thomas & Maraston 2003).Figure 9(a) displays the Hβ indices as functions of [MgFe] ′ indices along the major axis, with the primary and secondary disks in red and blue, respectively.The model-grid is the prediction from single stellar population models given by Thomas, Maraston, & Johansson (2011).Comparing the data points with the model-grid, the secondary disk has younger stellar population and poorer stellar metallicity than the primary disk.
Figure 9(b) displays the Dn4000 indices as functions of radii along the major axis.The circles represent the results from the pPXF two kinematic component fit, with the primary and secondary disks in red and blue, respectively.The stellar population of the primary disk is old, with Dn4000 of ∼1.8.Meanwhile, the stellar population of the secondary disk is young in the disk region, presenting a constant Dn4000 of ∼1.4.The Dn4000 of the secondary disk increases in the central region, where the emission is bulge-dominated.
Figure 9(c) displays the stellar metallicity as functions of radii along the major axis, which is estimated by comparing the position of each circle with the model-grid in Figure 9(a).The stellar metallicity of the primary disk is rich, and follows a constant distribution.Meanwhile, the stellar metallicity of the secondary disk is poor, and follows a negative gradient with metallicity decreasing on the outskirts.Given the bulge-dominated feature in the central region displayed in Figure 9(b), we don't overinterpret the negative gradient here.Moreover, we estimate the stellar population age for two stellar disks in Figure 9(a), and find the stellar population of the primary disk being older than that of the secondary disk, which is consistent with the Dn4000 indices in Figure 9(b).

FORMATION OF COUNTER-ROTATING STELLAR DISKS IN SDSS J0748+4441
Early studies proposed that internal structures such as bulge and bar (Evans & Collett 1994) are responsible for the formation of CRDs.
However, these scenarios failed to explain the presence of different stellar populations of two stellar disks (Coccato et al. 2013), which instead can be the results of external processes such as mergers and gas accretion (Corsini 2014).Major merger is known to be a destructive event that heats the system significantly, destroys the morphology of galaxies, and results in the formation of an elliptical galaxy.On the other hand, gas-rich minor merger and gas accretion are milder processes delivering the external gas that is counterrotating with the pre-existing gas, and have weaker influence on the galaxy morphology (Bournaud, Jog, & Combes 2005).
With the development of long-slit spectroscopic and IFS instruments, studies on both individual galaxies (Coccato et al. 2011;Johnston et al. 2013) and galaxy samples (Bevacqua, Cappellari, & Pellegrini 2022;Bao et al. 2022) confirmed that the external processes including mergers and gas accretion contribute to the formation of CRDs.The common point in these processes is that the lately forming stellar disk inherits the angular momentum from the external gas, which is counter-rotating with the pre-existing stellar disk.During the formation of CRDs, once the external gas with opposite direction is abundant enough to dominate the gas rotation, the misalignment between gas and stars will be observed in galactic scale (Chen et al. 2016;Xu et al. 2022).Davis et al. (2011) found that the gas and stars are misaligned in ∼42 percent early type galaxies in ATLAS 3D .However, assuming mergers as the only source of misalignment in simulation, Lagos et al. (2015) obtained only ∼2-5 percent early type galaxies with misaligned gas and stars.After adding gas accretion in their simulation, the misalignment between gas and stars was found in ∼46 percent early type galaxies.These studies proved that gas accretion is a key process of bringing the external gas with opposite direction.
Based on the pPXF two kinematic component fit as displayed in Figure 4(d), we calculate the flux contributions of two stellar disks in each spaxel.Figure 10(a) displays the flux contributions as functions of radii along the major axis, with the primary and secondary disks in red and blue, respectively.The primary disk is more luminous than the secondary disk across the whole galaxy.The average flux ratio between the primary and secondary disks is ∼1.8.We also estimate the masses of two stellar disks based on the decomposed spectra.For each disk model, we convolve the model spectrum with the SDSS u, g, r, i, z filters to get the AB magnitudes.The stellar mass-to-light ratio is then obtained by comparing the u, g, r, i, z magnitudes of each stellar disk with mimic galaxies generated in Chen et al. (2012).We refer readers to Chen et al. ( 2012) for more details of the methodology.Figure 10(b) displays the stellar masses of two stellar disks as functions of radii, with the primary and secondary disks in red and blue, respectively.The primary disk is more massive than the secondary disk across the whole galaxy, with mass ratio of ∼5.2.
A gas-rich minor merger in retrograde orbit can provide abundant gas that is counter-rotating with the pre-existing stars, and fuel the formation of secondary disk.Ji, Peirani, & Yi (2014) simulated mergers with different mass ratios, and found that the merger remnant features can be observed for ∼5.8 Gyrs at r-band magnitude limit of ∼25 mag arcsec −2 for a merger with a mass ratio of 6.The stellar population age of the secondary disk in SDSS J0748+4441 is ∼2 Gyrs in Figure 9(b), suggesting that the gas-rich minor merger should happen ∼2 Gyrs ago.On the other hand, the DESI g, r, z color image of SDSS J0748+4441 is displayed in Figure 1(d), which is much deeper than 25 mag arcsec −2 in r-band.We follow the method described in Li et al. (2021) to check whether there is any merger remnant feature, finding that it is an isolated galaxy.The gas-rich minor merger is inapplicable to explain the formation of CRDs in SDSS J0748+4441 due to the lack of merger remnant features.
We propose gas accretion in retrograde orbit as formation scenario for CRDs in SDSS J0748+4441.Given the stellar population age of the secondary disk in Figure 9(b), the gas accretion in this galaxy happened ∼2 Gyrs ago.The global stellar mass of this galaxy in the MaNGA DRP is ∼10 10 M⊙.Combining the mass ratio of ∼5.2 between two stellar disks (Figure 10b), accretion brought at least 1.6×10 9 M⊙ counter-rotating gas into SDSS J0748+4441.The gas accretion can originate from galaxy environment including the cosmic web (Katkov et al. 2016) or a gas-rich companion (Coccato et al. 2011).
Figure 3(b) displays the spatially resolved BPT diagram, where the color-codes are the same as Figure 3(a).The over-plotted red dashed ellipse is the same as that in Figure 2(a) with a half major axis of |R| = 2.8 ′′ .Regions inside the ellipse are dominated by AGN, while the outer regions including the 2σ peaks (black stars) are dominated by star formation.

Figure 1 .
Figure 1.Images and spatially resolved kinematics.(a) SDSS g, r, i color image.(b) & (c) Stellar velocity and velocity dispersion fields for spaxels with median spectral SNR per spaxel higher than 3.The black dashed lines show the photometric major axis in Figures 1 to 4. The black stars mark the 2σ peaks along the major axis in Figures 1 and 3. (d) DESI g, r, z color image.(e) & (f) Gas velocity and velocity dispersion fields for spaxels with Hα emission line SNR higher than 3.The grey dashed line shows the kinematic major axis.

Figure 2 .
Figure 2. Spatially resolved stellar mass surface density.(a) A map of stellar mass surface density for spaxels with median spectral SNR per spaxel higher than 3.The red dashed ellipse with a half major axis of |R| = 2.8 ′′ outlines the bulge region in Figures 2 and 3. (b) The black circles represent the stellar mass surface density at different radii along the major axis.The black vertical bar shows the stellar mass surface density ±1σ error at each radius.The grey area marks the bulge region with |R| ⩽ 2.8 ′′ in Figures 2, 5, 9 and 10.The blue solid profile shows the gaussian model for the data points inside the grey area.The blue dashed profile shows the gaussian model for the data points outside the grey area.

Figure 6 .
Figure 6.Spectra fit and LOSVD of bin 130 in orbit-based dynamical model generated by BAYES-LOSVD.Top left: the binning map with the position of bin 130 marked by the red id.Top right: the blue solid line and shadow are the LOSVD and ±1σ error extracted by fitting the MaNGA spectra with BAYES-LOSVD.Bottom: MaNGA spectra (black), fit (red) and residual (green) of bin 130.The emission lines covered by grey area are masked prior to fit.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Fitting LOSVD along the major axis in orbit-based dynamical model.First panel: the binning map showing the ids and positions (in arcsec) of all the bins shown in the following panels.Other panels: the blue solid line and shadow are the LOSVD and ±1σ error obtained with BAYES-LOSVD, the red dashed line are the best-fit model, and the black dotted line are the symmetrized data for comparison because the model is axis-symmetric.

Figure 10 .
Figure 10.Contributions of two stellar disks along the major axis.(a) & (b) The circles show the fluxes and stellar masses as functions of radii, with the primary and secondary disks in red and blue.