The Physical Properties of S0 Galaxy PGC 26218: The Origin of Starburst and Star Formation

We obtained optical IFU spectroscopic observations of PGC 26218 from CAHA on 2016 March 11 and 12. The 3.5 m telescope covers two optical overlapping setups. The red spectrograph covers the wavelength range from 3745 to 7500 Å with a low spectral resolution (R ∼ 850, V500), while the blue covers the wavelength range from 3400 to 4840 Å with a medium spectral resolution (R ∼ 1650, V1200). In order to obtain a filling factor of 100%, a three-pointing dithering scheme was used. The exposure time per pointing was 900 s for the red (taken on the first day) and 1800 s for the V1200 (split in two individual exposures of 900 s). We used a Python-based pipeline for reduction of the PPAK data based on an upgraded version of García-Benito et al. (2015) and Sánchez et al. (2016). The reduction process can be summarized in the following steps: identification of the position of the spectra on the detector along the dispersion axis; extraction of each individual spectrum; distortion correction of the extracted spectra; wavelength calibration; fiber-to-fiber transmission correction; flux-calibration; sky-subtraction; cube reconstruction; and finally differential atmospheric correction. To reduce the effects of vignetting on blue part (up to ∼4200 Å) of the V500 data in some parts of the field of view (FOV); (see Figure 11 of Husemann et al. 2013), we combine both setups producing a so-called COMBO data cube. The V1200 spectral resolution is degraded to match the V500 data. We then combine the data from both data sets averaging the V1200 cube and V500 cube in the overlapping wavelength range, weighted by the inverse of the error. For the remaining wavelength range, the blue part corresponds to the matched V1200 cube and the red part to the right of the overlapping region to the original V500 data cube. More details of the reduction process can be found in Husemann et al. (2013), García-Benito et al. (2015) and Sánchez et al. (2016). The final data cube, corrected for Galactic extinction, fully covers the optical range from 3700 to 7300 Å.

Preprocessed spectral data (78 × 73 spaxels) are stored in a 3D data cube. We first bin the 2D spectral data to a constrained signal-to-noise ratio (S/N) with the Voronoi binning method (Cappellari & Copin 2003). We adopt the average flux and the standard deviation of the flux as the signal and noise at the wavelength range from 5075 to 5125 Å. This range is not contaminated by emission and absorption lines. Spaxels below a lower S/N limit set by us are stacked together with neighboring spaxels until they satisfy the required minimum S/N. We tested several realizations with S/N = 5, 10, and 20 and finally adopted S/N = 10 in this work. Although higher S/N is better for the fitting, this will reduce the number of spaxels and smooth the details in space. This process produces 1370 Voronoi bins with S/N higher than 10, 80% of which include only one spaxel. Next, we use the pPXF code (Cappellari & Emsellem 2004; Cappellari 2017), a stellar population analysis method, to fit the stellar continuum, during which the emission lines are masked (Hα, [O iii], Hβ, etc.). The code adopts MILES simple stellar population templates (Sánchez-Blázquez et al. 2006; Vazdekis et al. 2010), assuming a Salpeter (1955) IMF and a Calzetti et al. (2000) dust extinction curve, covering 25 population ages (between 0.06 and 15.85 Gyr) and 6 metallicities (log[M/H] = −1.71, −1.31, −0.71, −0.4, 0.0, 0.22). All spectra are fitted in the wavelength range from 3800 to 7200 Å. We check visually all the fitting results and exclude some peculiar fits with stellar velocity dispersions higher than 400 km s⁻¹ and D_n4000 (the ratio of continuum 4000–4100 and 3850–3950 Å, Balogh et al. 1999) higher than 2.0, which results in 1200 spectra analyzed in this work. We emphasize that the elimination of spaxels by the visual inspection process does not significantly affect our results.

Figure 2 shows the distributions of fitted parameters derived from pPXF. Panel (A) represents D_n4000, an indicator for stellar population age (Kauffmann et al. 2003b). The small value of D_n4000 implies that there is ongoing star formation in the central region. A similar trend can be found from the light-weighted age (panel (B)), which is weighted by the fraction of light for each template. Panel (C) represents the stellar velocity map, showing a characteristic rotated disk-shape structure. It is interesting that the equivalent width (EW, positive value represents absorption) of the Hδ line is large along the major axis (panel (D)), suggesting that these regions might have experienced a starburst process several hundred Myr ago.

Zoom In Zoom Out Reset image size

For the emission lines, we subtract the estimated stellar continuum from the observed spectrum and obtain the line flux with the SHERPA IFU line fitting software (SHIFU; García-Benito, in preparation), based on the package CIAO SHERPA (Freeman et al. 2001; Doe et al. 2007). Small deviations with respect to the stellar continuum are taken into account by a first-order polynomial. Single Gaussians have been fitted for the emission lines and the width of the Gaussians were tied for ions of the same element.

PGC 26218 is a star-forming S0 galaxy and displays post-starburst (PSB) regions in the disk. In order to show the locations of regions with different star formation histories, we pick out these regions according to the EW of Hα and Hδ.

A traditional picture of post-starburst galaxies is that these objects have experienced violent starbursts and then rapidly quenched star formation within the last few hundred Myr. Their optical spectra are characterized by strong Balmer absorption lines and weak nebular emission lines. Therefore, the post-starburst galaxies are usually selected based on the deficiency of emission lines (such as Hα or [O ii]) and on strong Hδ absorption lines (Blake et al. 2004; Quintero et al. 2004; Goto 2005, 2007; Poggianti et al. 2009; Vergani et al. 2010). The strong high-order Balmer absorption lines are evidence of the existence of A-type stars, while the absence of emission lines indicates that there is no recent star formation in the past few Gyr.

The traditional definition for the post-starburst phase imposes a rigorous cut on Hα emission line (Balogh et al. 2005; Hogg et al. 2006). The restriction on Hα intensity may miss some post-starburst galaxies that reside in the early evolutionary phase. Recent studies show that the early stage is also important for understanding the evolution of post-starburst galaxies (Yesuf et al. 2014; Rowlands et al. 2015; Alatalo et al. 2016; Pawlik et al. 2018). In this work, the quiescent post-starburst regions (QPSB) are defined as EW (Hδ) > 4 Å and EW (Hα) > −3 Å, while the transitioning post-starburst regions (TPSB) are defined as EW (Hδ) > 4 Å and −10 Å < EW (Hα) < −3 Å. For QPSB, the definition is the same as in Yesuf et al. (2014). The upper limit on EW (Hα) ensures that the TPSB region has only residual star formation. The typical errors for the EW of Hα and Hδ are 1.0 Å and 1.4 Å, respectively. The corresponding regions are marked in panel (A) of Figure 2. We find that the QPSB regions (red crosses) are mainly located outside 1.5 effective radii (R_e ∼ 13''), while the TPSB regions (blue crosses) are mainly concentrated around the center of PGC 26218. Furthermore, we define the regions with EW (Hδ) < 4 Å and EW (Hα) < −3 Å or EW (Hδ) > 4 Å and EW (Hα) < −10 Å as the star-forming (SF) regions. We emphasize that the number of spaxels for QPSB accounts for 20% of the total spaxels of these three regions. This significant amount of QPSB spaxels shows that PGC 26218 has undergone a starburst several hundred Myr ago. Based on the optical features, millimeter observation coordinates (black cross) are selected between the two PSB regions. We will describe the CO(J = 1–0) emission in Section 3.

Figure 3 shows the stacked spectra for SF, TPSB, and QPSB regions, respectively. For the spectrum of SF, it shows the strongest Hα emission and weakest Hδ absorption. The spectrum of TPSB has similar Hδ absorption but stronger Hα emission than that of QPSB. This feature might indicate that PGC 26218 needs a period of transition for the emission lines to fade from TPSB to QPSB.

Zoom In Zoom Out Reset image size

The emission lines may be contaminated by active galactic nuclei (AGNs) and shocks, which would lead to the overestimation of the star formation rate (SFR) when we use Hα as tracer of star formation. To this end, we investigate the excitation mechanisms for SF and TPSB regions according to the traditional BPT diagram (Baldwin et al. 1981, left panel of Figure 4). In this diagram, we impose the S/N of these four emission lines (i.e., Hα, [N ii], Hβ, and [O iii]) to be greater than 3. For the spaxels with Hα and [N ii] S/N larger than 3, we use the relation between the EW (Hα) and [N ii]/Hα (right panel of Figure 4) as the diagnostic diagram of excitation mechanisms (Cid Fernandes et al. 2010). We find no spaxels in the Seyfert region, which indicates that PGC 26218 does not host an AGN. Although some spaxels locate at the composite and LINER regions, the effect from shock is ignorable considering that most of the spaxels locate at the SF region and the contribution from shocks to emission lines could be low in the SF region. Taking into account that all of the SF and TPSB spaxels have EW(Hα) < −3 Å and the percent light contribution (at λ = 5100 Å) of young populations in those spaxels is larger than 4%, we believe that most of the contribution is dominated by ionization due to young stars and only a tiny fraction by hot low-mass evolved stars/post AGB stars (González Delgado et al. 2016; Lacerda et al. 2018). At any rate, our SFR value is an upper limit.

Zoom In Zoom Out Reset image size

In order to estimate the spatially resolved SFR from extinction-corrected Hα emission lines, we calculate the extinction according the Balmer decrement assuming a Cardelli et al. (1989) extinction curve with case B condition for each Voronoi bin. The extinction-corrected SFR is calculated using the formula given by Kennicutt (1998a):

$$\begin{eqnarray}&&\mathrm{SFR}({M}_{\odot }\,{\mathrm{yr}}^{-1})=7.9\times {10}^{-42}L({\rm{H}}\alpha ),\end{eqnarray} \tag{ 1 }$$

where L(Hα) is the extinction-corrected Hα luminosity. For the total extinction-corrected SFR, we construct a spectrum by stacking all spaxels in Voronoi bins of SF and TPSB regions. We use the pPXF code again to model the stellar continuum. For the continuum-subtracted emission-line spectrum, we simultaneously fit Hα, [N ii], Hβ, and [O iii] lines using the single Gaussian model (see Figure 5) and calculate the extinction. It is worth noting that the MPA/JHU SDSS DR7 catalog (Brinchmann et al. 2004) provided the aperture-corrected SFR of PGC 26218 (SFR ∼ 0.20 M_⊙ yr⁻¹) based on a Kroupa (2001) IMF, whereas we adopt a Salpeter (1955) IMF and find the total SFR is ∼0.28 M_⊙ yr⁻¹. The conversion between the two IMFs is Kroupa IMF × 1.5 ≈ Salpeter IMF. Therefore our SFR is consistent with that given by Brinchmann et al. (2004).

Zoom In Zoom Out Reset image size

The inclination and axial ratio (see Table 2) derived from GALFIT (Version 3.0.5, Peng et al. 2002) show that the projection effect may affect the accuracy of surface densities. We calculate the inclination based on

$$\begin{eqnarray}&&i=\arccos \sqrt{\displaystyle \frac{{q}^{2}-{q}_{0}^{2}}{1-{q}_{0}^{2}}}\end{eqnarray} \tag{ 2 }$$

where q = b/a, the ratio of minor semiaxis and major semiaxis. We adopt q₀ = 0.25 (Sandage et al. 1970) instead of 0.2 (Hubble 1926) to correct the inclination for classical S0 galaxies. We find the inclination-corrected star formation rate surface density (Σ_SFR), also referred to as the intensity of the star formation, is lower by 0.3 dex than that of non-inclination-corrected Σ_SFR. However, the correction for the inclination will lead to an equal decrease in molecular gas surface density (Σ_gas). Thus, the inclination correction will only cause the parameters to move to the lower left in the K–S law, but will not significantly affect the results shown in Figure 9.

Table 2.  A Summary of the Physical Properties for PGC 26218

PGC 26218
Av (mag)	0.76 ± 0.06
log SFR (M☉ yr−1)	−0.55 ± 0.01
log ΣSFR (M⊙ yr−1 kpc−2)	−1.39 ± 0.01
Inclination (deg)	64.5
b/a	0.49
12 + log (O/H)	8.53 ± 0.13
SCO (J = 1–0) Δv (Jy km s−1)	6.78 ± 0.75
$\mathrm{log}\,{L}_{\mathrm{CO}\,(J=1-0)}^{{\prime} }$ (K km s−1 pc2)	6.96 ± 0.75
log ${M}_{{{\rm{H}}}_{2}}$ (M☉)	7.60 ± 0.15
log Σgas (M⊙ pc−2)	1.30 ± 0.14
log Σ1 (M☉ kpc−2)	8.25

Download table as:  ASCIITypeset image

The gas-phase metallicity of PGC 26218 is derived based on the Bayesian method⁵ (Blanc et al. 2015) with the Levesque et al. (2010) photoionization model, which allows us to input arbitrary sets of strong nebular emission lines to infer the probability density functions of gas-phase metallicity and ionization parameter. We use the spectra that are used to compute the SFR to estimate the gas-phase metallicity, which is ∼1σ below the relation derived by Tremonti et al. (2004) while consistent with the relation given by Sánchez et al. (2013) within the error range.

We use the average gas-phase metallicity (see Table 2) to derive the surface mass density of atomic gas (Σ_{H i}) according to the empirical formula given by Schruba et al. (2018). In this scaling relation, the optically thin conversion is adopted and a factor of 1.36 is taken into account to include the heavy elements. Furthermore, the effects of diffuse H i on the saturation column density is neglected. However, the assumption does not significantly affect Σ_{H i}.

PGC 26218 was observed with NOEMA located in the south of the French Alps on 2018 June 16 (Project S18BN001. PI: Xue Ge). The source was observed for 2.5 hr with the compact D configuration. This configuration is best suited for deep integration and coarse mapping (resolution ∼37 at 100 GHz and ∼16 at 230 GHz) and provides the lowest phase noise and highest sensitivity. The source 3C 273 is chosen to calibrate the bandpass, while the sources J0851+202 and J0923+282 are chosen to calibrate the phase. The CO(J = 1–0) rest frequency is 115.271 GHz (the redshifted frequency is 114.643 GHz). We observed the emission line with the receivers set to a 3 mm band tunable sky frequency between 70.4 and 119.9 GHz. The receiver band has dual-polarization capabilities and each of the two polarizations delivers a bandwidth of 7.744 GHz in the lower sideband and upper sideband simultaneously. Each sideband contains two adjacent basebands of ∼3.9 GHz width, called the inner and outer baseband. The spectral resolution is ∼2 MHz throughout the wide sidebands.

The CO(J = 1–0) line data are calibrated by CLIC, a module of the available GILDAS software package. We choose the channels from −860 to 400 km s⁻¹, which is a wide enough range to analyze the CO(J = 1–0) line. The final channels are smoothed to 10 km s⁻¹. The cleaning of the line image is done by the mapping module of GILDAS and the resulting channel map is shown in the left panel of Figure 6. The synthesized beam size is 1062 × 476. The uncertainty of calibration in the flux is ∼20% and the 1σ uncertainty is about ∼5.6 mJy/beam. The reference center of the map is set to be the same as the optical center so that we can compare the relative positions of Hα and CO(J = 1–0); (see the right panel of Figure 6). We can see that the emission of CO(J = 1–0) is mainly concentrated at SF and TPSB regions and basically overlaps the emission of Hα. We do not detect CO(J = 1–0) emission in the QPSB regions, but weak emission exists in the TPSB region. In addition, we find that there is an offset between the peaks of the optical and the millimeter CO emission (see the detailed description in Section 4.4).

Zoom In Zoom Out Reset image size

We extract the CO(J = 1–0) emission of each spaxel with enough S/N (3σ) at the peak and construct a stacked spectrum (Figure 7). The interesting point is that the CO(J = 1–0) line shows asymmetry in its profile. We fit the spectrum with two Gaussians. One of them represents the redshift component, while the other represents the rotation component. It is found that one of the Gaussians shows a redshift of ∼46 km s⁻¹, which might indicate the inflow of gas, while the other does not. The asymmetrical profile does not disappear even though we stack the spectra within an effective radius. The flux density of CO(J = 1–0) emission is estimated from the best-fit model. We use the relation given by Solomon & Vanden Bout (2005) to calculate the luminosity of CO(J = 1–0):

$$\begin{eqnarray}&&{L}_{\mathrm{CO}}^{{\prime} }=3.25\times {10}^{7}\,{S}_{\mathrm{CO}}{\rm{\Delta }}v\,{\nu }_{\mathrm{obs}}^{-2}\,{D}_{L}^{2}\,{\left(1+z\right)}^{-3}\end{eqnarray} \tag{ 3 }$$

where S_COΔv, ν_obs, and D_L are the CO integrated flux density in units of Jy km s⁻¹, the observing frequency in GHz, and the luminosity distance in Mpc, respectively.

Zoom In Zoom Out Reset image size

The mass of molecular hydrogen (${M}_{{{\rm{H}}}_{2}}$) is estimated using following formula:

$$\begin{eqnarray}&&{M}_{{{\rm{H}}}_{2}}={\alpha }_{\mathrm{CO}}\times {L}_{\mathrm{CO}}^{{\prime} }\end{eqnarray} \tag{ 4 }$$

where α_CO is the conversion factor between CO(J = 1–0) and H₂. Combining the Σ_{H i} derived in Section 2.4, the inclination-corrected Σ_gas (${{\rm{\Sigma }}}_{{\rm{H}}{\rm{I}}+{{\rm{H}}}_{2}}$) is estimated by adopting α_CO = 4.3 M_☉ (K km s⁻¹ pc²)⁻¹, for the inner disk of Galaxy (Bolatto et al. 2013).

Figure 8 shows the map of SFR in logarithmic space for PGC 26218. It is found that star formation occurs mainly in an effective radius, which is consistent with the distribution of CO emission. As mentioned in Figure 7, the redshift component might indicate gas inflow, which could have led to the star formation in the central region.

Zoom In Zoom Out Reset image size

Previous studies have shown that the stellar mass surface density within the central 1 kpc, Σ₁, strongly correlates with the star formation (Cheung et al. 2012; Fang et al. 2013; Luo et al. 2019). In order to compare our result to Fang et al. (2013), we also adopt the stellar mass from the MPA-JHU DR7 catalog ($\mathrm{log}{M}_{* }=9.15$) and a similar method to estimate Σ₁. Here, we estimate Σ₁ from the NASA-Sloan Atlas catalog providing the surface brightness profiles in a series of angular sizes. The extinction-corrected and k-corrected surface brightness profiles within 1 kpc in i-band are used to convert to Σ₁ according to the mass-to-light ratio of the i-band (see Fang et al. 2013 for detailed calculations). The value of log Σ₁ for PGC 26218 is 8.25 M_⊙ kpc⁻² (see Table 2), which is 0.34 dex (∼2σ) below the best-fitting relation of Σ₁–M_* obtained from the combination of blue and green valley galaxies (Fang et al. 2013). We also extract the mass value within 1 kpc from the stellar population analysis method and find that the value of Σ₁ does not change significantly. The position of PGC 26218 in the Σ₁–M_* relation suggests that the central star formation might not be suppressed significantly. With the current SFR, the remaining molecular gas can last for about 0.1 Gyr.

The SFMS (Elbaz et al. 2007; Salim et al. 2007) shows the relationship between stellar mass and SFR. S0 galaxies are often thought to have an absence of noticeable star formation regions. To investigate whether PGC 26218 resides in the SFMS, we put the source on the SFMS relation given by Elbaz et al. (2007), who obtained SFR and stellar mass from the MPA-JHU DR4 catalog.⁶ It is noted that the SFMS in Elbaz et al. (2007) was calibrated based on a Salpeter (1955) IMF. So, we adopt here $\mathrm{log}{M}_{* }=9.32$ to investigate the position of PGC 26218 in the SFMS. We find PGC 26218 deviates slightly from the best-fitting relationship (slightly below the SFMS), but it follows the SFMS relation well within the error range. Speagle et al. (2014) gave a redshift-dependent SFMS relation using a compilation of many studies from the literature. We find that PGC 26218 still follows the redshift-dependent SFMS relation considering 0.2 dex scatter.

The tight relationship between Σ_gas and Σ_SFR (i.e., the K–S law) has been found in different types of star-forming galaxies (Schmidt 1959; Kennicutt 1998b; Kennicutt & Evans 2012). Figure 9 displays the K–S law for ETGs, star-forming, and starburst galaxies. We find that the ETGs (green points, most of them are S0 galaxies) in Davis et al. (2014) deviate from this relationship, while PGC 26218 (red point) basically obeys the K–S law, which implies that some mechanisms enhance its SFR. It is noted that Σ_SFR is computed in the same region as the CO(J = 1–0) measurement (see Section 3.2). We find that the result will not change if we compute Σ_SFR and Σ_gas in an effective radius. Although the fraction of molecular gas for PGC 26218 (∼2%) is slightly lower than that of Davis et al. (2014) (∼3%), the star formation efficiency in PGC 26218 might be higher. Wei et al. (2010) presented the relationship between molecular gas and star formation in low-mass ellipticals/S0 galaxies. They found that most blue-sequence ellipticals/S0 galaxies show higher Σ_SFR at fixed ${{\rm{M}}}_{{{\rm{H}}}_{2}}$, which is similar to local dwarf galaxies. In addition, the fraction of molecular gas is a factor of 2 higher than that of PGC 26218. Kokusho et al. (2017) studied the star formation of the same sample as Davis et al. (2014) in the ${\mathrm{ATLAS}}^{3{\rm{D}}}$ survey. They found that the local ETGs have a similar star formation efficiency to the star-forming galaxies and follow the K–S law. The difference between them might be attributed to the methods of calculating SFR. Davis et al. (2014) calculated the SFR via the far-UV and 22 μm, while Kokusho et al. (2017) did so via spectral energy distribution fitting. Colombo et al. (2018) combined the CO observations made by the Combined Array for Millimeter-wave Astronomy from Extragalactic Database for Galaxy Evolution survey (Bolatto et al. 2017) with the 2D-spectroscopic observations made by CAHA from the Calar Alto Legacy Integral Field Area survey (Sánchez et al. 2012). They calculated Σ_SFR from the Hα map and investigated the dependence of Σ_SFR on Hubble types. They found that S0 galaxies have lower Σ_SFR than that of spirals, which is inconsistent with PGC 26218.

Zoom In Zoom Out Reset image size

Within the region of CO(J = 1–0) emission, we compare the kinematics of stars, Hα, and CO(J = 1–0). Figure 10 shows the distributions of the line-of-sight velocities for these three components. The major axis of stellar, Hα, and CO(J = 1–0) velocity fields (dashed lines in Figure 10) are measured using the KINEMETRY code (Krajnović et al. 2006). We define the position angles as the counterclockwise angle between north and a line that bisects the velocity field of gas or stars on the receding side. We find that the major axis of stellar, Hα, and CO(J = 1–0) velocities are 156° ± 2°, 178° ± 4°, and 133° ± 10°, respectively. The difference in rotation axes between Hα and CO(J = 1–0) is larger than 45°. In addition, we find that the velocity of CO(J = 1–0) displays an asymmetry, which indicates that the gas has a motion relative to the stellar components. In general, the kinematic misalignment between ionized gas and CO(J = 1–0) implies that the fuel for star formation might come from the external environment.

Zoom In Zoom Out Reset image size

George (2017) found that star-forming ellipticals/S0 galaxies display the features of gas-rich minor mergers that fuel the central star formation. The above results give us the motivation to explore the mechanisms causing PGC 26218 to undergo a burst of star formation.

The morphological disturbations of galaxies can give us clues into the origins of nuclear activities. In order to investigate the central structure of PGC 26218, we simply use a disk component (Sérsic index = 1) to model the galaxy disk of PGC 26218 in the r-band (Figure 11). We can clearly see from the residual image that PGC 26218 displays a peculiar structure at its center. It is worth noting that the peculiar structure in the residual image does not disappear even though we use two components (i.e., bulge and disk components) to decompose the r-band image. The Sérsic indices of bulge and disk are 1.2 and 1.0, respectively. In order to highlight the structure of the central region, we compute only the disk component assuming a pure disk model. Figure 12 highlights the disk-subtracted residual image. We find PGC 26218 shows extreme morphological peculiarities that may be triggered by a galactic merger. We mark the positions with the highest flux compared with surroundings as core a, core b, and core c in the r-band, respectively. It is found that the projected distances between core a and core b, core a and core c, and core b and core c are 32 (0.35 kpc), 60 (0.65 kpc), and 28 (0.31 kpc), respectively. The peak of CO(J = 1–0) nearly overlaps core c and it lies ∼0.6 kpc from the center of Hα (i.e., core a). However, we should keep in mind that the offset is comparable to the resolution of CO(J = 1–0). Higher resolution observations are needed to resolve this question.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Actually, many previous studies have presented the presence of multiple optical nuclei in PGC 26218. Mazzarella & Boroson (1993) presented the properties of over 100 Markarian galaxies with multiple nuclear structures or peculiar morphologies. They used the software IMAGES to measure the structure of nuclei and found that PGC 26218 (Mrk 1230 in their paper) shows extended nuclear structure along the major axis. Nordgren et al. (1995) studied the morphologies and kinematics of 16 Markarian galaxies. Their results suggested that PGC 26218 has three nuclei and the farthest and nearest projected distances between two nuclei are 1.2 kpc and 0.5 kpc, respectively. In addition, PGC 26218 is also included in Gimeno et al. (2004) as a double nucleus case. They found that the separation between the two nuclei is about 0.4 kpc. We note that if the two nuclei depicted by Gimeno et al. (2004) are the a and b components in our Figure 12, then the projected separation estimated by us agrees well with theirs.

There are several scenarios that can be used to explain the occurrence of central star formation activity. (i) Internal secular evolution, such as the turbulence of a bar (Pfenniger & Norman 1990). The presence of a bar causes the gas to inflow to the center, thus inducing the starburst and fueling the mass concentration. In this scenario, the vigorousness of the starburst induced by the bar depends on the stellar mass. Carles et al. (2016) suggested that massive (${M}_{* }\gt 2\times {10}^{10}{M}_{\odot }$) barred galaxies more easily consume gas by bar-induced starbursts than low-mass barred galaxies. We do not expect that the starburst that occurred several hundred Myr ago was triggered by a bar because we do not find a prominent bar component in such a low-mass S0 galaxy. (ii) Disk instability. The instability of a gas-rich disk may lead to the formation of star-forming clumps, which can lose angular momentum due to interactions and fall toward the center of the galaxy (Elmegreen et al. 2008; Dekel et al. 2009). The observational evidence shows that the dissipational processes triggered by disk instabilities are expected to occur in massive galaxies in earlier cosmic times where there is a higher gas fraction than that of local universe (Bouché et al. 2007; Daddi et al. 2010). However, we do not expect a large amount of disturbed gas to be conserved in such a low-mass, low-z S0 galaxy because the current fraction of molecular gas mass is only ∼2% of the stellar mass and the galaxy has relatively regular rotations in its stellar and gas disks. (iii) Major and minor mergers. Cosmological hydrodynamic simulations show that the frequency of major mergers declines with the cosmic time (Maller et al. 2006). They found that the average merger rate of massive galaxies is 0.054 Gyr⁻¹ at z ∼ 0.3, and is 0.018 Gyr⁻¹ for low-mass galaxies. It is believed that major mergers may be too destructive to preserve the inner components and regular disks. Minor mergers are expected to significantly increase the SFR (Saintonge et al. 2012; Kaviraj 2014) although they are not as violent as major mergers. At low redshift, the incidence of minor mergers is higher than that of major mergers and the minor mergers can build bulges without destroying disk structure. As suggested by Eliche-Moral et al. (2011), minor mergers can explain the existence of multiple inner components in unbarred galaxies, although this mechanism is more complex than other processes (such as bars and ovals, Kormendy & Kennicutt 2004). Eliche-Moral et al. (2018) have studied whether galaxy mergers can reproduce the features of S0-like remnants based on the GalMer simulated database. They suggested that the mergers can result in relaxed morphologies and inner subcomponents, such as ovals, lenses, and compact sources. Furthermore, the relics are more durable in minor mergers than in major mergers. For PGC 26218, a minor merger leads to the slow rearrangement of gas from the disk to the center of galaxy and this progress has not yet been completed (as indicated by the multiple nuclear structures we see).

Bellstedt et al. (2017) studied stellar kinematics of four S0 galaxies with ${M}_{* }\gtrsim {10}^{10}{M}_{\odot }$ utilizing the DEIMOS instrument on the Keck telescope. They found that these S0 galaxies generally resemble spiral progenitors more than merger remnants. Although they suggested that these S0 galaxies are likely formed via faded disks, merger events cannot be ruled out for the formation of S0 galaxies. Fraser-McKelvie et al. (2018) investigated the relative importance of gas stripping and mergers using 279 S0 galaxies from the MaNGA (Mapping Nearby Galaxies at APO) survey (Bundy et al. 2015). They found that the bulges of low-mass S0 galaxies are almost always younger than their disks, which is consistent with PGC 26218. However, they contributed it to bulge rejuvenation or disk fading. Temi et al. (2009) studied star formation using a sample of local ETGs from the Spectroscopic Areal Unit for Research on Optical Nebulae survey (Bacon et al. 2001). They found that the cold gas for star formation in S0 galaxies is created via stellar mass loss. A similar result has been found by Shapiro et al. (2010). Furthermore, they also suggested that the minor mergers might be building up the bulges of red sequence S0 galaxies. Kaviraj et al. (2009, 2010) studied the importance of minor mergers in low-level star formation of ETGs. They found that minor mergers are the main mechanism driving star formation at low redshift. Méndez-Abreu et al. (2019) studied the star formation of local S0 galaxies with M_* ≳ × 10¹⁰ M_☉ from the CALIFA survey. Their results show that the star formation in bulges is significantly lower than that in disks. This contradiction to our result might be attributed to the type of bulge because a larger bulge is more likely to quench star formation in the central region via the morphological, but related with mass, quenching scenario. This can be also understood by differences of peculiar objects (as their NGC 3773 case) from the general trend of low-mass objects (${M}_{* }\lt {10}^{10}{M}_{\odot }$) as seen in the dispersion in the derived stellar population properties in the low-mass range (García-Benito et al. 2017).

In summary, PGC 26218 is a good case to connect the relevance between minor mergers in the past and star formation in the present, together with the evolutionary trajectories of S0 galaxies. We suggest by combining the results in this work with those of previous studies that PGC 26218 could have undergone a gas-rich minor merger, which triggered the starburst at the outskirts a few hundred Myr ago and the central star formation activity as indicated by the multiple nuclear structures.