Probing the History of the Galaxy Assembly of the Counterrotating Disk Galaxy PGC 66551

Stellar counterrotation in disk galaxies directly relates to the complex phenomenon of the disk mass assembly believed to be driven by external processes, such as accretion and mergers. The study of such systems makes it possible to reveal the source of external accretion and establish the details of this process. In this paper, we investigate the galaxy PGC 66551 (MaNGA ID 1-179561), which hosts two large-scale counterrotating (CR) stellar disks identified in the Sloan Digital Sky Survey MaNGA data and then confirmed using deep follow-up spectroscopy with the 10 m Southern African Large Telescope. We measured the properties of ionized gas and stellar populations of both CR disks in PGC 66551. We found that the CR disk is compact, contains young stars with subsolar metallicity, and has a stellar mass of 5 × 109 M ⊙, which amounts to ≈20% of the galaxy’s total. Surprisingly, the main 8 Gyr old disk has a significantly lower metallicity of −0.8 dex than other CR galaxies. We developed a simple analytic model of the history of the metal enrichment, which we applied to PGC 66551 and constrained the parameters of the galactic outflow wind, and estimated the metallicity of the infalling gas that formed the CR disk to be −0.9... −0.5 dex. Our interpretation prefers a merger with a gas-rich satellite over cold accretion from a cosmic filament as a source of gas, which then formed the CR disk in PGC 66551.

1. INTRODUCTION One of the fundamental questions in galaxy physics is "how do galactic disks grow?"External gas accretion is one of the important processes governing disk growth.Many observational phenomena cannot be explained without gas accretion (see Combes 2014;Putman 2017, for a recent review).For example, star formation rates (SFRs) observed at high redshifts sug-gest a galaxy would run out of gas in a few gigayears at most.As a result, without external replenishment, most galaxies today would be red and dead, contradicting the substantial star formation (SF) activity observed in most low-redshift systems (Genzel et al. 2010;Kennicutt 1998).It also turns out that the total budget of cold (< 10 4 K) atomic gas present in galaxies does not decrease significantly as the stellar mass density evolves, which is expected to be due to the lack of accretion (Madau & Dickinson 2014;Neeleman et al. 2016).Accretion of gas onto a galaxy is also required by chemical evolution models to successfully reproduce the metallicity distribution in the Milky Way and other galaxies (Chiappini 2009;Gallazzi et al. 2005).Despite the importance of this process, there is a lot we do not know about it, for example, How much material can be acquired during accretion?How long does the accretion last and at what redshift does it typically occur?What are the sources of acquired material?And finally, how does this process shape galaxies?Furthermore, the direct detection as well as the clear consequences of the past accretion are very elusive.Hence, studies of galaxies with clear evidence of past accretion are very important.Some of the best examples of such systems are galaxies with two counterrotating (CR) stellar disks.
CR stars are thought to originate from the infall of external material, the sources of which might be (i) cosmological filaments supplying pristine (i.e., primodial or very-low metallicity) gas (Algorry et al. 2014) or (ii) minor mergers with gas-rich satellites (Thakar & Ryden 1996, 1998;Lu et al. 2021), followed by a burst of SF (Pizzella et al. 2004).For the new gas with misaligned angular momentum to have settled in the main plane of the galaxy, the preexisting gas must have been removed, for example, by means of strong active galactic nuclei (AGN) feedback or gas stripping in a dense environment (Starkenburg et al. 2019).Another possibility is that the infalling gas removes the preexisting interstellar medium (ISM) and then forms CR stars in situ (Khoperskov et al. 2021).Both scenarios will lead to different time lags in SF and age differences of stellar population components.
All galaxies accrete material in some stage of their evolution, but in most cases, we cannot distinguish between accreted and original material in the finally shaped galaxy because the angular momenta of the accreted and original gas are sufficiently similar so that we cannot kinematically separate these two components.However, this is not the case in CR galaxies, where they are sufficiently different to be able to be identified well after the merger occurred.Kinematical disentanglement of the two disks provides a unique opportunity to probe stellar populations formed from the accreted material and the preexisted stellar populations before the infall event.
The presence of two CR stellar disks can be determined by a symmetric double peak on the velocity dispersion map, the so-called two-σ feature (Krajnović et al. 2011;Rubino et al. 2021).Thanks to this simple diagnostics, numerous CR galaxies have recently been detected in the SDSS Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey data (Bao et al. 2022;Bevacqua et al. 2022), exhibiting a variety of kinematic properties.However, so far only about a dozen CR galaxies have been studied in detail so far using deep spectroscopic observations (Coccato et al. 2011(Coccato et al. , 2013(Coccato et al. , 2015;;Katkov et al. 2013Katkov et al. , 2016;;Johnston et al. 2013), and analyzed with a dedicated spectral decomposition technique (Coccato et al. 2011;Katkov et al. 2013) that allows one to determine the parameters of the stellar population of both the main and CR stellar disk.In this paper, we expand the sample of well-studied counter-rotating galaxies by the intriguing galaxy PGC 66551, a dominant system in a pair outside of any major cluster or group, whose stellar population properties are remarkably different from the other known cases.Throughout the paper, we have adopted a physical scale of 0.395 kpc arcsec −1 according to NASA/IPAC Extragalactic Database 1 for the redshift of z = 0.0195 and the Planck 2015 cosmology (Planck Collaboration et al. 2016).
In Section 2, we describe new observations and archival data.In Section 3, we analyze our spectroscopic data using different approaches.In Section 4 we present our simple model of the history of the metal enrichment.A discussion and summary are provided in Sections 5 and 6, respectively.

OBSERVATIONAL DATA 2.1. MaNGA integral field spectroscopy
The MaNGA survey (Bundy et al. 2015) survey contains Integral Field Unit (IFU) (Drory et al. 2015) observations on the Sloan 2.5 m telescope (Gunn et al. 2006) of some 10,000 galaxies in the nearby Universe covering a wide range of stellar mass and colors (Wake et al. 2017).Each galaxy was observed using 2 ′′ fibers packed in bundles that varied in diameter from 12 ′′ (19 fibers) to 32 ′′ (127 fibers) covering (1.5 − 2.5)× the effective diameter of the target.Fiber bundles feed the BOSS spectrographs (Smee et al. 2013) providing a wavelength coverage of λλ3,600-10,300 Å with a resolving power of R ∼ 2000, which corresponds to an instrumental dispersion of σ inst.≈ 75 km s −1 at 5100 Å (Law et al. 2016).To achieve full coverage, each galaxy was observed with three dithered exposures to fill in the gaps between the fibers in a bundle (Law et al. 2015;Yan et al. 2016).
We inspected data from the DR17 release and identified ≈60 galaxies with two distinct off-center peaks in the stellar velocity dispersion maps, indicative of two CR disks (Katkov I. et al. in preparation).During the inspection process, we intensively used the web-based interactive service for the visualization of MaNGA data (https://manga.voxastro.org2), which we are developing.The PGC 66551 (MaNGA ID 1-179561) galaxy has one of the best examples of a CR system identified in our sample and therefore was chosen for a separate detailed study.
The MaNGA Data Analysis Product (DAP) map of the stellar velocity (Figure 1b) clearly shows that the  ± 2 • .Taking into account the limited MaNGA field-of-view, these estimates appear to be consistent within ≈ 1σ.The similarity of the P A values supports the idea that the gas is settled and rotates in nearly the same plane as the stellar disk.However, it is not possible to directly compare the rotation of the gaseous disk with that of the stellar one, since the former is a combination of two stellar disks.

SALT Long-slit spectroscopy
To study the properties of CR stars across the entire galaxy disk at a higher spectral resolution and overcome the limitations of the MaNGA field-of-view, we carried out deep long-slit spectroscopic observations of PGC 66551.The observations were collected with the Robert Stobie Spectrograph (RSS; Burgh et al. 2003;Kobulnicky et al. 2003) of the Southern African Large Telescope (SALT, Buckley et al. 2006;O'Donoghue et al. 2006) during the period May 19-31, 2020, in which seven pointings of 2 × 1200s exposures each were made, totaling to 4h40m scientific exposure on target.We used two configurations of the RSS spectrograph.We spent 4 hours of total exposure to obtain a deep stellar contin-3 https://photutils.readthedocs.io/uum spectrum using the PG2300 grating with a 1.25 ′′wide slit, providing a spectral resolution of R ≈ 3200 in the spectral range of 4500 − 5500 Å.The remaining 40 minutes were used to trace emission lines with a lower resolution of R ≈ 2300 in the 4865 − 6900 Å range using a PG1300 grating and a 1 ′′ wide slit.The slit was aligned along the major isophotal axis of the galaxy with P A = 285 • .The primary, preliminary and long-slit reduction of the spectra was carried out using the procedure described in Kniazev (2022).We also corrected the science spectra for scattered light in the instrument using the approach as presented in Katkov et al. (2019).Reference arc spectra were used to trace the instrumental resolution as a function of wavelength, which is an essential step in the subsequent modeling of spectra.The reduced spectra for both spectral configurations have a spatial scale of 0.25 ′′ pixel −1 , and the spectral sampling for the higher-and lower-resolution setup are 0.32 Å pixel −1 and 0.64 Å pixel −1 respectively.In terms of instrumental velocity dispersion, the spectral resolution is σ inst = 37 km s −1 at 5300 Å (the location of the prominent Mgb absorption feature) for the higherresolution setup and 48 km s −1 for the lower-resolution setup at 6700 Å around the Hα emission line.

ANALYSIS
The main goal of our spectroscopic analysis is to determine the properties of the stellar populations of the CR disks and ionized gas, important for understanding the evolution of that galaxy.In general, we followed a similar analysis workflow applied earlier to our studies of the CR galaxies IC 719 (Katkov et al. 2013) and NGC 448 (Katkov et al. 2016).

One-component spectral fitting
Age and metallicity determine the appearance of a galaxy's stellar spectrum.The age of the stellar pop-ulation influences the width and strength of absorption lines, with younger populations exhibiting stronger Balmer lines.The metallicity primarily affects the number and strength of metal absorption lines, with higher metallicity populations showing more pronounced metal features.To determine the properties of the stellar population, we applied the full spectral fitting package NBursts (Chilingarian et al. 2007a,b) which models the entire galaxy spectrum and effectively uses the information about the stellar population encoded there.We used NBursts with a grid of simple stellar population models (SSP) built with the evolutionary synthesis code pegase-hr (Le Borgne et al. 2004) based on the R = 10, 000 empirical stellar library Elodie 3.1 (Prugniel et al. 2007).At each model evaluation during the χ 2 -minimization loop, the stellar population spectrum is interpolated from the model grid for a given age T SSP and metallicity [Z/H] SSP values.The spectrum is then broadened with stellar line-of-sight velocity distribution (LOSVD), in our case parameterized by a pure Gaussian, and multiplied by a polynomial continuum to account for the difference between the shape of the model spectrum and the observed one due to imperfect spectral sensitivity calibration and effect of dust extinctions.We used multiplicative Legendre polynomials of the 15th and 20th degrees for the PG2300 and PG1300 setups, respectively.The model also contains a set of strong emission lines (Hγ, Hβ, possessing the same Gaussian kinematics, but decoupled from the stellar kinematics.Emission lines are additive components whose weight, i.e. emission line fluxes, are determined as a linear problem solved in every step of a non-linear minimization loop.The emission line templates as well as a grid of stellar population templates were pre-convolved with wavelength-dependent instrumental resolution before the main minimization loop.This approach of simultaneously fitting the stellar continuum and emission lines with independent kinematics is similar to the gandalf (Sarzi et al. 2017) or recent versions of the ppxf code (Cappellari 2017).
To obtain spatially resolved parameter profiles with reliable measurements, we applied a straightforward adaptive binning scheme to our long-slit data.This involved increasing the bin size to achieve a minimum signal-to-noise ratio of SNR=10 in the continuum at 5100 ± 30 Å in the rest frame.The radial profiles derived by our single component analysis are shown in Figure 2. The kinematic profiles clearly show typical signs of the presence of two CR stellar components (Krajnović et al. 2011), a wavy velocity profile, where in the central region (R < 2 ′′ ) stars rotate in the opposite direction relative to the outer part (R > 5 ′′ ) of the galaxy and a doublepeaked velocity dispersion profile.The two peaks in the velocity dispersion profiles are located at the transition radius R ≈ 3 − 4 ′′ , where the internal rotation, dominated by one component, switches to the external one, which corresponds to a zero velocity amplitude at that distance.The same clear dichotomy is observed in the parameters of the stellar population.The central region is dominated by the younger stellar population, SSP-equivalent age T SSP ≈ 1.5 − 2 Gyr, with metallicity [Z/H] SSP ≈ −0.3 dex.The outer regions are dominated by a much older stellar population T SSP ≈ 7 − 8 Gyr which is metal-poor [Z/H] SSP ≈ −0.8 dex.

Ionized gas
The optical spectra of PGC 66551 exhibit strong emission lines of ionized gas.The NBursts technique used in this study incorporates the contributions of emission lines into the model and provides estimates of line fluxes.Since the spatial distribution of the emission line fluxes and the stellar continuum are not identical, we applied a different binning schema based on the total Hβ+[O iii] flux to analyze the emission lines.This allows us to measure the spatially resolved radial profile of the emission line fluxes and kinematics for the regions extending far beyond the stellar measurements in Section 3.1.
To identify the gas excitation mechanisms we used both the standard Baldwin-Phillips-Terlevich diagrams (BPT, Baldwin et al. 1981) with the demarcation lines by Kewley et al. (2001) andKauffmann et al. (2003) and W Hα versus [N ii]/Hα (WHAN) diagram (Cid Fernandes et al. 2011) for both MaNGA and SALT spectra (see Figure 3).Both methods indicate that the dominant excitation mechanism is ionization by young stars with a slight change toward the composite or LINER regions with increasing distance from the galaxy center.Since no obvious sources of shocks, such as irregular gas motions due to the bar or tidal interaction, are observed, this behavior may be attributed to the increasing relative contribution to ionization by post-Asymptotic Giant Branch (post-AGB) stars of the old stellar population (Yan & Blanton 2012;Singh et al. 2013).
The dominant excitation mechanism by young stars clearly indicates current star formation.We corrected emission line fluxes by using Balmer decrement and the Fitzpatrick (1999) extinction curve.Using the total extinction-corrected Hα flux enclosed by the MaNGA field-of-view covering ≈2R eff and the Kennicutt (1998) relation, we estimated the total star formation rate (SFR) as ≈0.89 M ⊙ yr −1 .This value compares well to 0.93 M ⊙ yr −1 derived from the broadband spectral energy distribution (SED) fitting and provided by the GALEX-SDSS-WISE Legacy Catalog (GSWLC, Salim et al. 2018).
We compared these SFR estimations with the expected values for a star-forming galaxy located on the star formation main sequence (Brinchmann et al. 2004;Salim et al. 2007).The expected SFR value for a starforming galaxy with stellar mass of M * = 2.8 × 10 10 M ⊙ (the mass measurement is described in Section 3.6), is  3 M ⊙ yr −1 using the relation from Salim et al. (2007) or 1.4 M ⊙ yr −1 for time-dependant main-sequence best-fit from Speagle et al. (2014).This comparison implies that the current star formation in the PGC 66551 is slightly quenched and the galaxy has already moved down from the main star formation sequence (if it was there before).
To discuss the origin of the accreted material that formed the counter-rotating disk, we need to know the gas-phase metallicity.To determine the ionized gas metallicity, we utilized the O3N2 and N2 calibrations from Pettini & Pagel (2004).We found that the gas metallicity along the slit and throughout the MaNGA field-of-view shows no variation within the typical uncertainty of 0.1 dex for used calibrations.The average values of 12 + log O/H4 for both the O3N2 and N2 calibrations are 8.54±0.01dex5 and 8.47±0.02dex, respectively.Assuming the solar metallicity to be 8.69 dex (Asplund et al. 2009), our estimates correspond to a slightly subsolar ionized gas-phase metallicity in the range between −0.2 dex and −0.1 dex.We also applied the IZI technique (Blanc et al. 2015) and found that estimates based on all photoionization model grids except d13_kappa20 and d13_kappaINF (Dopita et al. 2013) are in agreement with the O3N2 and N2 calibrations.
We also check whether PGC 66551 falls on the massmetallicity relation for star-forming galaxies.Various metallicity calibrations result in significantly different shapes and positions of the mass-metallicity relation varying by almost 1 dex (Kewley & Ellison 2008).To make a valid comparison, we checked the massmetallicity relation derived for the same calibrations.For a galaxy mass of 2.8×10 10 M ⊙ and the O3N2 and N2 calibrations, we expect a 12 + log O/H value of 8.76 and 8.72 dex, respectively.These values are 0.22-0.25 dex higher than the metallicity in the galaxy and significantly higher than the RMS of the relation itself, which is 0.1 dex (Kewley & Ellison 2008).

Non-parametric stellar LOSVD
To unveil the two-component kinematics of the stellar population in PGC 66551, we applied the nonparametric LOSVD recovery approach, which we have successfully used in our past studies (e.g., Katkov et al. 2011Katkov et al. , 2013Katkov et al. , 2016;;Kasparova et al. 2020).A spectrum of a galaxy logarithmically binned in wavelength can be modelled as a convolution of a high-resolution spectral template of the stellar population with the stellar LOSVD.In our approach, we treat the deconvolution problem as an ill-posed linear problem whose solution is the desired LOSVD.We solve the problem using a linear least-squares method with a regularization term that requires smooth solution of the LOSVD and minimiza-tion of artifact peaks away from the systematic velocity (so-called "tail" regularization).As a template, we used the unbroadened bestfit model from the NBursts one-component fitting (Section 3.1).Here we also utilized the latest modification of our method that implements an iterative approach (Gasymov & Katkov 2021).Rather than finding the full LOSVD solution, we reformulate the problem to search for corrections required on top of the initial approach of the LOSVD profile.For the first iteration, we use the Gaussian LOSVD shape from the one-component fitting.For subsequent iterations, we use the full solution from the previous step.In this approach, the regularization operates on the correction vector rather than the full LOSVD solution, requiring the corrections to be smoothed.
The resulting stellar LOSVD profiles reconstructed in all spatial bins along the slit are shown in the top panel in Figure 4 as a position-velocity diagram and clearly demonstrate the existence of two kinematically distinct components in the stellar population -the inner one dominates in the central part and co-rotates with the ionized gas, while the main stellar component counterrotates to the gas and extends to the galaxy outskirts.We have made the code of our method publicly available as a Python package sla6 .Note that a similar but Bayesian-based approach for non-parametric stellar LOSVD recovery bayes-losvd by Falcón-Barroso & Martig ( 2021) is also available.
To estimate the contribution of the secondary CR stellar component to the integral stellar light, we fit the reconstructed LOSVD.Assuming that at any given position of R, the LOSVD is parameterized by two Gaussians those parameters sigma σ, flux F and velocity v are functions of R. For flux and sigma we used exponential decaying laws F (R) = F 0 exp(−|R|/h f ), σ(R) = max{σ 0 , σ 1 + k|R|}, while for velocity profile we adopted following parameterization: v(R) = v max (tanh(π|R|/R 0 ) + c|R|/R 0 ).The best-fitting LOSVD model and its residuals are shown in the central and bottom panels in Figure 4.This model suggests that a secondary stellar counter-rotating disk contributes (42 ± 2)% to the integral light of PGC 66551 galaxy.
One limitation of our LOSVD reconstruction approach is that we use only one stellar template to reconstruct the complex shape of the LOSVD formed by two stellar disks.It is assumed that both stellar components have comparable stellar population properties and are reliably represented by the same stellar template.If this assumption is incorrect, this could lead to a distortion of the recovered LOSVD.To overcome the aforementioned problem of using an identical stellar population template and obtain the stellar population properties for both kinematically distinct components, we performed a two-component parametric fit, also known as a spectral decomposition (Coccato et al. 2011).In this method, the observed spectrum is modeled by two stellar populations, broadened by individual and independent LOSVDs.We previously implemented this approach based on the full pixel fitting technique NBursts when studying the stellar counterrotation phenomena in IC 719 (Katkov et al. 2013), NGC 524 (Katkov et al. 2011), NGC 448 (Katkov et al. 2016).In contrast to those studies, where we masked emission lines, here we incorporated a set of emission line templates into the model.A decomposition of the spectra from the region R ≈ 2.5 ′′ is shown in Figure 5.

Two-component spectral decomposition
The derived radial profiles of the recovered parameters are shown in Figure 6.The profiles are less extended compared to the single component analysis because, in this case, we applied a more stringent SNR requirement of 15 during the adaptive binning stage.They show obvious differences in the line-of-sight velocity, age, and metallicity.The main and counter-rotating stellar disks, as well as the ionized gas component, exhibit a similar rotation amplitude of 100 − 120 km s −1 (not corrected for disk inclination).However, the main disk counterrotates with respect to the gas, while the secondary component co-rotates.The results of the spectral decomposition clearly demonstrate that the counter-rotating stellar disk dominates the integral light in the central 3 ′′ (see right column in Figure 6).Our assumption that both components have disky morphology is also supported by the radial profile of stellar kinematics.Both disks have comparable contributions to the integral light at radii R = ±3 ′′ where two σ-peaks have been detected in the single component fitting profiles (see Figure 2).The main disk is significantly older T SSP = 5 − 10 Gyr compared to the counter-rotating disk of ≈ 1 Gyr, indicating that its stars formed much later.Finally, the metallicity distribution provides a clear distinction between the two disks.In contrast to most galaxies studied in detail using the spectral decomposition method (see Figure 7), the stellar population of the main stellar disk of PGC 66551 exhibits a significantly lower metallicity of [Z/H] SSP ≈ −0.8 dex in comparison to the counterrotating disk which has slightly sub-solar metallicity.Using the relative weights obtained from the spectral decomposition, we estimated that the contribution of the CR component to the integral light was approximately 50%, and its contribution to the stellar mass was estimated to be around 50%(M/L) CR /(M/L) main ≈ 20%, based on the mass-to-light ratios derived from the stellar population parameters (see Table 1).

Spectrophotometric fitting
In this Section, we propose a more complex and advanced approach to investigate PGC 66551 by simultaneous modeling of both the spectra and broad-band photometry aimed at overcomming the age-metallicity degeneracy and confirming the low metallicity of the main disk, as well as producing a self-consistent stellar population model needed for the analysis of the chemical evolution in Section 4.
An increase in the age and metallicity of the stellar population results in the reddening of SED and the effect is particularly strong in the ultraviolet and blue filters (Conroy 2013).We have compiled the spatially resolved SED including ultraviolet FUV, NUV bands from Galaxy Evolution Explorer (GALEX) obtained from the MAST archive7 , optical griz images from SDSS survey8 (u-band image was avoided since it is typically faint), grz images from Legacy Survey (Dey et al. 2019) Bin R = 1.5 − 2.8 arcsec tracted from cutout service 9 , i-band CFHT image extracted from CADC 10 , JHKs images from the The VISTA Hemisphere Survey (VHS) retrieved through the ESO Archive Science Portal 11 .We estimated and subtracted sky background level applying σ-clipping tech-9 https://www.legacysurvey.org/viewer 10 https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca 11 http://archive.eso.org/scienceportal/nique12 .We used azimuthally averaged surface brightness profiles constructed by using isophote module from Astropy-affiliated photutils package.
We used the spec+phot mode of NBursts (Chilingarian & Katkov 2012; Grishin et al. 2021;Katkov et al. 2022) to simultaneously model the optical spectra and broadband SED in two spatial bins R = ±2.5 ± 0.3 ′′ where the contributions of both disks are comparable.
We initially attempted to use a combination of two SSP models, as in Section 3.4.However, due to the complex star formation history (SFH) and the current star formation in PGC 66551 galaxy, no combination of two SSP models was able to describe the ultraviolet part of the broadband SED.Then we applied stellar population models with an exponentially decreasing star formation history (hereafter exp-SFH models) defined by the initial epoch of star formation T start , the exponential decay scale τ , and the mean stellar metallicity.We used and discussed the same set of stellar population models to analyze stellar population properties of galaxies in the RCSED catalog (Chilingarian et al. 2017).
While extinction A V is a crucial parameter for spec+phot modeling, it is an external parameter in the NBursts package, which means that it is not fitted, but fixed at some predetermined value.We performed the fitting for a range of A V values from 0.0 to 1.5 with a step size of 0.05 and found that the best χ 2 values for both spatial bins corresponded to A V = 0.65 mag.Thanks to strong emission lines in the analyzed spatial bins we estimated extinction as A V = 0.63 ± 0.24 mag from the Balmer decrement Hα/Hβ which is in remarkable agreement with that derived from the spectrophotometric fitting.
Another constraint on the model parameters we could make from the current SFR estimate, which can be recalculated into τ CR value for any given formation epoch of a CR disk T start,CR .The main version of NBursts can only minimize two-parameter grids of stellar population models.Hence, we fixed the initial epoch for the main disk as T start,main = 13 Gyr and performed a grid search along all nodes for the epoch of the CR component.The analysis of both spatial bins revealed that the optimal value for the formation epoch of the CR disk is T start,CR = 2 Gyr.The completely free parameters remaining in the fitting procedure were the metallicities of both components, the exponential decay scale of the main disk τ main , and the kinematics parameters, which we do not discuss here, but are generally consistent with the SSP analysis.The final parameters of the fitting in two spatial bins are given in Table 1. Figure 8 shows the bestfit model and model components for the bin R = 2.5 ′′ .
The proposed much more complex and comprehensive analysis preserves our conclusion that the main component of PGC 66551 does have a metall poor stellar population.Compared to the SSP analysis ([Z/H] SSP ≈ −0.8 dex), the metallicity [Z/H] expSFH ≈ −0.5 dex is not as extremely low, but remains way lower than the CR disk metallicity ≈ −0.1 dex in PGC 66551 and strongly lower than that in the main components of the other studied CR galaxies (see Figure 7).
Also in Figure 9 we illustrate how spectrophotometric fitting breaks the age-metallicity degeneracy.We draw slice of the χ 2 space on the plane of the main stellar disk parameters [τ main , [Z/H] expSFH ] for several fitting variants: i) where only spectral information is analyzed (α = 0.01), ii) our main model where spectral and photometrical data are used (α = 0.5) and iii) purely photometrical data fit (α = 1).The less elongated and more circular shape of 1σ-and 3σ-contours around minimum in the α = 0.5 map clearly suggests a lower degree of degeneracy between the age-related parameter τ main and metallicity.The SALT long-slit spectrum along the major axis shows that the emission lines trace the galaxy rotation to 20 ′′ ≈ 8 kpc, where it has already reached a plateau.We used this extended rotation curve to estimate the total dynamical mass of the galaxy and the fraction of the dark matter halo.First, we corrected the determined line-of-sight velocity measurement for the inclination effect by adopting the angle i = 30 • from the isophote analysis of the r-band SDSS image carried out by isophote module from Astropy-affiliated photutils package13 (Bradley et al. 2022).Since the decomposition of the rotation curve is ambiguous (see e.g.Saburova et al. 2016), which leads to a strong degeneracy between parameters, we minimized the number of free model parameters to stabilize the solution.We decomposed the rotation curve into two stellar disks (main and counter-rotating) and a dark matter halo.

Rotation curve decomposition
To make this analysis consistent with the stellar population analysis, we applied relative weights derived in the spectral decomposition (see Section 3.4) to the rband radial surface brightness profile to calculate individual profiles for both disks.The profiles were found to be similar to exponential profiles, and therefore, we fit them by exponents and used their parameters to simplify the calculation of their contribution to the rotation curve.To convert the light profiles into mass profiles, we explored the refined mass-to-luminosity (M/L) ratios for adopted stellar population parameters τ RC , [Z/H] RC from the spectrophotometric analysis (see parameters in Table 1).This also yielded in a total mass of 2.8 × 10 10 M ⊙ consisting of M main = 2.3 × 10 10 M ⊙ and M CR = 4.8 × 10 9 M ⊙ .In this way, the contribution from the stellar disks was totally fixed.We tried to add a bulge component to the model, but the central part of the rotation curve was already well described and there was no empty room for an additional component.The parameters of the dark matter halo remained free.Simi- lar to Saburova et al. (2016) we considered three types of DM density profiles: Navarro-Frenk-White (NFW) profile (Navarro et al. 1996), the density profile by Burkert (1995) and the pseudo-isothermal profile (hereafter, piso).All three profiles are parameterized by central density ρ 0 and radial scale of the halo R s .For convenience, we re-parameterized the profiles in parameters R s and V max , where V max is a peak halo velocity at radius R max .The example of the χ 2 -map for the NFW halo in the (V max , R s ) coordinates is shown in Figure 10 obviously demonstrates degeneracy even between V max and R s .In the last simplification -we constrained V max to be equal to the value of the rotation curve plateau of 250 km s −1 with only one remaining model parameter R s .Different halo profiles provided similar mass estimates: M NFW = (1.6 ± 0.4) × 10 12 M ⊙ , M burk = (0.8 ± 0.1) × 10 12 M ⊙ , M piso = (3.0 ± 0.6) × 10 12 M ⊙ .We found that these halo mass estimates are in perfect agreement with empirical calibration M halo − V flat (Marasco et al. 2021) built on the SPARC galaxy sample (Lelli et al. 2016) and rotation curve decomposition by Posti et al. (2019).This relation predicts M halo = 1.9 × 10 12 M ⊙ for velocity plateau of 250 km s −1 .Combining M NFW and stellar mass (from the GSWLC and our estimates), we obtained the stellar fraction f ⋆ = M ⋆ /f b M halo , where f b is a baryon fraction 14 , in the range of 0.04−0.06.In general, this value is in the cloud of the SPARC galaxies on the f ⋆ − M ⋆ plot from Posti et al. (2019) (Figure 2) and noticeably lower than the prediction from abundance matching study by Moster et al. (2013) (grey stripe on the same Figure 2).This comparison for PGC 66551 implies that, despite its kinematic peculiarities and very unusual way of building its stellar disk, the ratio between the stellar components and dark matter halo is not significantly abnormal.In this section, we construct a framework for modeling the history of the chemical enrichment, which can explain the properties of the stellar population of CR disks and ionized gas in PGC 66551.Since we found no significant gradients in the stellar population properties and gas metallicity, we can consider the chemical enrichment as a global process, assuming the entire galaxy as a single well-mixed zone with a uniform star formation history and metallicity.We exploited instantaneous recycling approximation (IRA, Schmidt 1963;Pagel & Patchett 1975;Tinsley 1980), which assumes that stars with m < 1 M ⊙ have infinite lifetime, while massive stars m ≥ 1 M ⊙ instantaneously die as they form enriching the ISM by newly produced metals, which immediately mix into the gas phase and are instantly recycled into new generations of stars.We constructed the IRA model of chemical evolution following a detailed description by Spitoni et al. (2017) and Chilingarian & Asa'd (2018).The model typically includes the following ingredients: SFH ψ(t), initial mass function (IMF), the total mass of the gaseous reservoir, and infall I(t) and outflow O(t) rates.For the convenience of formulating equations, the following parameters are introduced: R -the so-called returned mass fraction of gas that returns to the ISM and can be recycled for star formation and y Z -the yield per stellar generation that is a fraction of mass in heavy elements synthesized by stars and returned to the ISM through stellar winds and supernova explosions.Both R and y Z depend on the IMF and much less on metallicity (Vincenzo et al. 2016).We adopted the following values from Spitoni et al. ( 2017) R = 0.441 and y Z = 0.0631 for the Kroupa (2001) and Chabrier (2003) IMFs.
To calculate the evolution of metallicity in the ISM Z = M Z /M gas we can define a set of equations.In the closed-box model, the total mass of metals M Z in the ISM decreases as some metals are locked up into stars Ṁ − Z = Z Ṁgas = −Z Ṁstars , while it increases due to enrichment from the evolution of massive stars Ṁ + Z = y Z Ṁstar .The balance of mass of heavy elements can be expressed as ṀZ = Ṁ − Z + Ṁ + Z = (−Z + y z ) Ṁstar .We defined a set of equations for the evolution of the galaxy mass M tot , gas mass M gas and metallicity, taking into account terms for infall I(t), outflow O(t) and Ṁstar = (1 − R)ψ(t), as follows: (1) where Z inf and Z out are metallicities of the infalling and outflowing gas, respectively.Combining the above equations we have that We used a simple model where the outflow rate is proportional to the SFR O(t) = λψ(t) (Matteucci & Chiosi 1983;Matteucci 2012).The mass-loading factor λ represents the efficiency with which the gas is expelled from the galaxy by galactic winds.Following Dalcanton (2007) we consider the metallicity of the outflowing gas as follows: Z out = ϵZ + (1 − ϵ)Z SN .The outflow material comes completely from the ISM in the case of the entrainment fraction ϵ = 1.ϵ = 0 means that the outflow is dominated by supernova ejecta.The metallicity of the supernova ejecta is Z SN = ηy Z , where η = 4.5 − 5.2 for the Kroupa ( 2001) IMF (Dalcanton 2007).For our calculations, we adopted η = 4.85.For the infall rate we used the most common exponential formulation I(t) = A exp (−t/τ inf ) (Matteucci 2012).We applied a fourth-order Runge-Kutta method for Equation 2 with initial conditions [Z/H](0)=-3 (Bromm & Yoshida 2011), M star (0) = 0, M gas (0) = M tot (0) = M 0 .Despite the simplicity of such an approach to chemical evolution, it allows one to obtain the physical properties of galaxies at different stages of their lifetime and to model and retrieve from spectra their star formation histories (Grishin et al. 2019), as well as additional information on infall and galactic winds (Chilingarian & Asa'd 2018;Grishin et al. 2021).

Application for PGC 66551
Now that we have defined our chemical enrichment framework, we can adapt it to the context of the studied galaxy.We are aiming to constrain the properties of the outflowing galactic wind and the metallicity of the infalling gas by reproducing the observables in PGC 66551.We consider two options for constructing the model: 1) first reproducing the properties of the main disk, and then the CR disk, and 2) first reproducing the CR disk, and then the main disk.
Using the exp-SFH scales τ main , τ CR from the spectrophotometric analysis (see Section 3.5) and stellar mass estimates of the disks we can fully constrain their SFHs ψ main and ψ CR .We assumed that the main stellar disk was formed from a primary pristine baryonic reservoir of mass M 0 .The formation began 13 Gyr ago and proceeded according to exp-SFH with τ main = 4.8 Gyr and resulted in the disk of stellar mass M main = 2.33 × 10 10 M ⊙ .
The fact that the externally accreting material formed a large-scale CR disk means that the host galaxy did not contain a significant amount of gas at the epoch T start,CR = 2 Gyr ago when the accretion started.Otherwise, the gas clouds would collide, lose angular momentum, and fall toward the galaxy center without forming a CR disk.
Formation of the CR disk started T start,CR = 2 Gyr ago with the accretion of gas with an exponential rate I(t) ∝ exp (−t/τ inf ) and metallicity Z inf .From spectrophotometric analysis, an exponential SFH with a scale τ CR = 0.98 Gyr is required for the current star formation to be at the estimated level of 0.93 M ⊙ yr −1 .Initially (at the moment T start,CR ), the CR SFR is highest ≈10 M ⊙ yr −1 and no gas, therefore, in order to have a physically meaningful situation, it is necessary to ensure a sufficient gas content.Therefore, we delayed star formation from starting moment until 50 Myr after the gas infall began.The change in this time offset does not have a significant effect on the results of modeling.In our baseline model, we adopted that the exponential scale of infall τ inf is equal to exp-SFH scale τ CR .Later we discuss variations of this choice.Fixed τ inf allowed us to calculate the normalization of the infall exponential parameterization using the condition that the residual mass of the gas is equal to the mass of atomic hydrogen M HI = 1.35 × 10 9 M ⊙ from the Hi-MaNGA survey (Masters et al. 2019).Here, we assumed that externally accreted stars comprise an insignificant fraction.
The generation of stars formed at any given time inherits the ISM metallicity Z(t).For any given CEH model of the ISM, a combination of ψ(t) and the massto-light ratios from the grid of stellar population models can be used to calculate the luminosity-weighted average metallicity of each disk component.Requiring these values to be equal to those found from the spectrophotometric fit, [Z/H] main = −0.55dex and [Z/H] CR = −0.2dex, we can obtain constraints on the parameters of the galactic wind and the metallicity of the infalling gas Z inf formed the CR disk.For a given λ, to reproduce the low metallicity of the main component, one might want to have a metal-rich galactic wind (low ϵ), but a higher metallicity of the CR disk requires preserving rather than expelling metals (higher ϵ).Therefore, the properties of both disks set the limits on the entrainment fraction ϵ of the galactic wind for a given mass-loading factor λ.
We assumed that the outflow process is identical in terms of λ and ϵ during the formation of both disks.We repeated the below calculations for every λ value from the set of 0.1, 0.15, 0.2, 0.25, 0.3, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5.First, we focused on the properties of the main disk.ϵ is only one free parameter that we can vary to reproduce [Z/H] main = −0.55dex.Low metallicity of the main disk controls the lower boundary in the ϵ − λ diagram shown in blue solid line in Figure 11 (top left panel).Then, determined ϵ was used for the history of the enrichment of the CR disk, where we varied the infalling gas metallicity Z inf to reproduce [Z/H] CR = −0.2dex (solid blue line in the bottom left panel in Figure 11).
Then we built a model of the metal enrichment focusing on the properties of the CR disk.In the first step, we considered a variant of the accretion primordial pristine gas assuming Z inf = 10 −3 Z ⊙ (Bromm & Yoshida 2011).Reproducing the stellar metallicity of the CR disk, we obtain ϵ, which is the maximum possible value.In all other scenarios of pre-enriched gas infall, the ϵ will be lower, i.e. more metal-rich wind is needed to get the observed stellar metallicity.This defines the upper boundary shown by the dashed black line in the ϵ − λ diagram.In this inefficient wind case (higher ϵ), in order to derive the low metallicity of the main disk, we have to increase the initial reservoir of metal-poor gas.This results in residual gas at the look-back time T start,CR = 2 Gyr which needs to be expelled before the massive accretion of gas forming the CR disk.Exploring Illustris simulations Starkenburg et al. (2019), identified two channels responsible for a significant episode of gas removal: an intense burst of AGN feedback and gas stripping occuring during a flyby passage through a massive group environment.Another possibility, as described by Khoperskov et al. (2021) and also based on cosmological simulations (IllustrisTNG), suggests gas removal through the capturing and mixing of pre-existing gas with infall.However, this scenario implies that the mass of infalling gas substantially exceeds the pre-existing one.Below, we will attribute the episode of gas loss to AGN, implying that this could be any physical process resulting in the effective expelling of residual gas by the time the CR disk formation occurs.
The oxygen abundance in the ionized gas [O/H] = −0.15dex is another observable quantity available to us.To use it, we added the oxygen enrichment history calculation to the model by solving an identical Equation 2 assuming the oxygen yield y O = 0.0407 (Spitoni et al. 2017) instead of y Z .We found that in the latter scenario, the oxygen abundance appeared to be much higher than the value measured in the galaxy15 (see grey dashed line the bottom right panel in Figure 11).To reconcile this discrepancy, the winds must blow metals out more efficiently (lower ϵ).We minimally reduced ϵ requiring that the model-predicted oxygen abundance be equal to the observed value, this corresponds to the upper boundary shown by the solid red line in the ϵ − λ diagram.The metallicity of the CR disk is adjusted according to a variation in Z inf (red line in bottom left plot in Figure 11), while the main disk metallicity is fitted by variations in M 0 .In this scenario, the mass of residual gas at T start,CR = 2 Gyr is significantly lower compared to the case of primordial accretion with Z inf = 0 without reproducing oxygen [O/H] abundance in the ionized gas.Evolution of the total gas mass and metal enrichment history for different λ and scenarios are shown in the middle and bottom right panels in Figure 11.
We chose a very simple exponential infalling rate with τ inf = τ CR as a baseline value.However, there is no physical reason for this, therefore we tested a set of different infall rates.We found that the range of possible τ inf values is physically limited.Indeed, a very high gas infall rate (small τ inf ) would result in a rapid accumulation of gas in the initial stage and a slow enrichment track, whose averaged value is final metallicity at T = 0 can be lower than the measured ones.In the case of prolonged accretion (large τ inf ), the incoming gas can be totally exhausted and SFR cannot be maintained at the predicted level.Depending on λ both types of problems occur at different τ inf .We chose a secure range of τ inf between 0.6τ CR and 1.3τ CR for all considered λ up to 2.5.For the lower λ = 1 the range is slightly broader; τ inf can be between 0.5τ CR and 1.5τ CR .We found that different infall regimes result in different shapes of enrichment tracks.Namely, the short-scale infall corresponds to gradual enrichment, while the long-scale one to rapid metall enrichment.In the scenario without expelling gas by AGN at look-back time T start,CR , the luminosityweighted metallicity changes insignificantly resulting in very little variation in [Z/H] inf = 0.45 dex.In the scenario with AGN feedback, [Z/H] inf is calculated based on the current ISM metallicity match, which is the last point on the track and is much more sensitive to the shape of the Z(t).In this scenario, the [Z/H] inf values may vary between −0.9 dex for τ inf = 1.5τCR and −0.5 dex for τ inf = 0.5τ CR .Given the difficulties in measuring absolute metallicity values in the ISM (Kewley & Ellison 2008), we tend to consider [Z/H] inf ≈ −0.5 dex as more reliable estimate for infall metallicity, keeping in mind that possibly it can be as low as −0.9 dex.
Above, we adopted the atomic hydrogen mass M HI = 1.35 × 10 9 M ⊙ as a proxy for the total residual gas mass at T = 0 Gyr.However, the unknown contribution of molecular gas can significantly increase this value.Considering that M H2 /M HI can vary widely (Boselli et al. 2002;Keres et al. 2003) we also re-computed our modeling for an extreme case where the residual mass is equal to M gas,now = 10 10 M ⊙ corresponding to M H2 /M HI = 7.4.We observed that this modification implies a more massive gas infall and, consequently, a slower enrichment process in the CR disk.Therefore, to achieve the observed metallicity [Z/H] CR , one needs to start with a higher metallicity of the accreting material, specifically [Z/H] inf = −0.3dex in the no-AGN scenario and −0.7 dex in the AGN scenario.Additionally, the upper bound of entrainment fractions (ϵ) has increased, making this parameter less constrained.In addition, we explored a scenario in which residual gas from the main disk formation was expelled (by an AGN event for instance) before the beginning of the CR disk formation (T start,CR ).For example, the gas expelling at epoch T = 4 Gyr ago resulted in almost unnoticeable changes in wind parameters, with only slight variations in [Z/H] inf decreasing by 0.05 dex.Both experiments indicate the stability of our simple modeling to various modifications, with one of the crucial output parameters for interpretation, the metallicity of the infall gas [Z/H] inf , showing only marginal variations.

DISCUSSION
Comparing PGC 66551 with other CR galaxies studied in detail using the spectral decomposition method, the most notable features are the ultra-low metallicity of the main disk and the exceptionally large difference in the ages of the stellar disk populations (see Figure 7).The first emerged idea was that before the galaxy accreted material to form the CR disk, something stopped the formation of the main disk and kept the disk chemically undeveloped.Indeed, the exceptionally intensive episode of AGN activity might cause such a dramatic effect (Silk & Rees 1998;Springel et al. 2005).Also, AGN-driven winds inevitably interact with inflows and can drastically reduce accretion rates, affecting star formation in the galaxy (van de Voort et al. 2011;Nelson et al. 2015).The AGN feedback was explored as one of the mechanisms for gas removal, allowing the subsequent formation of the kinematically misaligned components (Starkenburg et al. 2019;Duckworth et al. 2020).
The majority of detailed studies of counter-rotating galaxies have addressed the question of formation sce-nario and source of the externally accreted material that formed the CR disks.Two main possibilities are usually considered -accretion from cosmological filaments (Algorry et al. 2014;Khoperskov et al. 2021) and merger with gas-rich satellites (Thakar & Ryden 1996, 1998;Lu et al. 2021).In our study, we complemented the spectral decomposition approach (and its spectrophotometric extension) by simple chemical evolution modeling which gives some results relevant to this discussion.
In the configuration with AGN feedback episode, our model suggests the metallicity of the infalling gas [Z/H] inf to be in the range from −0.9 dex to −0.5 dex depending on the adopted model of the gas accretion.The scenario without AGN feedback gives evidence for −0.45 dex and is dependent on neither the infall rates nor the wind parameters.In both scenarios, the expected infall metallicity is higher than that of the nearpristine gas, which we would naively expect in case of accretion from the cosmological filaments.Direct observational evidence of gas accretion is still extremely challenging and remains elusive (Fox & Davé 2017).However indirect measurements of gas reservoirs in the circumgalactic medium (CGM), probed by absorption features from background quasars, suggest a certain degree of metal enrichment in the CGM, ranging from 2.5% to 10% of solar metallicity Prochaska et al. (2013); Lehner et al. (2013); Prochaska et al. (2014).Numerical simulations also suggest that the gas accretion should be metalpoor, but not pristine because of the interaction of cold accretion flows with already enriched galactic halos and outflowing material (van de Voort & Schaye 2012; Hafen et al. 2017).We have a preference towards [Z/H] inf = −0.45dex because this value is based on a comparison of the stellar metallicity, which is an average over the entire enrichment history.[Z/H] inf = [−0.9,−0.5] dex seems less reliable, as it depends on the accretion mode and is based on a comparison of the current metallicity in the ISM, which measurements strongly depend on the applied methods (Kewley & Ellison 2008).
We also considered gas depletion time t dep = M gas /SFR suggested by different infall scenarios.More prolonged accretion caused a gradual increase of gas content in the initial stage of CR disk formation which corresponds to a fairly short depletion time 30 − 130 Myr depending on the wind parameter λ.Typical depletion time in spiral galaxies is 1 − 2 Gyr (Bigiel et al. 2011).
More extreme values can be achieved locally in the star forming regions down to 50 − 100 Myr (e.g., Evans et al. 2009Evans et al. , 2014;;Lada et al. 2012).More rapid infall leads to faster gas accumulation and as a result, to higher depletion times by factor of two.Despite still being quite short, the tendency to longer and more realistic depletion time gives support to the rapid rather than prolonged infall.In symmary, both infall metallicity and regime arguments seem to be more naturally interpreted as a merger with a gas-rich satellite galaxy with the subsequent in situ chemical enrichment, rather than a prolonged gas accretion from a cosmic filament.
6. SUMMARY In this study, we used deep long-slit spectra obtained with the 10 m class SALT telescope to investigate the galaxy PGC 66551 hosting two large-scale counterrotating stellar disks identified in the SDSS MaNGA survey.Despite the kinematic peculiarity of PGC 66551 which directly points to a very unusual process of building stellar disk, the galaxy generally follows the main scaling relations for disk galaxies.We found that (i) the metallicity of ionized gas is only slightly lower than expected from the mass-metallicity relation, (ii) the current star formation rate is also slightly suppressed compared to that for star forming galaxies of a given mass residing on the 'main sequence of star formation', and (iii) the baryon-to-total mass ratio is comparable to other galaxies but slightly lower than the prediction from simulations.We employed a novel spectral decomposition technique and its extension for the simultaneous fitting of spectra and broad-band photometry to determine stellar population parameters for both the main and counter-rotating disks.We demonstrated that simultaneous fitting of spectra and SED can effectively mitigate the age-metallicity degeneracy.
We found a remarkable feature of PGC 66551 among the counter-rotating family.While the young CR disk with sub-solar metallicity contributing ≈20% to the total stellar mass is quite typical for CR galaxies, the main disk populated by old and significantly metal-poor stars has not previously been detected.
We have developed a simple framework for modeling the history of the metal enrichment, which has revealed interesting results for the interpretation of PGC 66551 formation.The metallicities of both disks and current oxygen abundance in the ISM can be reproduced in scenarios with and without AGN feedback.The assumption of the same properties of the galactic wind during the formation of the main disk and the CR disk makes it possible to constrain wind parameters.This model also indicates a relatively short time scale of the gas infall, the metal content in it appears to be significantly pre-enriched which is difficult to interpret as evidence for accretion from cosmological filaments.Therefore, we concluded that the formation of stellar and ionized gas counter-rotation in PGC 66551 galaxy is the result of a merger with a gas-rich satellite.

ACKNOWLEDGMENTS
We are grateful to the anonymous referee for useful comments and suggestions.We also thank Anatoly Zasov, Mark Krumholz and Andrea Macciò for their valuable comments and discussion on this study.Authors acknowledge the support from the Russian Scientific Foundation grant No. 21-72-00036 and the Interdisciplinary Scientific and Educational School of Moscow University "Fundamental and Applied Space Research".ER also acknowledges the RScF grant No. 23-12-00146 for supporting the optimization of the NBursts method for this study.Spectral observations reported in this paper were obtained with the Southern African Large Telescope (program 2020-1-SCI-002) supported by the National Research Foundation (NRF) of South Africa.AYK acknowledge the Ministry of Science and Higher Education of the Russian Federation grant 075-15-2022-262 (13.MNPMU.21.0003)This work made use of Astropy:16 a communitydeveloped core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration et al. 2013, 2018, 2022).
This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.This research has made use of NASA's Astrophysics Data System Bibliographic Services.

Figure 1 .
Figure 1.A composite RGB image of the PGC 066551 galaxy taken from the DESI Legacy Imaging Surveys (Dey et al. (2019), legacysurvey.org) is shown in (a).SALT long slit orientation is shown by dashed yellow line while the hexagon presents MaNGA field of view.The next three panels show stellar velocity (b), stellar velocity dispersion (c) and ionized gas velocity based on Hα line (d), all taken from MaNGA DR16 Data Analysis Products (DAP).stellar rotation direction reverses as one moves from the central part to the outer region of the disk.Two off-center peaks in the velocity dispersion map are another evident feature of the stellar CR disk (Krajnović et al. 2011) (Figure 1c).The kinematic position angle, estimated based on the Hα velocity map using the XookSuut code (López-Cobá et al. 2021), yields P A kin = 102.5 • ± 0.2 • .Meanwhile, the photometrical positional angle derived from the r-band SDSS image using the isophote module from photutils package 3 (Bradley et al. 2022) is P A phot = 105• ± 2 • .Taking into account the limited MaNGA field-of-view, these estimates appear to be consistent within ≈ 1σ.The similarity of the P A values supports the idea that the gas is settled and rotates in nearly the same plane as the stellar disk.However, it is not possible to directly compare the rotation of the gaseous disk with that of the stellar one, since the former is a combination of two stellar disks.

Figure 2 .Figure 3 .
Figure 2. Result of one-component analysis of the long-slit spectrum.From left to right panels show line-of-sight velocities, velocity dispersions, SSP-equivalent ages and metallicities.Brown symbols correspond to kinematics of ionized gas.

Figure 4 .
Figure 4. Recovered non-parametric stellar LOSVD (top panel) based on the SALT long-slit spectrum.The recovered LOSVD is normalized to unity in each spatial bin.The two counter-rotating component model and its residuals are shown in the middle and bottom panels.Two counterrotating components are clearly detected, and the CR disk is clearly visible within ≈ 5 ′′ .

Figure 5 .Figure 6 .
Figure 5. Example of the two-component spectral decomposition in the spatial region R ≈ 2.5 ′′ .This figure shows only a fragment of the analyzed SALT spectrum (black line) to demonstrate the fine structure of the absorption lines, in particular Mgb.Stellar best-fit spectral model shown by red line.Orange lines highlight the emission line model.Main and counterrotating stellar components are shown in blue and green colors, correspondingly.Residuals are presented by small dark grey dots bordered by the error level (black lines).

Figure 8 .
Figure 8. Spectrophotometric fitting at spatial bin R = 2.5 ± 0.3 ′′ .Top and bottom right panels are identical to that in Figure 5, while the left bottom panel shows broad-band photometrical SED fitting.Black circles are photometrical data for a given spatial bin.Red, blue, and green circles are bestbit model, main and CR components.

Figure 9 .
Figure9.χ 2 scan in the model parameter space for spectrophotometric analysis at spatial bin R = 2.5 ± 0.3 ′′ .Maps show χ 2 values in the plane of main stellar disk parameters for a set of different fitting modes: purely spectral fit (α = 0.01), spectrophotometric fit (α = 0.5) and purely photometrical fit (α = 1).Less elongated and more circular contours on the α = 0.5 map suggest lower degeneracy between parameters.The red circle shows bestfit solution from fitting procedure described in Section 3.5.

Figure 10 .
Figure 10.Rotation curve decomposition.The top panel displays the rotation curve (magenta stars) corrected for inclination along with the best-fit model (dashed red), disk components (green and blue lines), dark matter halo (cyan line), and residuals (gray dots and line).The bottom panel shows the χ 2 -map in the halo parameters coordinates of Vmax and Rs.The red line indicates Vmax = 250 km s −1 , which is adopted from the rotation curve.The red star represents the position of the best-fit value.the black solid lines depict isolines of the χ 2 values, and the black dashed lines demonstrate equal mass of halo.

Figure 11 .
Figure 11.Modeling of metal enrichment history.Top left panel shows ϵ − λ diagram.ϵ is entrainment fraction parameter controlling the metal content of the outflowing galactic wind.Red and blue colors correspond to scenarios with and without episode of AGN activity expelled remaining gas at look-back time T = 2 Gyr.Bottom left panel presents metal abundance of the infalling gas formed CR disk.Top right panel demonstrates the evolution of total gas mass for different scenarios and mass-loading factors λ.Tracks for λ = 0.5 are shown by thick lines, thin lines show other λ values.Tracks before T = 2 Gyr corresponds to evolution of the main disk, then the evolution of the CR disk.Bottom middle and right plots shows metal enrichment tracks for main and CR disks.Horizontal yellow and purple lines represent metallicities of main and CR disks, which are reproduced in our modeling as luminosity-weighted average values derived from the Z(t) tracks and mass-to-light ratios.
. This material is based upon work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CASS.Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah.The SDSS web site is www.sdss.org.SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.