Formation of S0s in extreme environments II: the star-formation histories of bulges, discs and lenses

Different processes have been proposed to explain the formation of S0s, including mergers, disc instabilities and quenched spirals. These processes are expected to dominate in different environments, and thus leave characteristic footprints in the kinematics and stellar populations of the individual components within the galaxies. New techniques enable us to cleanly disentangle the kinematics and stellar populations of these components in IFU observations. In this paper, we use buddi to spectroscopically extract the light from the bulge, disc and lens components within a sample of 8 S0 galaxies in extreme environments observed with MUSE. While the spectra of bulges and discs in S0 galaxies have been separated before, this work is the first to isolate the spectra of lenses. Stellar populations analysis revealed that the bulges and lenses have generally similar or higher metallicities than the discs, and the $\alpha$-enhancement of the bulges and discs are correlated, while those of the lenses are completely unconnected to either component. We conclude that the majority of the mass in these galaxies was built up early in the lifetime of the galaxy, with the bulges and discs forming from the same material through dissipational processes at high redshift. The lenses, on the other hand, formed over independent timescales at more random times within the lifetime of the galaxy, possibly from evolved bars. The younger stellar populations and asymmetric features seen in the field S0s may indicate that these galaxies have been affected more by minor mergers than the cluster galaxies.


INTRODUCTION
The lenticular, or S0, classification is often used for galaxies whose morphology is intermediate between ellipticals and spirals, such that they display the same discy structure seen in spirals but with Contact e-mail: evelynjohnston.astro@gmail.com the redder colours and old stellar populations typically found in ellipticals. S0s are now considered, in most cases, to be evolved spirals in which the star formation has been truncated and the spiral arms have faded (e.g. Bedregal et al. 2006;Moran et al. 2007;Laurikainen et al. 2010;Cappellari et al. 2011;Prochaska Chamberlain et al. 2011;Kormendy & Bender 2012;Johnston et al. 2014). Further evidence for this scenario lies in the morphology-density relation, which has shown that the frequency of S0s increases with increasing environmental density while the frequency of spirals decreases proportionally (Dressler 1980;Dressler et al. 1997;Cappellari et al. 2011). Many mechanisms have been proposed to explain this transformation, such as the gas being stripped out via rampressure stripping (Gunn & Gott 1972), harassment (Moore et al. 1998) or strangulation (Larson et al. 1980), or used up in a rapid star formation event following a minor merger (Ponman et al. 1994;Arnold et al. 2011) or successive minor mergers (Bekki & Couch 2011). Furthermore, simulations have shown that major mergers can create a bulge around which a disc is able to form to ultimately create an S0 galaxy (Spitzer & Baade 1951;Tapia et al. 2017;Diaz et al. 2018;Eliche-Moral et al. 2018;Méndez-Abreu et al. 2018), or even a spiral galaxy (Fabricius et al. 2014). All of these processes require an external influence, such as an interaction with a neighbouring galaxy or the intracluster medium, but the existence of S0s in isolated environments indicates that other evolutionary processes must be possible. For example, Eliche-Moral et al. (2013) proposed that S0s could simply be remnants of spirals after the gas reservoirs were exhausted and the spirals arms have faded, while Saha & Cortesi (2018) suggested that low-mass S0s may have formed through disc instabilities and fragmentation.
Many studies have attempted to determine the conditions under which these different processes are more dominant in transforming a spiral galaxy into an S0, such as dependencies on the mass of the galaxy and the local environment. For example, evidence of a mass dependence in the formation of S0s has been found by Fraser-McKelvie et al. (2018b). They studied the stellar populations within the bulges and discs of a sample of S0s from the MaNGA survey, and concluded that S0s with masses > 10 10 M transformed via morphological or inside-out quenching, while lower mass counterparts have undergone bulge rejuvenation or disc fading. Another study by Domínguez Sánchez et al. (2020) concluded that S0s with masses > 10 10 M have formed mainly through mergers, while those with < 10 10 M are closer to spirals in their formation, and thus in-situ star formation has played a more significant role than mergers.
The first paper in this series (Coccato et al. 2020, hereafter Paper I) investigated the effect of both the mass and environment on the formation of S0s through studies of their kinematics and global stellar populations. The kinematics revealed that the cluster galaxies are more rotationally supported while the isolated galaxies are more dispersion supported, thus suggesting different formation mechanisms, i.e. through gas removal (e.g. starvation, ram-pressure stripping) and minor mergers respectively. Another study of the kinematics of S0 galaxies in the SAMI Survey (Green et al. 2018) by Deeley et al. (2020) also found the same result. The integrated stellar populations studied in Paper I, on the other hand, revealed no significant differences as a function of environment, but a tentative correlation was found with the stellar mass of the galaxies such that more massive galaxies contained older and more-metal-rich stellar populations. However, this stellar populations analysis was limited to the global properties within 1R e , where the bulge dominates the light, and so a more detailed stellar populations analysis across the galaxies is therefore needed to better understand this trend.
The different transformation processes would leave characteristic signatures within the physical structure and stellar populations of the different components within the final galaxy, and so by studying these properties we can better understand how S0s formed. S0s were historically seen as featureless discs with, in many cases, a spheroidal bulge dominating the light at the core. However, more recent studies have found that their structures are not so straightforward. For example, bulges are no longer seen as simple spheroids, and instead can be classified as classical or pseudo-bulges based on the steepness of their surface brightness profiles and whether they are dissipationally or rotationally supported structures (Kormendy & Kennicutt 2004).
Additionally, S0s are now known to host many distinct substructures, such as bars, ovals and lenses (e.g. Laurikainen et al. 2005Laurikainen et al. , 2009). Kormendy (1979) defined a lens as an independent component with distinct surface brightness profiles, which are morphologically intermediate between bulges and discs with their flatter surface brightness profiles and sharp edges. Ovals have similar surface brightness profiles to lenses, but are more elongated, while bars display even higher ellipticities and can be detected through distortions in the kinematics across the galaxy (e.g. Bureau & Freeman 1999;Merrifield & Kuijken 1999;Gadotti & de Souza 2005). The true nature of lenses and ovals is still uncertain, although many current theories identify them as either remnants of dissolved bars (Kormendy 1979) or artefacts from major mergers that are associated to stellar halos or embedded inner discs (Eliche-Moral et al. 2018). It was originally thought that lenses and ovals were related to the host galaxy morphology, where lenses are found in galaxies of S0 and S0/a morphologies while ovals reside in later type disc galaxies. However, Laurikainen et al. (2006) found that S0s can host both lenses and ovals, reflecting that the distinction between these classifications based on the galaxy morphology is not so straightforward.
While lenses have been seen in galaxies of all discy morphologies, they appear most frequently in S0s (Laurikainen et al. 2009), both barred and unbarred (Laurikainen et al. 2005(Laurikainen et al. , 2007. The fraction of S0s reported to display lenses and ovals varies from ∼ 15% (Nair & Abraham 2010) up to 97% (Laurikainen et al. 2009). Buta et al. (2010) proposed that lenses and ovals are the remnants of bars that have become disrupted during the transformation from spirals to S0, and thus are created through secular evolution. As further evidence towards this scenario, the frequency of lenses from spirals to S0s increases proportionally to the decrease in the frequency of bars between the same two morphologies (e.g. Laurikainen et al. 2009). Additionally, bars in S0s have been seen to be in general more evolved than in spirals, being longer, more massive and more often with ansae morphologies (bars with bright enhancements at each end; Elmegreen & Elmegreen 1985;Laurikainen et al. 2007;Martinez-Valpuesta et al. 2007;Díaz-García et al. 2016). Gas infall has also been proposed by Pfenniger & Norman (1990) and Laurikainen et al. (2010) in connection with the secular evolution theory, whereas Eliche-Moral et al. (2018) showed that major mergers can also create these structures using data from the dissipative merger simulations of the GalMer database (Chilingarian et al. 2010).
Many imaging surveys have been carried out to study the structural parameters of the different structures within S0s. For simplicity, the larger, more statistical surveys have focussed on modelling the galaxies using a single or double light profile, where the single profile was usually a Sérsic profile and gives information on the global structure of the galaxy while double profiles usually model an exponential disc with either a de Vaucouleurs or Sérsic bulge (e.g. Simard et al. 2011;Mendel et al. 2014;Vika et al. 2014;Head et al. 2014;Bottrell et al. 2019). However, failures to account for all the components within a galaxy can affect the results. For example, excess light from a lens or bar may be attributed to the bulge, thus skewing the B/T light ratio to higher values and affecting the measurements for the effective radius and Sérsic profile for that component. Additionally, Head et al. (2015) showed that discs are not always perfectly exponential, but instead show a wide range of Sérsic indices. As a result, limiting the fit to the discs to this profile may also result in light in the inner regions of the galaxy being mis-attributed to the bulge or disc. As a result, a few recent studies have expanded their fits beyond such simple models for their galaxies. For example, the Near-IR S0 Galaxy Survey (NIRS0S, Laurikainen et al. 2011) modelled K s -band images of S0s and classified their morphologies in terms of bulges, discs, bars, lenses, barlenses etc, and Head et al. (2015) modelled g, r and i-band images of galaxies in the Coma Cluster with up to three Sérsic components to take into account lenses, bars and breaks or truncations in the discs. This latter study found that the fraction of pure Sérsic profiles increases towards the centre of the cluster while the fraction of exponential discs is higher in the outskirts, but found no significant variation in the fraction of multi-Sérsic systems with radius within the cluster.
In order to better understand the different formation scenarios of bulges, discs and lenses in S0s, and thus their roles in the transformation from spiral to S0, one needs to spectroscopically separate the light from each of these components. With the development of Integral Field Unit (IFU) spectrographs with wide fields of view and high spatial resolution, it is now possible to model the light profile of a galaxy in each image slice in order to cleanly extract the spectrum for each component included in the model with minimal contamination from the superposition of light from the other structures. One approach to carry out this modelling and extraction of spectra is Bulge-Disc Decomposition of IFU data . In this paper, we apply this novel approach to a sample of 8 S0s observed with the Multi-Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010) at the Very Large Telescope (VLT) in order to identify the physical structures within these galaxies and derive estimates of their stellar populations and star formation histories from their spectra. While several studies in the literature have extracted spectroscopic stellar populations for bulges and discs of S0s (e.g. Sil'chenko et al. 2012;Johnston et al. 2012Johnston et al. , 2014Coccato et al. 2018;Tabor et al. 2017Tabor et al. , 2019Méndez-Abreu et al. 2019a,b;Oh et al. 2020), this work achieves the first such analysis of lenses. Through this analysis, we use the combined photometric and spectroscopic information to build up a better understanding of how each component has formed and evolved, and their role in shaping the galaxies we see today.
This paper is laid out as follows: Section 2 describes the sample selection and data reduction; Section 3 studies the the stellar populations across each galaxy; Section 4 outlines the spectroscopic decomposition technique used to cleanly model and extract the spectra of each component; Section 5 discusses the components included in the models; Section 6 presents the stellar populations analysis for each component through line strength analysis and full spectral fitting, and finally Section 7 discusses the results and presents the conclusions. Throughout this paper we assume a Hubble constant of H o = 70 ± 1.3 km s −1 Mpc −1 (Carroll et al. 1992).

DATA SAMPLE AND REDUCTION
This study uses observations of a sample of 8 S0 galaxies that were observed with MUSE in Wide-Field Mode without adaptive optics as part of ESO Programme 096.B-0325 (PI: Jaffé). The sample consists of four galaxies from the Centaurus Cluster and four galaxies identified as lying in isolated environments according to the 2MASS Isolated Galaxies Catalogue (2MIG, Karachentseva et al. 2010). The galaxies were selected to have morphological Hubble types of −3.5 < T < −0.5 from LEDA, clear bulge+disc structure with no evidence of bars or tidal disruptions, and to have inclinations above 40 • to reduce the chances of including misclassified ellipticals in the sample. The Centaurus Cluster galaxies lie at a distance of 48.7 Mpc (z ∼ 0.011) while the isolated galaxies have distances of 60−160 Mpc (z ∼ 0.014−0.037), and all galaxies have masses within the range log(M/M ) = 10.6 − 11.6. Full details of the observations and data reduction can be found in Paper I.

RADIAL STELLAR POPULATIONS
In Paper I, the global mass-weighted stellar populations were measured within 1 R e for each galaxy. Due to the small sample, no clear differences were detected as a function of environment. However, a dependence on mass was detected, where higher mass galaxies tend to have higher ages and metallicities, and built up 50% of their total mass earlier on in their lifetimes. Similar trends were found for ETGs in McDermid et al. (2015), and for bulges and discs of S0s in Fraser-McKelvie et al. (2018b).
In this section, the stellar populations across each galaxy will be studied in more detail. As described in Paper I, six of the galaxies in this sample contain Hα emission, either concentrated in the core or distributed throughout the galaxy. The results presented in this paper focus on the analysis of the underlying stellar populations as measured from the stellar continuum and absorption features.

Luminosity-Weighted Ages and Metallicities
Since the light from a galaxy is dominated by the light emitted from the youngest stars, even if they contribute towards only a small fraction of the total galaxy mass, the luminosity-weighted stellar populations are often more of an indicator of how long ago the most recent episode of star formation was truncated. Thus, through the analysis of the luminosity-weighted stellar populations across each galaxy, one can determine where in the galaxy the star formation was truncated most recently, which in turn can provide clues as to the process which caused the gas to be stripped out and/or used up.
The datacubes for each galaxy were first binned using the Voronoi Tesselation technique of Cappellari & Copin (2003) in order to extract spectra from across the galaxy with a higher Signalto-Noise (S/N> 50). The Lick indices for the hydrogen (Hβ), magnesium (Mgb) and iron (Fe5270 and Fe5335) were then measured in these spectra using the definitions of Worthey et al. (1994) with the software of Cardiel (2010), and used as age (Hβ) and metallicity ([MgFe] = Mgb (0.72 × Fe5270 + 0.28 × Fe5335)) indicators. This software estimated uncertainties in these measurements through a series of simulations using the errors in the lineof-sight velocities and the S/N measured from the spectrum itself, details of which can be found in Cardiel et al. (1998). Where emission lines were detected in the spectra, the affected spectra were modelled again using to fit both the stellar absorption and gas emission spectra. The best fit model for the emission lines was then used to correct for the emission in the original spectrum.
The measurements were plotted onto the Single Stellar Populations (SSP) model grids of Vazdekis et al. (2010) in order to convert them into estimates of the luminosity-weighted stellar populations. These models grids use the MILES stellar library (Sánchez-Blázquez et al. 2006) Figure 1. Age, metallicity and [α/Fe] ratio maps for each galaxy, colour coded according to the colour bar at the top of each column. The isolated galaxies are in the left columns while the cluster members are in the right columns, and the label for each galaxy is in the age map. The scale bar at the bottom-right of the [α/Fe] shows 10 and the corresponding physical size in kpc, and the log 10 of the mass is given in the bottom left of these maps. All maps are orientated with North up and East left and are plotted to the same scale.
information that normally occurs when degrading the data to match lower-resolution models.
The results for the luminosity-weighted ages and metallicities are plotted as maps in Fig. 1, and as a function of radius along the major axis in Fig. 2 after correcting for the inclination of the galaxy. The error bars in Fig. 2 represent the effect on the line strengths due to uncertainties in the line-of-sight velocities and the S/N of each spectrum. The ages show no clear systematic trend, with most galaxies showing a flat age gradient and generally very old ages where small differences are hard to detect. The cores of many galaxies show hints of younger stellar populations in their cores (e.g. 2MIG 131, CCC 43 and CCC 137) or immediately outside their core (2MIG 445), which may indicate that the final episode of star formation within these galaxies occurred in their innermost regions. However, in most cases, even these younger stellar populations are above 10 Gyr, where age differences leave only very small effects on the line strengths, and so one must take care when interpreting these trends.
Two galaxies, 2MIG 445 and CCC 158, show significantly younger stellar populations than the rest of the sample. While S0s are generally observed to contain older stellar populations, younger luminosity-weighted stellar populations have been detected in sev-eral studies (e.g. Kuntschner & Davies 1998;Poggianti et al. 2001;Denicoló et al. 2005). The age plot for CCC 158 shows no strong gradient, but the map appears to show that the younger stellar populations lie along the major-axis in the inner regions of the galaxy, while the older stellar populations are seen in the outermost bins. 2MIG 445 on the other hand shows a ring of younger stellar populations circling a core of intermediate-age stars with an outer radius of ∼ 12 (∼ 9.1 kpc), with intermediate ages throughout the rest of the galaxy. These results indicate more recent star formation within these galaxies, potentially indicating more recent transformation from spirals, and the properties of these galaxies will be revisited throughout this paper.
The metallicity maps and plots show evidence that all galaxies have higher metallicities in their cores than in the outskirts. Similar metallicity trends have been seen before in studies of S0s, e.g. (Ogando et al. 2005;Bedregal et al. 2011;Prochaska Chamberlain et al. 2011). Zibetti et al. (2020) suggested that the steeper negative metallicity trend in the inner regions of galaxies was created in the early dissipative collapse in the core of the galaxy, while the flatter gradients in the outer regions are an effect of the accretion of exsitu stars from low-metallicity satellites. Interestingly, the maps in Fig. 1 show that in CCC 158 and 2MIG 445 the regions of higher metallicity coincide with the youngest stellar populations, indicating that the most recent episode of star formation was fuelled by gas that was likely enriched by star formation throughout the rest of the galaxy.

Luminosity-Weighted α-enhancement
While the luminosity-weighted stellar populations presented in the previous section give details about the most recent episode of star formation within each component, the α-element enhancement can provide information on the star-formation timescale of that event and the origin of the gas that fuelled it. The α-enhancement is measured using the ratios of α-element-sensitive indices, in this case Mgb, to Fe = (Fe 5270 + Fe 5335 )/2. α elements are released into the interstellar medium by Type II supernovae, which start occurring shortly after the onset of star formation, while Fe originates in Type Ia supernovae, which only start ∼ 1Gyr after the onset of star formation. Consequently, the ratio of these elements can be used as an indicator of the star-formation timescale, such that shorter episodes of star formation will give an α-enhanced stellar population due to the enrichment of magnesium from the SNII, and the αenhancement will begin to drop after SNIa appear due to the dilution of magnesium with iron in the interstellar medium. For example, it has been found that the highest [α/Fe] ratios are achieved in galaxies with the shortest half-mass formation time (<2 Gyr; de La Rosa et al. 2011), and that the [α/Fe] decreases with increasing star-formation timescale out to ∼ 10 Gyr (Wiersma et al. 2009).
To measure the α-enhancement, the Mgb and Fe line strengths measured from the binned spectra were plotted onto the SSP model grids of Thomas et al. (2011, hereafter TMJ models). These models give estimates of the relative metallicity and [α/Fe] ratio for a range of ages of stellar populations, and for each galaxy the model was selected that most closely matched the mean luminosityweighted age for galaxy. Since the TMJ models are based on the MILES stellar library, the spectral resolution is similar to that of the MUSE data, and so no corrections were necessary to account for the difference in resolution. However, with velocity dispersions of between 130−300 km s −1 for each galaxy, the line strengths needed to be corrected for this broadening. The correction was carried out by fitting each spectrum with and using the weights for each template spectrum used in the fit to obtain a best fit spectrum at the resolution of the MILES library, i.e. without including for the velocity dispersion of the data. This model spectrum was then convolved with a series of Gaussians to simulate the broadening for a range of velocity dispersions, and the strengths of each line measured in the same way from these broadened spectra. Finally, the necessary correction was calculated for the measured velocity dispersion from the galaxy spectrum via interpolation. The corrected line strengths were then plotted onto the grids to obtain estimates of the α-enhancement, and the results are given in the right columns of Fig. 1 and 2.
It can be seen that the [α/Fe] ratios are relatively flat with radius across most galaxies, with 2MIG 131 and 2MIG 1546 showing a weak negative gradient with radius. Additionally, variations in the [α/Fe] ratios are detected in the cores of several galaxies, with an increase in the inner 2 (0.5 kpc) of CCC 122 that coincides with the older, more metal-rich stellar populations, and a drop in the inner 2 (0.5 kpc) of CCC 137 where the significantly younger, more metal-rich populations are found. These results hint again at different star-formation histories in the cores of these galaxies compared to their outer regions.
The [α/Fe] ratio for 2MIG 445 is particularly interesting, with a similar shape to the age plot and map within 5 . Together, these results indicate that the light in the very core of this galaxy is dominated by stars created during a short episode of star formation, surrounded by a ring of younger stars that were formed during a more extended period of star formation.
However, the strength of all of these features and their uncertainties are of a similar magnitude to the scatter in the measurements, and little can be said for certain without looking at the properties in each component within the galaxies.

EXTRACTING BULGE AND DISC SPECTRA WITH BUDDI
Clues to how S0s have transformed from their progenitor galaxies and which processes led to the truncation of star formation lie in the stellar populations of their bulges and discs. However, in the inner regions of S0 and spiral galaxies, the light from the disc is always superimposed upon the bulge, and so the radial stellar populations represent the light from both components in varying fractions (e.g. Kennedy et al. 2016). As a result, cleanly extracting the spectra from either of these components is difficult, and if an incorrect model is used, the stellar populations analysis can suffer from contamination of light from other components. One way to cleanly extract the spectra of individual components within galaxies is to create a wavelength-dependent model of the light profile of each component with (Bulge-Disc decomposition of IFU data; Johnston et al. 2017).
uses a modified form of G (Peng et al. 2002(Peng et al. , 2010 called G M Vika et al. 2013) to model the light profile of multi-waveband images of a galaxy simultaneously. User-defined Chebychev polynomials are implemented to account for variations in the structural parameters of each component as a function of wavelength, thus using information from all the images for each fit. As a result, the S/N of the data set is boosted over that of any individual image, thus allowing reasonable estimates for the structural parameters to be made for images with lower S/N. While this technique was initially developed for large scale surveys with a handful of wavelengths/bands, it is particularly powerful for IFU datacubes, which can be considered as a series of narrow-band images of the target at each wavelength step in the spectral direction, as it can produce reliable fits to image slices with low S/N or which are affected by sky lines. While was originally written to model the bulges and discs of galaxies, it has been successfully applied to fits to extreme galaxy components, from extended stellar haloes (Johnston et al. 2018) to the nuclear star clusters within dwarf galaxies (Johnston et al. 2020).
It is important to note that uses only the light profile information of the galaxy within the datacube to create the wavelength-dependent models for each component and thus extract their clean spectra. No information is provided or extracted by regarding the spectra or the line strengths. Consequently, the spectra extracted for each component is independent of any previous assumptions or measurements of the stellar populations, flux fractions or mass fractions of each component. The main bias that can be introduced by the fits is the number of components and their starting parameters, and to combat this issue, users are encouraged to repeat steps 2 and 3 below to ensure the best fit has been achieved.
A full description of can be found in Johnston et al. (2017), with a particularly useful flow chart in Fig. 1 of that paper, and a brief overview is presented below.

Step 0: Create the mask and PSF profile
Before modelling a galaxy with , a few preparations are needed for the data, such as creating a bad pixel mask to mask out pixels within the image that fall outside of the MUSE field-ofview (FOV). Some of the datacubes also contain light from other sources, such as foreground stars or background galaxies. Generally it is preferable to fit these objects simultaneously to better model for the light from the target galaxy. However, in this study it was decided to mask out these objects since including models for even just the brightest objects within the FOV would have significantly slowed down the fits. Furthermore, since the kinematics corrections described in the next step are only applied to the dominant kinematical component, any distortion in the kinematics from these objects would likely lead to contamination of the spectra extracted from the target galaxy.
In order to achieve a good fit to the image of a galaxy, the model must be convolved with a point-spread function (PSF) profile that matches the image quality of the data. Consequently, to achieve a good fit with , a PSF datacube must be provided. Ideally this model datacube would be created from non-saturated stars within the FOV of the science image, but the large size of the majority of the galaxies within the MUSE FOV and the lack of isolated stars in the fields meant that this option was not always possible. In these cases, between 1 and 6 stars were identified in the corresponding sky exposures, which were observed as part of the same observing block as the galaxy exposures and thus display similar image quality to the science data. In rare cases where the sky exposure contained no suitable stars, the standard star field was used. In such cases, the standard star was found to be very bright, with excess light spreading horizontally in the images within the slicers, and so other, fainter stars within the FOV were used to create the PSF datacube, generally between 1 and 4 stars. Since the standard stars were observed at different times during the night and in different parts of the sky, a single standard star exposure was selected for each target datacube that best matched the image quality of the corresponding data, as measured from the FWHM of a point source in the white light image of the target datacube. These model PSF profiles were tested by modelling a star or other point source within the FOV of the science datacube and evaluating the residual images. While this approximation of the PSF is not perfect, it was found to provide a satisfactory fit for each datacube.

G
M can only fit symmetric models to an image of a galaxy, and with the narrow wavebands of image slices in an IFU datacube, one must be careful when modelling image slices at wavelengths of strong spectral features. For example, the rotation of a galaxy can mean that part of the image slice might reflect the flux of the continuum while another part may be in an absorption feature, resulting in an asymmetric image. Additionally, the velocity dispersion of a galaxy may decrease with radius, meaning that spaxels from the centre of the galaxy are from within the wings of a spectral feature while those at larger radii are within the continuum. Both these effects can seriously affect the fits to the affected image slices, resulting in messy spectra. Therefore, the first step carried out by is to measure and smooth out the kinematics across the galaxy in order to reduce these effects.
For this step, the datacubes were binned using the Voronoi tessellation technique of Cappellari & Copin (2003), and the kinematics of the binned spectra were measured using the penalised Pixel Fitting software (pPXF) of Cappellari & Emsellem (2004) using the template spectra of Vazdekis et al. (2010), which are based on the MILES stellar library of Sánchez-Blázquez et al. (2006) and consist of 156 template spectra ranging in metallicity from −1.71 to +0.22, and in age from 1 to 17.78 Gyr. The spectra in each spaxel within the datacube were then shifted to correct for the line-of-sight velocity, and broadened to match the maximum velocity dispersion measured within the galaxy. The kinematics across this sample of galaxies are discussed in detail in Paper I.

Step 2: Fit the white-light image
The first step in the actual light-profile fitting process was to model the white-light image of the galaxy with G M. This step allows the user to quickly assess the fit in terms of the number of components necessary and the starting parameters before moving onto the more time-consuming steps to fit the images throughout the rest of the datacube. Each fit started with a uniform sky background, hereafter referred to as the sky component, and a single Sérsic profile, with additional Sérsic profiles added successively to the model until the best fit was achieved. The best fit was decided by looking at the residual images, which were created by subtracting the best fit model from the input image, to find the balance between including enough model components for a good fit while at the same time using the minimum number of components necessary to create a good fit.

Step 3: Fit the narrow-band images
The datacube was then rebinned into a series of 10 narrow-band images with high S/N spanning the full spectral range of the MUSE datacube, and these images were modelled with G M using the best fit to the white light image as the starting parameters. The first fit was completely unconstrained, simulating the results obtained by running G on each narrow-band image independently and giving an idea of the variations in the parameters as a function of wavelength. The fit was then repeated using Chebychev polynomials to constrain the parameters for effective radius (R e ), Sérsic index (n), axis ratio (q) and position angle (P A) as a function of wavelength while leaving the integrated magnitude (m tot ) free. This step is repeated if necessary until the Chebychev polynomials give a good representation of how each parameter varies with wavelength and the residual images show a good fit at all wavelengths. This step defines the polynomials for all parameters but the magnitudes in the next step of the fitting process. In general, the P A and q were modelled using polynomials of order 1 (flat with wavelength) for all galaxies, and R e and n were modelled with both polynomials of order 1 and 2 (linear variation with wavelength). Examples of these polynomial fits to the parameters from the free fits are given in Fig. 3 for CCC 122 and 2MIG 1546, showing cases where the fits using orders 1 and 2 result in small and large differences in the final spectra. For example, on can see that the two fits to 2MIG 1546 result in up to a ∼ 10% difference in the light fractions and SEDs for the bulge and disc components in the blue end of the spectrum, while these differences are almost negligible in CCC 122. In 6 of the eight galaxies, the fits using these three different orders converged on similar fits, while in 2MIG 131 and 2MIG 1814 the fits using a polynomial of order 1 resulted in very different structural parameters In the case of 2MIG 131, the fits using a polynomial order 1 resulted in a more compact bulge with a lower Sérsic index, and the light fractions for the spectra from each component varied by up to ∼ 40%. On the other hand, the flat fit to 2MIG 1814 converged on a model with a very extended, very faint component which was ∼ 4−5 magnitudes fainter than the disc component in any other fit. Tests were carried out on these two galaxies, modelling their narrowband images with different combinations of initial parameters and polynomial fits, and it was found that these different fits were derived only when all the parameters except the magnitudes were modelled with polynomials of order 1. Consequently, it was decided that this order of polynomial induced too much restriction in the fit, and that the final fits were untrustworthy. Consequently, for these galaxies, the flat fits were carried out with polynomials of order 2 for R e and order 1 for n, q and P A.
Ideally, at this step, one wants to select the lowest order of polynomial that provides the best match to the results from the free fit. In Fig. 3, an order of 1 appears to provide a good fit to CCC 122, while 2MIG 1546 clearly requires a higher order for the disc R e and bulge and disc n. However, in this study, the small number of galaxies meant that the fits to the image slices described in the next step could be carried out using both orders, thus allowing a comparison of the spectra and stellar populations for each component before deciding which order to use for the analysis.

Step 4: Fit the image slices and extract the decomposed spectra
The final step of the process is to find the best fit model for every image slice within the datacube, and to use the results from these fits to derive the clean spectra for each component. To reduce the required computing time for this step, the datacube was split up into batches of 10 consecutive image slices, which were modelled by G M simultaneously, with the parameters for the R e , n, q and P A held fixed as a function of wavelength according to polynomials derived in the previous step. This constraint on the model ensures a reliable fit to the galaxy in each image slice since the S/N of the individual images are lower than that of the stacked images used in the previous two steps, while also preventing artefacts appearing in the fit parameters at wavelengths where each batch of 10 image slices meet. The final parameter included in the fit is the integrated magnitude of each component, which is allowed to vary with complete freedom in each image slice. Consequently, the integrated magnitudes derived reflect total flux of each component at the wavelengths of the image slices used in that fit. By plotting these values for the total flux as a function of wavelength for each component, their clean 1-dimensional spectra are created, which are integrated to infinity. As described above, uses no information on the spectral properties as part of the fit, simply seeing a series of images to model. As a result, the spectra created by are free from bias on the stellar populations or line strengths within the galaxy. Additionally, the combination of the 2-dimensional fit and the polynomials used for the structural parameters lead to consistent fits to each image slice, and thus very little noise in the final spectrum, especially for large, bright galaxies like those used in this study.
An example of these spectra extracted for each component in CCC 137 is given in Fig. 4, showing the fluxes of each component relative to the total flux from the galaxy. These final spectra represent the clean model spectra for all the components included in the fit, and can be used for stellar populations analysis of those structures. This figure also plots the sky spectrum, which was created from the sky component included in the fit, and the residual spectrum, which reflects the remaining flux after subtracting the combined best-fit spectrum from the total flux of the galaxy. The residual spectrum consists of light from faint asymmetric features within the galaxy and background objects within the FOV, which will be discussed in more detail in Section 5.4. Figure 4 also shows a zoom-in of the spectral region covering the Mgb and Fe lines used for the luminosity-weighted stellar populations analysis in Section 6.1, showing the normalised spectra for the bulge, disc and lens of CCC 137 superimposed. It can be seen that the bulge spectrum shows stronger lines, indicating a higher metallicity, followed by the disc and then the lens. Additionally, a difference in the shape of the Mgb triplet can be seen between the bulge and lens spectra due to the different strengths of the Mg 5183.6 line to the right of the triplet.
By this stage, two spectra had been extracted for each component within each galaxy, as derived using the fits with polynomials of order 1 and 2. The analysis described in the following sections was carried out with both sets of spectra, and while small variations were found in the individual measurements, the same trends were seen with both data sets. In the end, the final decision on which set of spectra to use for each galaxy was based on a combination of the better fit to the structural parameters from the free fit (see Fig. 3), fewer artefacts in the residual images (e.g. see Section 5.4), and the stellar populations. For example, where the line strengths fell off the Single Stellar Population (SSP) model grids used in Section 6.1, the order that led to the smallest offset between the line strengths in the spectrum and the nearest point in the grid was selected as the better fit. If these steps failed to identify the better fit, then the simpler model was used for the analysis. In this way, the spectra derived using the flat fits (order 1) for 2MIG 131, 2MIG445, CCC 43 and CCC 122, and the linear fits (order 2) for 2MIG 1546, 2MIG 1814, CCC 137 and CCC 158 were used for the analysis presented in the rest of this paper.

The number and nature of the components in each galaxy
The models used for each galaxy were created with no prior information from fits to photometric data in the literature to reduce bias during the fitting process. Unlike many other studies where galaxy light profiles are modelled (e.g. Simard et al. 2011;Mendel et al. 2014;Vika et al. 2014;Bottrell et al. 2019), the fits were not restricted to single or double Sérsic profiles, where the latter is taken to represent a bulge+disc fit. Instead, the fits were started using a single Sérsic profile, and the complexity of the model was built up by adding additional Sérsic or PSF components until the best fit is obtained. The best fit was determined first by comparing the χ 2 , and then by studying the residual images of each fit to identify the model that gave a good fit to the data with the minimum number of components. For example, one can compare the flux levels in particularly bright or dark regions of the residual image, which represent areas where the light was under and over subtracted by the light of the model, or look at the standard deviation of the all the unmasked flux values in the residual image. Additionally, broad rings of alternating light and dark regions can be a characteristic signature of several issues, such as poor initial estimates for the structural parameters, that the light profile of the galaxy is cored and that a core-Sérsic profile is required instead of a standard Sérsic profile (e.g. Dullo & Graham 2014), or that an additional component that is required for the fit. Finally, when the inclusion of an additional component showed negligible differences or improvements the residual image, it was considered unnecessary and the best fit was accepted as the fit with one less component. It was found that all the galaxies in this sample achieved a better fit using at least three components, with an additional fourth component in the form of a PSF profile in the cases of 2MIG 1814, CCC 122 and CCC 137. The fits to all the galaxies are given in Appendix A, with the residual images from the 2-component fits given in Appendix B for comparison. It can be seen that in the 2-component fit, the fits to the light profile along the major axes are generally quite good, but the residual images reveal more artefacts, such as the rings of alternating light and dark regions outlined above. Thus, these images emphasise the importance of using the two-dimensional spatial information wherever possible. Similar findings were presented by Laurikainen et al. (2009) for their fits to S0 galaxies, where better results were obtained when fitting 2-dimensional images over a one-dimensional light profile in cases where the galaxy contains more than simply a bulge and disc.
All galaxies were found to contain a clear extended disc structure alongside two more compact components with Sérsic light distributions. Similarly complex models for the bulges have been seen before -composite bulges, which are made up of two or more structures such as pseudobulges, central discs or boxy/peanut bulges, have been detected in small samples of galaxies by Erwin et al. (2003), Nowak et al. (2010), de Lorenzo-Cáceres et al. (2012) and Erwin et al. (2015), and a study by Méndez-Abreu et al. (2014) found that 70% of their sample of barred disc galaxies host composite bulges. The images of each component in Appendix A demonstrate the diversity of these inner components, with some showing rounder profiles while others are more elongated, and with a range of Sérsic indices between 0.5 and 4.5. These structures likely represent a mix of classical and pseudobulges, as well as structures such as lenses, bars and thin discs.
Lenses have been detected in photometric studies of S0s by Laurikainen et al. (2005Laurikainen et al. ( , 2009, and are distinguished from bulges by their flatter surface brightness profile and sharper edges. They have been found to have generally lower surface brightnesses than pseudobulges (Kormendy & Kennicutt 2004). Bars and discs on the other hand can be identified through their kinematic signatures. For example, Krajnović et al. (2011) found that the distribution of h 3 against v/σ can be used to distinguish between rapidly rotating galaxies with and without a bar or ring present in an IFU datacube. Furthermore, Guérou et al. (2016) demonstrated that a clear trend in this distribution appears when the spaxels or bins corresponding to a nuclear disc are highlighted.
Having identified the disc in each galaxy, two more compact Sérsic components remained. Within these components, the bulges were selected by eye as the more extended component of the two with a steeper surface brightness profile and the softer edge (i.e. less sharp drop in the flux in the outskirts), and in 6 of the eight galaxies this component was found to be the brightest structure at the core of the galaxy, with the exceptions being CCC 43 and CCC 158. The remaining component was found to be elongated in around half of the galaxies-2MIG 445, 2MIG 1546, CCC 122 and CC158-and relatively circular in the remaining galaxies. This component may represent either a bar, an inner embedded disk or a lens. None of the galaxies show evidence of a strong bar in their visual appearance, and no kinematical signatures were detected in Paper I. In the case of an inner disc, one would expect it to have a similar inclination and position angle as the outer disc, and would expect to find a kinematical signature of a rotating disc within the dispersion-supported bulge. Simulations by Eliche-Moral et al. (2018) of S0 galaxies containing lenses viewed at different inclinations found that the lenses can show a range of shapes when viewed edge-on, with the majority appearing vertically thin and disc like (∼ 64%) while the rest displayed vertically thick lentil and spheroidal shapes. Consequently, the inclination angles were not used in this study, and instead the kinematics were used to distinguish between a lens and an inner disc.
The h 3 Gauss-Hermite polynomial parameter gives a measure of the skewness in the line-of-sight velocity distribution (LOSVD), and is generally found to anti-correlate with the velocity curve in rotating disc galaxies. If another kinematical component is present, its distinct velocity distribution will further affect the shape of the LOSVD of the galaxy as a whole in that region. Therefore, by plotting h 3 against the v/σ, one can look for evidence of additional kinematic components.
Using the kinematics information from Paper I, the plots for h 3 against v/σ are given in Fig. 5 for each galaxy, where the measurements from within the disc, bulge and lens/inner disc regions have been highlighted in blue, red and green respectively. For clarity, the measurements marked as bulges and lenses indicate only those binned spectra where the mean flux of each spaxel in that bin are greater than half of the peak flux of that component in the white-light image of the galaxy. In general, all galaxies show an anti-correlation between h 3 and v/σ, as expected, but the red and green points appear to show distinct trends within the plots for all the cluster galaxies that are similar to that seen for the nuclear disc in NGC 3115 by Guérou et al. (2016). These trends indicate that the kinematics are dominated by a fast rotating disc, but have a significant tail towards lower velocities due to the presence of the slower rotating bulge in the inner regions of the galaxy. If the kinematic The plots for v/σ versus h 3 for each binned spectrum, with those within the disc, bulge and lens components highlighted in blue, red and green respectively. The isolated galaxies are on the left while the cluster galaxies are on the right. signature of an inner disc was present, one would expect to see distinct trends in this figure between the bulge and this disc due to the bulge being larger and higher above the plane of the disc. Since no clear differences can be seen in the kinematics of the bulges and lenses/inner discs in these plots, we can only say that if an inner disc is present, its kinematic signatures are masked by those of the bulge and outer disc. The fact that these trends are only seen in the cluster galaxies may reflect a dependence on environment, but could also be a result of these galaxies being more rotationally supported, making such trends easier to see.
At this point it is interesting to note the unusual relationship between h 3 and v/σ in 2MIG 1546, which shows a correlation at larger radii and an anti-correlation at lower radii. This trend indicates that in the outskirts of the galaxy an additional kinematic component is present that is rotating faster than the main disc component but 2MIG131 2MIG445 2MIG1546 2MIG1814 CCC43 CCC122 CCC137 CCC158 Figure 6. White-light images of each galaxy alongside the unsharp masked images. The galaxies on the left are the isolated galaxies while those on the right are the cluster galaxies. never dominates the light in those regions. Such a phenomenon would result in a tail towards higher velocities in the LOSVD. In the inner regions, on the other hand, the tail in the LOSVD switches to lower velocities, which reflects the contribution of the slowerrotating bulge kinematics to the LOSVD. The velocity and velocity dispersion maps for this galaxy in Paper I appear relatively smooth, though a distortion in the h 3 and h 4 parameters can be seen to the west of the centre of the galaxy along the major axis. While not discussed in that paper, this feature could reflect this additional kinematic component, which may be an artefact of a minor merger in the history of the galaxy Another technique used by Laurikainen et al. (2005) to distinguish between lenses and inner discs is to apply unsharp masking to an image of the target, which can reveal structures such as spiral arms and edge-on discs. This technique was thus applied to the white-light images of each galaxy using a Gaussian with a kernel of 4 pixels. These images are given in Fig. 6, and in CCC 158 it is possible to see evidence of an edge-on disc in the core of the galaxy. Further evidence for a nuclear disc would be a drop in the velocity dispersion and an anticorrelation in the h 3 moment within that region. The kinematics maps for each galaxy are presented in the Appendix of Paper I, and it is possible to see both these signatures in the maps for CCC 158. These kinematic signatures correspond approximately with the light profile of the third component in this galaxy, which is also brighter than the bulge in the core of the galaxy, thus indicating that the trend seen in Fig. 5 may instead reflect the kinematics of this component as opposed to the bulge. If this component is an inner disc, it may contain distinct stellar populations relative to the bulge and outer disc. Revisiting the stellar population maps in Fig. 1 show particularly young and metal-rich stellar populations in this same region, though in this case the ellipse of younger stellar populations has a length of ∼ 15 or 3.6 kpc while the effective radius of the third component is ∼ 10 . Since younger stars tend to dominate the light from a galaxy, it is possible that this third component appears more extended in the stellar population maps because it contains younger stars than the other two components. Consequently, it is possible that this component in CCC 158 may be an embedded inner disc seen at high inclination, while in all the other galaxies in the sample this component is likely to be a lens. However, for the rest of this paper, we will refer to this component in all galaxies as a lens unless discussing CCC 158 explicitly.

Structural parameters of each component
Having identified the discs, bulges and lenses, it is now possible to look at their structural parameters in more detail. In recent years, bulges have been separated into two categories: classical and pseudobulges, where classical bulges are typically dispersion-supported spheroids with higher Sérsic indices while pseudobulges are rotationally supported, discy and with lower Sérsic indices. Kormendy & Kennicutt (2004) use the information for the bulge vs disc ellipticity and the kinematics of the galaxy in the form v/σ to distinguish between these classifications, while Fisher & Drory (2008) define a Sérsic index of n = 2 to separate the morphologies, and Fabricius et al. (2012) use a combination of the kinematics and the properties derived through bulge-disc decomposition. Since no standard system is common throughout the literature, we have chosen to adopt the Fisher & Drory (2008) technique in this work. Figure 7 plots the measurements for the R e and n for each bulge and lens component against that of the corresponding disc. Where the R e was modelled with a Chebychev polynomial of order 2, the value of R e was taken from the fit to the narrow-band image closest to the centre of the MUSE spectral range. It can be seen that 6 of the galaxies have bulges with n > 2.0, classifying them as classical bulges while the other two, CCC 43 and CCC 158, are pseudobulges. It is interesting to note that these two galaxies are the only cases in which the lenses dominate the light at the core of the galaxy, which is inconsistent with the findings of Kormendy & Kennicutt (2004) that lenses tend to have lower surface brightnesses than pseudobulges. Using the kinematics information derived in Paper I, it was also found that the central velocity dispersion in each galaxy correlated with the bulge morphology, such that the classical bulges have σ 0 > 200 km s −1 while the psuedobulges have σ 0 < 200 km s −1 .
Furthermore, all the lenses and inner discs have n < 2.0, and with the exception of CCC 43 and CCC 158, the bulges all have higher Sérsic indices than the lenses within the same galaxy. This trend is likely to be the result of the identification of the bulge components having the steeper light-profile, characteristic to higher values for the Sérsic index (see Section 5.1). In the case of CCC 158, both the bulge and the lens have similar Sérsic indices, whereas in CCC 43 the bulge has a significantly lower Sérsic index. In both cases however, the bulge was selected to be the more extended and less elliptical component.
The discs all have Sérsic indices n 2.5, which is consistent with the fits to disc galaxies in the Coma Cluster by Weinzirl et al. in which they found that galaxy discs are not universally exponential. Another feature that is worth mentioning in Fig. 7 is that the discs of the isolated galaxies have higher Sérsic indices than those of the cluster galaxies. A grey dashed line has been added to this plot at n = 1.35 to highlight this result. This trend may indicate that the structure of discs is different in field and cluster S0s, once the contamination from other components has been properly removed. However, with the small sample used in this study, we cannot explore this phenomenon further, nor rule out target selection bias.

Colours
The polynomials used to fit the structural parameters for the fits to the narrow band images, as described in Section 4.4, provide information on how those parameters vary as a function of wavelength. From this information, the colours of each component can be derived, which, in turn, can be used to compare their stellar populations.
The absolute AB magnitudes of each component are plotted against wavelength from the fits to the narrow band images in Fig. 8, and clearly show that all the components included in the models to these galaxies are brighter at redder wavelengths. It can be seen that the absolute magnitudes of the discs (top panel) of the cluster galaxies are generally fainter than those of the isolated galaxies. This effect may simply be a selection bias due to the small sample of galaxies and the limited range in their properties. Paper I showed that the isolated galaxies in this sample have slightly higher masses in general than the cluster galaxies, with mean masses of log 10 (M/M ) ∼ 11.4 as opposed to ∼ 10.9, which may explain their brighter discs. The magnitudes of the bulges (middle) and lenses (bottom), on the other hand, show no real correlation with environment.
The luminosity fractions of each component relative to the total integrated luminosity of the galaxies are also plotted in the middle panel of Fig. 8, and show that in all galaxies, the lenses are fainter than their corresponding bulges and discs, with luminosity fractions of 0.2. The light fractions of the lenses are generally consistent with the light fractions of bars in S0s in the Coma Cluster, where Lansbury et al. (2014) found light fractions of 30%. Interestingly, the lenses within the cluster galaxies appear to show higher luminosity fractions than those in the isolated S0s, particularly at redder wavelengths, though again, this trend may be a selection effect. Additionally, the bulges and discs in many galaxies have relatively flat luminosity fractions with wavelength. This trend may indicate that the bulges and discs have similar SEDs, which could result from poor fits where the bulge and disc have been improperly fitted or where a single component has been modelled with two Sérsic profiles. In such cases, one would expect to see similar structural parameters for both components. However, upon checking the parameters and the images of each component included in the model (see Appendix A), this scenario could be ruled out, and it is therefore more likely that both the bulge and disc in these galaxies have similar colours.
In two of the the isolated galaxies, the bulges are more luminous in the red while the discs are correspondingly fainter in the red and brighter in the blue. These trends may indicate more recent star formation in the discs or higher metallicities in the bulges of those galaxies. Consequently, they may hint at a dependence on the environment or the total mass of the galaxy, but with the small number of galaxies within this sample we cannot rule out that this effect is simply a coincidence.

Residual features
Another benefit of using to model these galaxies is that a residuals datacube is created by subtracting the best fit model created by G M from each image slice. While in each individual image slice the residuals are almost negligible, the white-light image created from this datacube reveals faint structures that would otherwise be hidden in the data. These white-light images created from the residual datacubes are shown for each galaxy in Fig. 9, having been scaled to highlight these features, and as part of Appendix A, using the same flux scaling as the white-light image for the same galaxy to allow a better comparison of the relative luminosity of these structures. These images reveal faint structures along the line of sight of each galaxy, including brighter globular clusters within the discs and background galaxies, and thus show the potential to use to cleanly model the light from foreground objects and cleanly extract the spectra of more distant objects within the field.
One can also see that all galaxies within this sample show evidence of faint features within the galaxies themselves, such as dust lanes or remnants of spiral arms, the latter of which are not visible in the unsharp-mask images in Fig. 6 nor in the datacubes and white-light images. These spiral arms are particularly prominent in both the inner and outer regions of CCC 122. A comparison with the components used to model this galaxy in Appendix A shows that the inner spiral arms lie within the component identified as the lens. While no clear spiral arms were detected in the unsharpmask image of this galaxy in Fig. 6, the inner region of this galaxy does appear to have distinct kinematics in Fig. 5. Consequently, it is still possible that the lens component in this galaxy is actually an inner disc, based on the definition used by Laurikainen et al. (2005). CCC 158 on the other hand also shows evidence of distinct kinematics in the lens component and a possible inner edge-on disc in the unsharp-mask image, but no clear spiral structure is detected 10" 2MIG131 2MIG445 2MIG1546 2MIG1814 CCC43 CCC122 CCC137 CCC158 Figure 9. The residual white-light images for the isolated (left) and cluster (right) galaxies, showing the remnants of tidal tails, dusty discs, and spiral arms. The images have been scaled to best display these structures.
in the residual image for this galaxy. However, the residual image for this galaxy does contain a lot of structure, which may be created by dust extinction or spiral arms in the outer disc seen edge-on, which may act to hide spiral arms within the lens if present.
The faint spiral arms that can be seen in some of the residual images, particularly for the cluster galaxies, could be evidence of their transformation from spiral galaxies and provide a link between passive spirals and S0s (e.g. Fraser-McKelvie et al. 2018a;Pak et al. 2019). The residual features in the isolated galaxies, on the other hand, appear more asymmetric, which may reflect different formation processes in each environment. These findings are consistent with those of Paper 1, where the cluster galaxies are more rotationally supported and are thought to have formed through gas removal, while the isolated galaxies are more dispersion supported, which is thought to reflect a transformation more dominated by minor mergers which would disrupt the spiral structure. Attempts were made to try to extract the spectra from the brighter features in the residual datacubes in order to study the different stellar populations across the arms and inter-arm regions, but the faintness of these structures and the consequent low S/N in the extracted spectra made the subsequent analysis unreliable.
Of particular interest is the residual white-light image for 2MIG 445, where the faint features appear to be remnants of shells or tidal tails following a merger. No clear remnant of the infalling galaxy can be seen, suggesting that if these features were created during a merger, it occurred long enough ago that the merging galaxy has been disrupted completely, while sufficiently recently that the tidal tails have not yet dissipated entirely (∼ 1-2 Gyr ago, cf. Eliche-Moral et al. 2018). Alternatively, the merger could have been minor, thus preventing significant disruption to the structure of the galaxy. While the stellar population maps for this galaxy show no clear trends that could be interpreted as the differences between those of the progenitor and infalling galaxies, the kinematics maps in Fig. A1 of Paper I show a clear distortion to the North-West of the centre of the galaxy, particularly in the velocity dispersion and h 4 parameters. No clear optical counterpart can be seen in Figures 6  and 9, thus indicating that this disturbance in the kinematics may be a remnant of a merger where the infalling galaxy has been destroyed and the kinematics progenitor galaxy are still somewhat disrupted.

Luminosity-Weighted Stellar Populations
Having extracted the spectra for the bulges, discs and lenses within each galaxy, estimates for their luminosity-weighted stellar populations were derived using the same techniques outlined in Sections 3.1 and 3.2. The results are given in Fig. 10, showing the correlations between the ages, metallicities and α-enhancements between the discs, bulges and lenses, and Fig. 11, which presents the trends between these properties for each component separately. The spectra for the PSF components were found to have very low S/N and their line strength measurements fell significantly off the grids in the direction of old, metal-rich stellar populations. As a result, they have been omitted from the stellar populations analysis presented in this work. It can be seen that with the exception of 2MIG 445 and CCC 158 all of the discs are generally old, with ages > 8 Gyr, while the bulges are all older than ∼ 6 Gyr and the lenses show a wide range of ages. Due to the old ages of the components, only large differences in the ages (i.e. several Gyr) between each component would be detected in these plots. Consequently with the large uncertainties in these age measurements, no clear trend can be seen between the ages of the discs with either those of the bulges or lenses or as a function of environment. This result is consistent with the relatively flat age gradients and generally old ages for all the galaxies in Figures 1 and 2. However, the lenses of 2MIG 445, CCC 43 and CCC 158 do show evidence of significantly younger stellar populations, indicating that these components hosted the most recent star formation within these galaxies.
The results for the metallicities are more interesting, showing that the metallicities of the lenses are generally consistent with or higher than the discs, with the exception of CCC 137, which may reflect that the lenses have formed in situ from the disc material. Furthermore, with the exception of CCC 158, the bulges generally display consistent or higher metallicities than the discs as well, though no clear trend is seen between the metallicities of the lenses and bulges. These higher metallicities in the more compact components compared to the discs is similar to the findings of Johnston et al. (2012Johnston et al. ( , 2014 that the bulges of S0s in the Virgo and Fornax Clusters generally contain more metal-rich stellar populations, though the younger stellar populations in the same regions is not replicated in this work. This difference in results may be due to the increased information available in IFU spectroscopy compared to long-slit spectroscopy, or due to the additional components that can be more reliably modelled with this data set. Figures 1 and 2 also showed higher metallicities in the inner regions of all the galaxies, and so it now appears likely that these gradients are a result of the superposition of varying fractions of light from the bulges, discs and lenses at each point in the galaxy. For example, the younger, more metal-rich stellar populations detected at small radii CCC 43, CCC 158 and 2MIG 445 in Figures 1 and 2 can be attributed to the lens component in Fig. 10. Similarly, the uniformly old stellar populations measured at all radii in 2MIG 1814 and CCC 122 in Figures 1 and 2 are also reflected in the ages of the individual components. Together, these plots show the strength of at extracting stellar populations parameters for individual structures within these galaxies despite using no information on the stellar populations in the fits. In the plots for the α-enhancement, it can be seen that the [α/Fe] ratio for the bulges are correlated with those for the discs, indicating that the most recent star formation activity in both components were connected. Only 2MIG 1546 differs, with a higher α-enhancement in the bulge than the disc. A similar trend between the bulge and disc α-enhancement was seen in Johnston et al. (2014) for S0s in Virgo, in which it was attributed to a final episode of star formation within the bulge, fuelled by gas that originated in the disc, which took place after the star formation in the disc had ceased. However, in that study the bulges were found to contain younger, more metal-rich stars than their discs, which is not the case in this sample. Instead, almost all the bulges and discs show evidence of old stellar populations, with the differences in their ages being smaller or similar to the uncertainties in the estimates. Consequently, the consistent [α/Fe] ratios between these two components suggests that they also formed over similar timescales.
The [α/Fe] ratios of the lens components on the other hand show no clear correlation with either the bulges or discs, indicating that they were formed through a distinct and more random star-formation episode. The metallicities of the lenses however are generally consistent with or higher than the discs, with the exception of CCC 137, which may reflect that the lenses have formed from the disc material. This scenario therefore rules out the theory that the lens components are simply artefacts resulting from a poor bulge+disc fit to S0 galaxies.

Mass-Weighted Stellar Populations
The stellar populations derived through analysis of line strengths only provides the luminosity-weighted properties, which can be dominated by stars produced in a recent star formation event and therefore may not reflect the true age of an old galaxy. Consequently, they are best used to determine how long ago the most recent star formation activity occurred. However, through full spectral fitting using template spectra for a wide range of ages and metallicities, it is possible to obtain estimates of the mass-weighted stellar popula-  Figure 11. A comparison of the luminosity-weighted ages, metallicities and α-enhancement between the discs (left), bulges (middle) and lenses (right) of each galaxy tions and star-formation histories for each component within these galaxies to better understand their relationship.
The mass-weighted ages and metallicities were derived from the spectrum of each component using to apply a regularized fit to the spectrum. A multiplicative Legendre polynomial of order 10 was also included in the fit to correct the shape of the continuum, reducing the sensitivity to dust reddening and omitting the requirement of a reddening curve (Cappellari 2017). The fit used a linear combination of template spectra of known relative ages and metallicities from the MILES library spanning an age and metallicity range of 0.06 − 18 Gyr and [M/H]=−2.32 to +0.22 respectively. The regularization process within works to create a smoothed variation in the weights of the templates with similar ages and metallicities, where the degree of smoothing is controlled by the regularization parameter determined by the user. The process started with an unregularized fit to each spectrum to measure the χ 2 value for the fit, and then the noise spectrum was scaled appropriately until χ 2 /N DOF = 1, where N DOF is the number of degrees of freedom in the fit, which corresponds to the number of unmasked pixels in the input spectrum. The fit to the spectrum was then repeated using this scaled noise spectrum and increasing values for the regularization parameter until the χ 2 of the fit increased by ∆ χ 2 = √ 2 × N DOF . This value represents the limit between a smooth fit that still reflects the star-formation history of the galaxy and one that has been smoothed excessively. It should be noted at this point that this smoothed fit may not reflect the true starformation history of the galaxy, which is likely to vary over shorter timescales than the models allow, but instead acts to reduce the age-metallicity degeneracy between spectra, thus allowing a more consistent comparison of systematic trends in their star-formation histories.
Since the template spectra are modelled for an initial birth cloud mass of 1 Solar mass, the final weight of each template reflects the 'zero-age' mass-to-light ratio of that stellar population. Consequently, the smoothed variation in the weights of neighbouring template spectra in the final regularised fit represents a simplified star-formation history of that part of the galaxy by defining the relative mass contribution of each stellar population (see e.g. McDermid et al. 2015). These star-formation histories for each component are given in Fig. 12, showing both the mass fraction associated to each step in age and metallicity within the models, and the mass fraction created as a function of lookback time. It is immediately clear in these plots that all the components are old and metal-rich, with only the disc and lens of 2MIG 445 and the bulge and disc of CCC 158 showing evidence of significant recent star formation. These results are in agreement with the luminosity-weighted stellar populations, which show evidence of younger stellar populations within these same components. However, these plots for the star-formation histories go beyond to show that all of the galaxies in this sample have underlying old, metal-rich stellar populations, indicating that the majority of their mass formed long ago and that any younger stellar populations detected represent a later star-formation episode that contributed very little towards the mass of that component.
The mean values for the mass-weighted ages and metallicities of each component can be calculated from the weights of the template spectra using and respectively, where ω i represents the weight of the i th template (i.e. the value by which the i th template stellar template is multiplied to best fit the galaxy spectrum), and [M/H] template,i and Age template,i are the metallicity and age of the i th template respectively. These results are plotted in Fig. 13, and the associated uncertainties were calculated through a series of simulations in which random noise added to the best fit model to each spectrum to replicate the S/N of the original spectrum. One hundred simulated spectra were created for each component in this way, and the mass-weighted ages and metallicities were measured from these spectra using the same values for the regularization. The final uncertainties were calculated as the standard deviation in the measurements from these simulations. In general, these results also reflect that the majority of the mass within these structures is made up of of old, metal-rich stellar populations, with only the lens and disc of 2MIG 445 showing evidence of a significant fraction of young stars. It is possible to see in Fig. 13 that the mass-weighted metallicities of the bulges and lenses of the isolated galaxies show much smaller distribution than the same components in the cluster galaxies, whereas the opposite is true for the disc components. Additionally, the bulges of all the isolated galaxies and CCC137 again show higher metallicities than the discs, while the bulge of 2MIG 445 is significantly older than the disc. These stellar populations could also explain the redder bulges and bluer discs in isolated galaxies compared to the cluster galaxies in Fig. 8. However, the low numbers of galaxies from each environment in this study makes it difficult to determine if this trend is true or coincidence, and it will be followed up in a future work with a larger sample of galaxies.

DISCUSSION AND CONCLUSIONS
We have carried out a detailed analysis of the stellar populations and star-formation histories of the independent components within a sample of unbarred S0s galaxies in isolated and cluster environments. Not only were the models of these galaxies expanded beyond the typical assumptions of an exponential disc surrounding a Sérsic or de Vaucouleurs bulge, but we also incorporated the first spec-troscopic analysis of lenses within S0s, based on our search of the literature. Through analysis of the spectra for each component, we derived estimates of their stellar populations in order to determine their contributions towards the evolution of the galaxies as a whole. Using , the light profiles of each galaxy were modelled in each image slice of the MUSE observations in order to cleanly model the luminosity and extract the spectrum of each component, allowing measurements of their stellar populations with minimal contamination.
The light profile fits revealed that all eight galaxies required more complex models than a simple bulge+disc fit, consisting of three extended components in all cases with a further three galaxies requiring a PSF profile at their cores as well. In each case, the disc was identified by eye as the most extended component, while the bulge was distinguished as the inner spheroidal component with the steeper surface brightness profile and softer edge, which in the majority of the sample also corresponded to the brightest component at the core of the galaxy. In each galaxy, the third extended component satisfied the criteria for being a lens, having a flat surface brightness profile with a sharp outer edge (Kormendy 1979), but they did display a range of elongations, with axes ratios ranging from 0.41 to 0.96. The kinematics, unsharp mask images and stellar populations maps of each galaxy indicated that this component may instead be an inner edge-on disc in CCC 158, but without conclusive evidence it has been included it under the lens umbrella for the purpose of this work.
The stellar populations within each component were analysed through both line strength measurements and regularized full spectral fitting to obtain estimates of the luminosity and mass-weighted properties respectively. It was found that in general, the bulges and discs of all galaxies are very old, with only the discs of 2MIG 445 and CCC 158 showing significantly younger stellar populations. This result is inconsistent with the findings of Johnston et al. (2012Johnston et al. ( , 2014, where the bulges of S0s in Virgo and Fornax were found to have younger and more metal-rich stellar populations than their discs. This difference could be a result of the different evolutionary pathways of the samples of galaxies used, but may also be an artefact from the fits to the one-dimensional light profile in the long-slit spectra in that work and the inclusion of the lens component in this work. Laurikainen et al. (2009) found that modelling the twodimensional light profile of a galaxy gave superior results than using the one dimensional light profile, particularly when the galaxy contains an additional component beyond the bulge and disc. However, the older bulges detected in this study are consistent with other studies of S0s. For example, Fraser-McKelvie et al. (2018b) found that the bulges in MaNGA S0s with similar masses to those used in this study are generally older than their discs, and Katkov et al. (2014Katkov et al. ( , 2015 identified that S0s in isolated environments tend to have similar ages in both the bulges and discs. Méndez-Abreu et al. (2019b) also found that the bulges in S0s are older than their discs, which they interpreted as evidence that star formation occurs throughout the disc, and only in the disc, not in the bulge. Another study by Rodríguez Del Pino et al. (2014) of cluster 'k+A' galaxies, which are thought to be an intermediate phase in the transformation of spirals to S0s, found that the final episode of star formation was located within the extended disc, though more centrally concentrated than the earlier star formation.
The grids showing the mass fractions created as a function of time and metallicity further expand on this result, showing that nearly all the bulges and discs have an underlying high mass of intermediate to old stars, with some more recent star formation occurring that dominates the light in these cases but only contributes towards a small fraction of the mass of that component. In general, it can be seen that the majority of the mass in the discs was built up over a longer timescale than in the bulges, which is consistent with spiral galaxies, which have star-forming discs and red bulges, being the progenitors of these S0s. However, in CCC 122 and CCC 137 the star-formation timescales in the bulges and discs appear to be similar, possibly indicating that the star formation was truncated soon after the progenitor spiral was formed. With the exception of CCC 158, the luminosity and mass-weighted metallicities of the bulges are all consistent with or higher than those of the discs, indicating that they formed from the same material. Furthermore, the [α/Fe] of the bulges and discs show a strong correlation, indicating that their most recent episodes of star formation are somehow connected. Morelli et al. (2008) and Morelli et al. (2016) studied the stellar populations within bulges of cluster and isolated S0s respectively, finding that galaxies in both environments display null age and negative metallicity gradients with radius, and that the isolated galaxies show negative [α/Fe] gradients while the cluster galaxies show no slope in this property. With these results they concluded that the bulges all formed through gas dissipating inwards and fuelling star formation in the inner regions of the galaxy, with successive mergers only playing a significant role in the cluster galaxies, where they have flattened out the gradients in [α/Fe]. On the other hand, Katkov et al. (2015) proposed that the similar [α/Fe] ratios in the bulges and discs instead point to a scenario where star formation in both components started simultaneously and ceased at around the same time. Furthermore, Fraser-McKelvie et al. (2018b) identified that the process that truncated the star formation is more dependent on the galaxy mass than environment, though their sample of galaxies didn't include S0s from the densest environments, and are best used to compare to the results for the isolated galaxies in this work. Our results for the bulges and discs are consistent with these studies. The similarly old mass-weighted ages and similar or higher metallicites in the bulges relative to the discs present a scenario where the majority of the mass in the bulges and discs was created long ago, and that these components have survived without experiencing any major mergers that sufficiently disrupted their discy morphology since then (Prochaska Chamberlain et al. 2011). The wider range in luminosity-weighted ages however indicate that most galaxies have experience a small amount of more recent star formation after the majority of their mass was built up, with the consistent α-enhancements between these components hinting that this most recent episode of star formation occurred over similar timescales in the bulge and disc of each galaxy. However, the youngest components and more extended star-formation histories are observed mainly in the field S0 of our sample (see Fig. 12), suggesting that minor events (e.g. minor mergers) might have occurred after the formation of the main bulk of stars in less dense environments; this scenario is consistent with the conclusions of Paper I, in which the kinematic properties of field S0s suggest that they have experienced minor mergers. Furthermore, the residual images of 2MIG 445 revealed shells and tidal tails, which would further support this theory.
The lenses were found to have generally higher metallicities than the discs, indicating that they may have formed in situ from the same material. However, where the bulges and discs show a strong correlation in the [α/Fe] ratios, the lenses show no correlation with either of these components. Furthermore, the lenses in CCC 43 and CCC 158 showed younger luminosity-weighted stellar populations, while the mass-weighted ages of all the lenses were found to be the same as or older than the discs. The lens of 2MIG 445, however, shows a mass-weighted age of ∼ 2 Gyr and thus has been formed more recently. Simulations by Eliche-Moral et al. (2012,2018) found that the typical timescale for a galaxy to become a relaxed elliptical or S0 following both minor and major mergers is 1−2 Gyr, though Borlaff et al. (2014) found a timescale of 3.5 Gyr to create S0s with truncated discs. These timescales are comparable to the age measurement for the lens in 2MIG 445, and thus are consistent with the scenario where this galaxy has undergone a recent merger that may have contributed towards the formation of the lens. Morelli et al. (2008Morelli et al. ( , 2016 proposed that the flat radial gradients in [α/Fe] ratios are produced by successive mergers mixing the stellar populations. Therefore, if true, the lack of correlation in the [α/Fe] ratios in the lenses with those of the bulges and discs could mean that the most recent star formation in the lenses occurred after the mergers, while an alternative scenario is that the most recent episode of star formation in the bulges and discs occurred over similar timescales and were fuelled by the same material, while the lenses formed independently over more random timescales within the lifetimes of their host galaxies. Another open question about the origin of lenses is their connection to bars. Kormendy (1979) proposed that lenses in unbarred galaxies are the remnants of dissolved bars that have spilled out into a lens shape. A study by Laurikainen et al. (2013) found that the stellar populations derived from photometry of bars and lenses in S0s are similar, which they believe is consistent with this theory that lenses are partially destroyed bars. On the other hand, models by Athanassoula et al. (1983) showed that galaxies that are too dynamically hot to form bars may still be cool enough to form lenses through disc instabilities, while Bosma et al. (1983) found evidence that lenses in unbarred galaxies are redder than their discs, and so likely formed through truncated star formation early in the life of the galaxy. If lenses were connected to bars, one would expect to see similar stellar populations properties in bars and lenses relative to the disc. In the literature, bars have been found to be more metal-rich and less α-enhanced than their surrounding discs by Neumann et al. (2020), while Carrillo et al. (2020) found that the bar in NGC 2903 has similar metallicity and a higher fraction of younger stars than the outer disc, which they interpreted as evidence that gas flowing along the bar has triggered more recent star-formation activity. Further evidence for ongoing star-formation along the leading edges of bars in spiral galaxies has also been found by Peterken et al. (2019) andFraser-McKelvie et al. (2020). In the sample of unbarred S0s used in this work, we found that the lenses tend to have higher metallicities than discs, while the discs display more recent or longer episodes of star formation than the lenses. While these results are consistent with the theory that lenses are partially dissolved bars, they cannot rule out the scenarios proposed by Athanassoula et al. (1983) and Bosma et al. (1983).
Together, the results presented in this work present a scenario in which the bulges, discs and lenses in these S0 galaxies formed in situ from the same material, but while the most recent star formation activity in the bulges and discs were truncated at around the same time and after similar star formation timescales, the lenses appear to have formed at more random times over the lifetime of the galaxy and over different timescales. One explanation for these results could be that the bulges and discs have formed through dissipational processes, such as major mergers or gas accretion, at high redshift, with later evolution through accretion of less massive satellites and internal mechanisms induced by a bar. The presence of younger stellar populations, asymmetric features in the residual images and the dispersion supported kinematics reported in Paper I in the field S0s may indicate that these galaxies have been affected more by minor mergers than the cluster galaxies, but the low numbers used in this study mean that no firm conclusions can be derived as to the factors affecting the formation of the bulges, discs and lenses as a function of environment. However, this scenario is consistent with the theories for the formation of high-mass S0s by van den Bergh (2009), Barway et al. (2013) and Méndez-Abreu et al. (2019b). The lack of correlation in the ages and α-enhancements of the lenses with the bulges and discs may indicate a later formation scenario for these components from evolved bars, which fits with the scenario that the frequency of bars drops in S0s while they also appear more evolved in these galaxies (Laurikainen et al. 2009).

APPENDIX A: SELECTED LIGHT PROFILE FITS
This section gives an overview of the fits to each galaxy in the sample, showing the white-light image, the best fit model and the residual image from left to right in the top row, and the models for the disc, bulge and lens component in the bottom row. Where a PSF component was included in the fit, it is not shown in these images. All images for the same galaxy are plotted to the same flux scale to allow a direct comparison of the strength of any features visible in the residual images. In general, the fits appear relatively good, with the residual light representing mainly faint asymmetric features within the galaxy such as remnant spiral arms or dust lanes within the disc.
On the right of these figures is the one-dimensional light profile fit for that galaxy. The black points show the light profile from the white-light image from a slit of width 1 along the major axis, and the blue, red, green and light blue lines show the light profile fits for the disc, bulge, lens and PSF components (where applicable). The purple line shows the combination of these components, and the grey dashed line gives the sky background level. The darker grey box represents the seeing FWHM, which is quantified on the bottom left of each plot.

APPENDIX B: 2-COMPONENT LIGHT PROFILE FITS
This section presents the fit to each galaxy using a simple double Sérsic profile to represent only the bulge and disc components. In the top row are the white-light and residual images, and below the one-dimensional light profile is plotted in the same way as in Appendix A. All images have been scaled to the same greyscale as in Appendix A to allow for a direct comparison between the 2 and 3 or 4 component models.. One can see that the fits to the one-dimensional light profiles are relatively good in some cases, but in general the residual images clearly show more artefacts. In particular, one can see a bright core (indicating under subtraction of the model in that region) surrounded by alternating dark and bright rings, which are a characteristic signature that an additional component is required to improve the fit. These results are reminiscent of Laurikainen et al. (2009), who found that better results were obtained when fitting 2dimensional images of their sample of S0s over a one-dimensional light profile in cases where the galaxy contains more than simply a bulge and disc. This paper has been typeset from a T E X/L A T E X file prepared by the author.  Figure A1. The fits to 2MIG 131 and 2MIG 445, showing in the top row the white-light image, the best fit model and the residual image, and below are the models for bulge, disc and lens components. The images of the PSF components are not included. On the right of these figures is the 1-dimensional light profile along the major axis of each galaxy, along with the light profiles of each component included in the model and the best fit model. The sky background is given by the horizontal dashed line, and the grey box represents the seeing FWHM measured form the datacube, with the numerical value at the bottom-left of the plot. The additional noise and sharp cut-off at high radii indicate the measurements from close to the edge of the FOV, where the S/N is lower due to fewer exposures covering that region as a result of the dithering and rotations between exposures.   Sky Data Figure B2. As for Fig. B1, but for the Cluster galaxy sample.