Formation of S0s in extreme environments I: clues from kinematics and stellar populations

Despite numerous efforts, it is still unclear whether lenticular galaxies (S0s) evolve from spirals whose star formation was suppressed, or formed trough mergers or disk instabilities. In this paper we present a pilot study of 21 S0 galaxies in extreme environments (field and cluster), and compare their spatially-resolved kinematics and global stellar populations. Our aim is to identify whether there are different mechanisms that form S0s in different environments. Our results show that the kinematics of S0 galaxies in field and cluster are, indeed, different. Lenticulars in the cluster are more rotationally supported, suggesting that they are formed through processes that involve the rapid consumption or removal of gas (e.g. starvation, ram pressure stripping). In contrast, S0s in the field are more pressure supported, suggesting that minor mergers served mostly to shape their kinematic properties. These results are independent of total mass, luminosity, or disk-to-bulge ratio. On the other hand, the mass-weighted age, metallicity, and star formation time-scale of the galaxies correlate more with mass than with environment, in agreement with known relations from previous work such as the one between mass and metallicity. Overall, our results re-enforce the idea that there are multiple mechanisms that produce S0s, and that both mass $and$ environment play key roles. A larger sample is highly desirable to confirm or refute the results and the interpretation of this pilot study.


INTRODUCTION
Lenticular galaxies (S0s) outnumber galaxies of other morphological types in the local Universe (Bernardi et al. 2010), and their number density has increased with time in clusters and groups since z ∼ 1 at the expenses of spirals (e.g., ⋆ E-mail: lcoccato@eso.org Dressler 1980;Poggianti et al. 2009). Therefore, understanding their formation is a key element in understanding galaxy evolution.
S0s have long been thought of as quenched spirals since they share the disky morphology of spirals, but have the redder colours of older stellar populations. Evidence for this evolutionary scenario comes from studies such as Dressler (1980); Dressler et al. (1997); Cappellari et al. (2011b), which show that the fraction of S0s increases towards higher density environments and lower redshift, while spirals show the opposite trend. A popular explanation is that spiral galaxies falling into clusters lose their gas via ram-pressure stripping (e.g., Gunn & Gott 1972;Steinhauser et al. 2012 or tidal interactions (e.g., Merluzzi et al. 2016), thus enhancing the population of lenticulars. Quenching of star formation in spiral galaxies is also a possible mechanism: among spiral galaxies hosting a classical bulge, the fraction of quenched galaxies increases with denser environments, which then lead to an increase in the fraction of S0s (Mishra et al. 2018). However, observations of nearby isolated S0s (van den Berg et al. 2009) have suggested that there must be alternative evolutionary pathways to form S0s. One possibility is that isolated S0 galaxies are the end product of past minor mergers that used up the gas in the disc, concentrating it to the central parts of the galaxy to form a classical bulge, analogous to "fossil groups" (Ponman et al. 1994;Arnold et al. 2011). Similarly, mergers could induce the formation of a bar in the disc that can, in turn, build up a pseudo-bulge. N-body simulations (Bournaud et al. 2005;Eliche-Moral et al. 2012) have indeed suggested that dry intermediate (mass ratios of 1:4−1:7) and minor (<1:7) mergers can induce global structural evolution resulting in S0 systems. Alternatively, S0s could simply be faded field spiral galaxies that exhausted their gas reservoirs and lost their spiral structures through disc instabilities (starvation, Eliche-Moral et al. 2013). Finally, low mass S0s could also have originated from primordial galaxies formed at redshift ∼ 2 through violent disk instability and fragmentation (Saha & Cortesi 2018).
Much observational effort has been made to distinguish between these scenarios. For instance studies, of nearby S0s (either isolated or in groups) have revealed that their discs are dynamically hotter than those in spirals of similar luminosity (Cortesi et al. 2013), and that their blue globular cluster population has a wide range of ages (Chies-Santos et al. 2011;Lee et al. 2019), suggesting past interaction with other galaxies, and therefore favouring minor-merger scenarios for the formation of S0s. In contrast, other recent studies on globular cluster kinematics (Bellstedt et al. 2017) do not favour mergers as formation mechanism for lenticular galaxies, although they cannot entirely rule them out. Moreover, using 3D spectroscopy, Katkov et al. (2014) found that the number of field S0s with counter-rotating gas kinematics is higher than in denser environments, implying that this gas could have been accreted from dwarf satellites. However, the S0 Tully-Fisher relation (Tully & Fisher 1977) has been found to be systematically fainter at fixed rotational velocity than the spiral relation (Williams et al. 2010;Cortesi et al. 2013;Jaffé et al. 2014), suggesting that S0s are faded spirals that consumed their gas reservoirs. Other studies have further shown that the luminosity of bulges in S0s, relative to their associated discs, are brighter than expected from a simple cessation of star formation in the disc (Christlein & Zabludoff 2004). In addition, Johnston et al. (2012Johnston et al. ( , 2014 found evidence of younger stellar populations in the bulge regions of Fornax and Virgo S0s, triggered by residual disc gas that was channelled into the centre of the galaxy, inducing star formation and thus increasing the luminosity of the bulge relative to the disc. Similar behaviour is observed in 13 post-starburst spiral galaxies in the Abell Cluster S1077 (AC114), which had their last episode of star formation towards the central parts of the galaxies (Rodríguez Del Pino et al. 2014), suggesting indeed a link between S0s and spirals.
Broadly speaking, the main formation mechanisms of lenticular galaxies discussed above related either to mass accretions (e.g. merger) or gas-related effects (e.g. ram pressure stripping, gas accretion, quenching of star formation). Undoubtedly, different environments can contribute in different ways to these mechanisms (see for example Aguerri 2012 and Vollmer 2013 for reviews). It is therefore fair to ask whether or not S0s in the field have different properties from S0s in clusters. The best way to characterise the global properties of galaxies is to combine morphological and kinematic information, with a spatially resolved analysis of their stellar kinematics and populations out to large radii. To this end, we have embarked on a project aimed at studying the properties of lenticular galaxies in extreme environments. Recent attempts (e.g. Fraser-McKelvie et al. 2018;Rizzo et al. 2018) found no significant dependency with the environment; in particular Fraser-McKelvie et al. (2018) indicates mass as the main driver for S0 formation. However their sample did not contain galaxies in cluster, whereas in this study, we seek to optimise any signal by select disk galaxies from the extremes of environments, field and cluster. The main scientific drivers of our study are: i) investigate whether or not S0s living in extreme environments have different properties. Then, if differences are indeed present, ii) establish a link between these differences and the formation mechanisms, and hence iii) identify the dominant formation mechanism that acts in a given environment. This paper presents the first results of the project, based on tailored pilot MUSE (The Multi Unit Spectroscopic Explorer Bacon et al. 2010) observations of a small sample of galaxies and complementary integral-field data from the literature. Here, we concentrate mainly on the spatiallyresolved kinematic analysis (stellar v/σ and specific angular momentum), presence of ionised gas, and on the average properties of the stellar population (mass-weighted age, metallicity and star-formation time-scale). Future works will present a detailed analysis of the stellar populations of the galaxy structural components (such as bulge and disk). The sample of galaxies studied in this paper is described in Section 2. Data reduction and analysis of new galaxies are discussed in Section 2.1. The measurement of global properties both for new and literature data and the results of their comparison is discussed in Section 3. Finally, Section 4 presents the initial conclusions from this study.

THE SAMPLE OF GALAXIES
Galaxies in our sample were chosen from a pilot MUSE program that includes dedicated observations of 8 S0s in extreme environments (i.e., field and cluster) 1 , plus 13 galaxies from publicly-available data of the ATLAS3D survey (Cappellari et al. 2011a) to augment the sample.
One part of the sample comprises isolated galaxies, drawn from the 2MASS Isolated Galaxies Catalogue (2MIG, Figure 1. Distribution of the basic properties of the sample galaxies. The plots show the total mass versus the effective radius (upper plot) and the total luminosity (lower plot), the histogram of their distribution on the sides of the axis. Blue indicates field galaxies, whereas red indicates cluster galaxies. The size of the symbols is proportional to the Disk-to-total light ratio, disky galaxies are represented by bigger symbols. Galaxies in the MUSE sample are highlighted by a black circle. Karachentseva et al. 2010). For the MUSE campaign we selected 4 2MIG targets that were visually classified as S0s by the team, and that were observable from Paranal. To these were added the 8 S0 galaxies studied by the ATLAS 3D survey that are present in the 2MIG catalogue or that are in an environment with galaxy three-dimensional density 2 lower than log ρ 10 /[M pc −3 ] ≤ −2.5. The rest of the sample galaxies belong to dense galaxy clusters. For the MUSE campaign we selected 4 S0s in the Centaurus cluster, from the sample of Jerjen & Dressler (1997). Tho these were added the 5 galaxies in the ATLAS 3D survey that live in an environment with a galaxy threedimensional density greater than log ρ 10 /[M pc −3 ] ≥ 1 3 .
All the galaxies in this study have morphological Hubble type −3.5 < T < −.5, from LEDA and do not show clear evidence of strong bars or tidal interactions in their visual appearance. These morphological criteria were the same as those used in the ATLAS 3D survey to disentangle S0s from ellipticals. However, because the classification of S0s in the past was very much influenced by the limited dynamic range of photographic plates, we have also visually inspected the sample galaxies. In addition, for the MUSE sample, the presence of a disk-like rotating stellar structure is also reflected in the kinematics (see Appendix A). For the ATLAS 3D sample, their S0 nature was also confirmed by Krajnović et al. (2013), with the exception of 8 galaxies (marked with asterisks in the last column of Table 2). The morphological structure of these 8 galaxies was uncertain to Krajnović et al. (2013) and therefore these galaxies were classified as "single Sersic component objects" (although some of them, e.g. NGC 6548, are clearly dominated by a disk component). We therefore performed an independent image decomposition on archival images (VLT-VIMOS or 2MASS archives) to verify the presence of a disk component in the photometric profile 4 . The inclusion or exclusion of these 8 objects and their impact on our results will be discussed in Section 4.1.1.
The total number of lenticular galaxies considered in this study is 21, 12 isolated and 9 in cluster (although not the densest environments in the Universe, the Virgo and Centaurus Clusters from which these objects were selected are very much cluster environments, and have the key benefit of observability and plentiful data from the literature). Table 1 and Figure 1 summarise their main properties. One can immediately see from Fig. 1 that the sample in this pilot study is not homogeneous: massive and luminous galaxies are under-represented in the cluster sample. The effects of selection biases will be discussed in Sect. 4.1. The scatter on the relations shown in Fig. 1 is mainly dominated by measurement errors (uncertainty on galaxy distance has been included in the luminosity error computation) and by the degeneracy between baryonic matter and dark matter that is present in the dynamical models.
In the following we describe the MUSE observations and data reduction.

The MUSE observations and data reduction
Observations were executed with the Multi Unit Spectroscopic Explorer (MUSE) mounted on Unit Telescope 4 a the ESO La Silla Paranal observatory (Chile). MUSE was configured in WFM-NOAO-N mode (wide field of view, no galaxies inside a sphere centred on the galaxy and containing the 10 nearest neighbours. 3 The adopted thresholds on ρ 10 (1 and -2.5) represent the densest and less dense environment bins of the ATLAS3D survey (e.g. Figure 8 in Cappellari et al. 2011b) For galaxies in the ATLAS 3D sample, magnitudes and distances are as reported in Cappellari et al. (2011a). Error on magnitude includes the contribution from the error on the distance. Column 7: Total apparent "face-on" magnitude corrected for galactic and internal extinction, and for redshift (from de Vaucouleurs et al. 1991). Error on magnitude includes the contribution from the error on the distance. Column 8: WISE (Wright et al. 2010) W1 magnitudes at 3.4 µ, corrected for internal and galactic extinction, and with aperture and k corrections (Sorce et al. 2012). Error on magnitude includes the contribution from the error on the distance. Column 9: log 10 of the effective radius. For galaxies in the MUSE sample, the value is fitted with galfit (Peng et al. 2002) on the reconstructed image obtained from the datacubes its associated error is ∼ 10%. For the ATLAS 3D sample the value of R e is as reported by Cappellari et al. (2013), its associated error is ≈ 22%, according to Section 4.1 of Cappellari et al. (2013). Columns 10 and 11: values of V rot /σ and the λ parameter within 1 effective radius, as computed in Sections 3.1 and 3.2. Column 12: flag indicating the presence of ionised-gas in the spectra (y=yes, n=no, i=only in the innermost regions). Column 13 (for MUSE sample only): observational set-up, indicating the exposure time in seconds for each exposure, the number of exposures on target (O) and on sky offset (S, if acquired).
adaptive optics, nominal wavelength coverage) that ensured a spatial sampling of 0.2 arcsec per spaxel and a spectral coverage of 4750 -9350Å with a nominal resolving power of R ∼ 2000.
Observations were executed in service mode between October 2015 and February 2016 (Period 96) and organised in a series of observing blocks. Each observing block contains several exposures on target that were slightly dithered and rotated by 90 degrees with respect to each other to minimise the signature of the IFU geometry in the final reduced product. For the targets that had a size comparable to the field of view, dedicated offset sky exposures were included in each observing block.
Data reduction was performed with the MUSE pipeline (Weilbacher et al. 2012) version 2.2 executed under the Es-oReflex environment (Freudling et al. 2013). Depending on the projected size of each target, the contribution of the sky background was computed either from dedicated exposures executed close in time to the science observations (within the same observing block) or from the edges of the field of view, where the contribution of the galaxy was negligible. After removing the sky background, the exposures were then aligned and co-added using bright sources as reference. Whenever available, the cubes of dedicated sky exposures were also sky subtracted and co-added following the same dither pattern as the science observations, creating "master cubes" of sky residuals. Sky residuals were further reduced using the ZAP algorithm (Soto et al. 2016); this method models the sky residuals with a series of principal components and then it fits and subtracts them from the science cubes. The skyresidual components were evaluated on the master sky residual cubes or on the border of the field of view, depending on the dimension of the target.
Adjacent spaxels in the datacube were co-added using Voronoi tessellation as implemented by Cappellari & Copin (2003), seeking a target signal-to-noise ratio of 50 per pixel. Spectra in each spatial bin were fitted using the ppxf (Cappellari & Emsellem 2004) and gandalf (Sarzi et al. 2006) spectral fitting procedures to extract stellar and ionized-gas kinematics. Stellar templates from the MIUS-CAT spectral library (Vazdekis et al. 2012) at the nominal MUSE spectral resolution of FW H M = 2.51Å were used in the fit. The spectral range running from 4750-7000Å was used, with regions affected by high residuals from sky or telluric lines masked and excluded from the fit. The two- Table 2. Best fit parameters of the kinemetry and Jeans axisimmetric models.

ANALYSIS
In this section we define the global galaxy properties we use to study the differences between cluster and field S0s and the analysis carried out to measure them. In particular, we will focus our attention on V rot /σ and the proxy for specific angular momentum λ. Then, we compute the Tully-Fisher relation between the circular velocity and the total luminosity in K-band, B-band and W1 (3.4 µm) band. Finally, we derived the mass-weighted values of age, metallicity and star formation time-scale within 1R e .

V rot /σ radial profiles
The V rot /σ radial profile is computed as the ratio of the stellar rotation velocity and velocity dispersion calculated along the kinematic semi-major axis. We performed a harmonic first-term expansion of the observed two-dimensional velocity and velocity dispersion maps using the kinemetry code by Krajnović et al. (2006). At each radius a we fit the amplitude of rotation V rot (a) (corrected for inclination), the kinematic position angle P A(a), and the projected axial ratio q(a). Errors are directly computed by the Levenberg-Marquardt least-squares minimisation algorithm used by the kinemetry code. The systemic velocity (V syst ) is assumed constant with radius. For each galaxy we then define the median kinematic position angle P A KIN and axial ratio q . Their errors are computed as standard deviation of the computed values, divided by the square root of the number of radial bins. Values for each galaxy are reported in Table 2. The velocity dispersion profile σ(a) is obtained as weighted average of the measured velocity dispersion computed on concentric ellipses with constant position angle and ellipticity (from P A KIN and q ), as determined by the fit to the velocity two-dimensional field. Figure 2 shows the radial profiles of V rot /σ as function of distance along the semi-major axis expressed in units of effective radii, and the distribution of V rot /σ interpolated at fiducial radii. The figure shows differences in the distributions of V rot /σ radial profiles of the two families of galaxies; in particular, there is indication that lenticulars in the field reach lower values of V rot /σ. If we use 1R e as reference radius, the mean V rot /σ for field galaxies is: 1.2 ± 0.2, whereas for cluster galaxies is: 1.8 ± 0.2. The mean values differ by ∼ 3σ; if we consider the errorbars, the null-hypothesis that the difference of mean values is consistent with 0 is discarded at ∼ 1.5σ level (i.e the errobars touch each other if multiplied by ∼ 1.5). If we use 1.5R e as reference radius, the separation of V rot /σ is more evident, but error bars get slightly larger because the number of galaxies for which we obtain Figure 2. Comparison between the V rot /σ radial profiles of field (blue) and cluster (red) galaxies in our sample. Thick lines identify the MUSE sample, thin lines the ATLAS3D sample. The right-hand side of the figure shows the histograms of the value of V rot /σ at 1 R e and 1.5R e . Figure 3. Comparison between the radial profiles of the proxy for specific angular momentum (λ) in field galaxies (blue) and cluster galaxies (red) in our sample. Thick lines identify the MUSE sample, thin lines the ATLAS3D sample. The right-hand side of the figure shows the histograms of the value of λ at 1 R e and 1.5R e . The green symbol shows the mean value for spiral galaxies at 1R e (diamonds: Graham et al. 2018, triangles: Falcón-Barroso et al. 2019, the error bar is the standard deviation of the values. The pink symbol shows the mean value of ellipticals (diamonds: galaxies in Table B1  kinematics out to that radius is smaller. While the V rot /σ of cluster lenticular galaxies increases from 1 R e to 1.5 R e , it decreases for field lenticulars. The difference in the V rot /σ profiles of cluster and field S0s suggests that cluster S0s are more rotationally supported than field ones, which could in turn indicate different formation mechanisms (see discussion in Section 4).

λ radial profiles
We compare the radial profile of the λ parameter, which is a proxy for the specific angular momentum as defined in Emsellem et al. (2007). The λ profile is computed as function of semi-major axis a of concentric ellipses that have the same flattening and orientation as those used for the computation of V/σ in Section 3.1. Note that λ is a cumulative and luminosity-weighted quantity, whereas V/σ is a "local" kinematic measurement. Figure 3 shows the λ(a/R e ) profiles for the sample galaxies and the distribution of values at fiducial radii. The distributions of λ profiles of the two families of galaxies look different, suggesting that lenticulars in the field tend to reach lower values of λ. Indeed, if we use 1R e as reference radius, the mean λ of field galaxies is 0.36 ± 0.04, whereas it is 0.46 ± 0.03 for cluster galaxies. The mean values differ by ∼ 3σ; if we consider the errorbars, the null-hypothesis that the difference of mean values is consistent with 0 is discarded at ∼ 1.43σ level (i.e the errorbars touch each other if multiplied by ∼ 1.43). If we consider 1.5R e as reference radius, the separation between field and cluster is ∼ 1σ (i.e., the error bars of the two measurements touch). Both field and cluster population show an increase of λ if computed at 1R e or 1.5R e . As in Section 3.1, the different radial profiles of the angular momentum parameter λ found in this section indicate that the field S0s tend to be less dominated by rotation than their cluster counterparts. From inspection of the values of λ(R e ) available in the literature for ellipticals and spirals, we note that our field S0s have specific angular momentum closer to the typical values of ellipticals. By contrast, our cluster S0s have angular momenta closer to, although not as high as, the typical values of spirals (see Fig. 3). This suggests that field S0s have formation processes more simi-lar to those of elliptical galaxies (e.g., Bournaud et al. 2005), and that cluster S0s have formation processes more similar to those of spirals (e.g., Eliche-Moral et al. 2013).

Tully-Fisher relation
The Tully-Fisher relation between the circular velocity and the luminosity of disk galaxies is known to be a very tight correlation for spirals, especially at near-infrared wavelengths (e.g., Sorce et al. 2013). It has been successfully used not only as an extra-galactic distance indicator(e.g., Tully et al. 1992;Courtois et al. 2011), but also to study morphological transformation of disk galaxies (e.g., Williams et al. 2010;Cortesi et al. 2013) and environmental effects on such transformations (e.g., Jaffé et al. 2011Jaffé et al. , 2014. In this section we derive the Tully-Fisher relation of field and cluster S0 galaxies in our sample, and compare it to that of other S0s and spirals from the literature. To construct the Tully-Fisher relation we first had to derive the circular velocity in our galaxy sample. For this we constructed axisymmetric models by fitting the second moment Jeans equations (e.g., Binney & Tremaine 2008) to the data. To do that, we adopted the procedures set out in the Multi-Gaussian Expansion (MGE, Cappellari 2002) and Jeans Axisymmetric Model (JAM, Cappellari 2008) packages. Details of the procedure are outlined in Section 3.3.1. Only those galaxies for which we reached a relatively flat regime in the circular velocity curve are included in the analysis (see end of Section 3.3.1). The error on V circ accounts for errors in the JAM best fit parameter plus a 10% uncertainty due to differences in various methods of evaluating V circ (via Jeans or Schwarzschild models, asymmetric drift, or CO rotation curves), as advocated by Leung et al. (2018). Figure 4 shows the Tully-Fisher relations obtained for field (blue) and cluster (red) S0 galaxies in K, B, and 3.4 µm bands. B-band is more sensitive to star formation, while M k and W1 are better tracers for the underlying older stellar populations. Overall, there is a good agreement within error bars with the relation determined for S0s by other authors (Williams et al. 2010). Lenticular galaxies lie below the Tully-Fisher relation found for spiral galaxies at all bands: at similar V circ , S0s are less luminous that their spiral counterparts. In the figure, we fit the data by fixing the slope of the Tully-Fisher relation to that determined for spiral galaxies. The scatter on V circ , which is mostly dominated by the uncertainties in the model, does not allow us to determine any significant difference in the relations obtained for field and cluster environments. Also, the measured scatter is different at different bands; however not all the same galaxies are used in different bands, therefore a direct comparison of the scatter is not straightforward. A number of outliers are observed in the B-band relation (NGC 4425 and NGC 4706 are over-luminous with respect to the rest of S0s of similar V circ ) and the 3.4 µm-band relation (NGC 6548 and NGC 4429 are under-luminous with respect to the rest of the S0 of similar V circ ). Uncertainties in photometry and distance alone do not explain these outliers; the most plausible explanation for the discrepancies is an additional source of error in their catalogue magnitudes. These outliers are excluded from the fit. , and 3.6 µm (lower panel) for field (blue) and cluster (red) galaxies. The thick blue and red lines show the linear fit to the data; shaded area indicate the 1-σ error of the fit. The slope of the S0 Tully-Fisher relation is fixed to the one determined for the Spirals, whose reference is given in each panel. MUSE data are identified by a black circle. As comparison, the Tully-Fisher relations for spirals (green) and S0s (cyan) from Williams et al. (2010) and Neill et al. (2014) are shown. Outliers are identified with crosses in the plots and removed from the fit.

Calculation of circular velocity
We constructed axisimmetric two-dimensional maps of V rms = √ V 2 + σ 2 by folding observed kinematic quantities. For the calculation of V rms we repeated the ppxf analysis by fitting only velocity and velocity dispersion, not the higher moments, in order to align our analysis of the MUSE observations with the one for the ATLAS3D data. The twodimensional maps were rotated to align the kinematic major axis (determined in 3.1 and listed in Table 2) along the x-axis. The fit to the folded V rms maps was performed using the JAM software following the same methodology adopted in Williams et al. (2009), which we summarise here.
Two different JAM models are considered. One in which the potential of the dark matter follows the light distribution (model A), and one in which we the dark matter has a spherical distribution (model B). In both models, the mass density distribution is parametrised by a sum of two-dimensional Gaussians aligned along the galaxy photometric major axis. Each Gaussian has a given surface density (in M ⊙ pc −2 ), dispersion along the major axis (in arcsec) and axial ratio.
In model A, light traces the total mass; we applied the MGE procedure on the reconstructed images obtained by integrating the MUSE datacubes over the SDSS-r bandpass; these images are deeper than other images available in the literature for the same band. Foreground and background sources were masked to avoid biasing the fit. The light profiles of non-saturated foreground stars were used for the determination of the point-spread functions. For the galaxies in the ATLAS 3D sample, we used the MGE parameters as determined by Scott et al. (2013) in the same band. In model B, we added the contribution of a dark matter halo, which was assumed spherical and following the Navarro et al. (1997) profile. We refer to (Williams et al. 2009) for further details on the equations used.
The total number of free parameters in the fit is five: stellar mass-to-light ratio (M/L, which represents the total mass-to-light ratio in model A, including both luminous and dark matter), inclination (Incl), orbital anisotropy (β), total dark matter content (M DM , untied from the light distribution, which is set to 0 for model A) and scale radius r s (not considered in model A). The influence of a central supermassive black hole is negligible given the spatial resolution of the data.
For the majority of the sample galaxies, the data quality was not sufficient to separate the degeneracy between stellar and dark matter distribution. We therefore adopted best model A for the galaxies in which the M DM and r s are consistent with 0, and best model B for the remaining systems. For the ATLAS3D sample, in which model A has been chosen, we fixed the inclination as in Table 1 of Cappellari et al. (2013), which lists values for the same model A.
Once the best-fit models parameters are determined, we computed the radial profile of the circular velocity using the mge_circular_velocity tool in the JAM software distribution.
We report in Table 2 the best-fit parameters related to the chosen model (A or B), and the "flat" value of V circ , which we then use to construct the Tully-Fisher relation. We classify a circular velocity curve as "flat" if the variation in the last half part of the curve is less than 15%. We then compute < V circ > as the average of the values of V circ (R) in the outermost 75% of the radial range.

Mass-weighted stellar populations
We inferred the mass-weighted ages and metallicities ([Z/H]) of the stellar population in the sample galaxies within 1R e . For the ATLAS3D sample we used the ages and metallicities determined by McDermid et al. (2015), which were obtained by integrating the spectra within 1R e and fitting them with stellar templates from the MILES library (Vazdekis et al. 2012). The mass-weighted values of age and metallicity were determined as weighted sum of the age and metallicity of the templates used (Equations 1 and 2 in their paper). The weights were determined by the fitting routine ppxf executed with regularization. Emission lines where fitted and removed from the spectra using gandalf (Sarzi et al. 2006). Errors were computed by means of Monte Carlo simulations. For the MUSE sample, we adopted the same strategy as in McDermid et al. (2015), and limited our spectral interval to theirs (4750 − 5400Å) to avoid systematic differences in the analysis between the two datasets. From the mass weights of the individual templates, we determined also t 50 , the time in Gyr needed to form 50% of the current-day stellar mass within 1R e , following the prescriptions of Section 5.1 in McDermid et al. (2015). The quantity t 50 offers a useful proxy for the star formation time-scale.
The inferred mass-weighted values of age, [Z/H], and t 50 are listed in Table 3 and shown in Figure 5. There is no apparent difference between the properties of field and cluster S0s. However, there is a tendency for more massive galaxies to have higher ages and metallicities, and lower t 50 . The stronger dependency of stellar population parameters on mass suggested by Figure 5 is not new, and is particularly evident in larger samples from the literature. McDermid et al. (2015) (their Figure 13) showed that the mass-weighted stellar populations in early-type galaxies strongly depend on mass, whereas the correlation with galaxy density is weak though present (galaxies in less dense environments are younger, more metal poor, and have a longer star formation time-scale). In Fraser-McKelvie et al. (2018) only the dependency of stellar populations with mass is evident (their Figure 5), and no dependency on environment is seen (their Figure 9, although they show only the results for the bulgedominated region).
In a forthcoming paper (Johnston et al, in preparation), we will investigate the two-dimensional maps of the stellar populations, exploiting the full wavelength range of the MUSE sample, and separating the contribution of young and old components, and of the disk and the bulge by means of spectral decomposition analysis (e.g., Coccato et al. 2011Coccato et al. , 2018Johnston et al. 2012Johnston et al. , 2017, to determine if there are any subtler dependencies on environment.

DISCUSSION AND CONCLUSIONS
We have investigated the kinematics and stellar population properties in a sample of 21 lenticular galaxies that reside in two extremely different environments (field and cluster) with the aim of detecting any dependence on environment in their formation and evolution. In particular, we compared the radial profile of the ratio between stellar rotation and velocity dispersion (V r ot /σ), and the λ parameter (a proxy for the specific angular momentum), the M k , B and 3.4µm Tully-Fisher relations, and the integrated massweighted age, metallicity, and half-mass formation time scale t 50 at 1R e in 12 field and 9 cluster S0s. The kinematic analysis revealed (at a ∼ 1.5σ level) that S0s in the field tend to have lower V/σ(R e ) and λ(R e ) than their cluster counter parts. In the Tully-Fisher relations however we do not find any quantifiable difference with environment.

Effects of observational biases
Before embarking on the interpretation of these results, we investigate the effects of observational biases in our sample selection.

Morphological mis-classification
The morphological distinction between elliptical and lenticular galaxies could be ambiguous, especially if the object is nearly face-on. As discussed in Section 2, 8 galaxies in the ATLAS 3D sample, classified as S0 according to LEDA, were indicated as "pure" Sersic objects by Krajnović et al. (2013). They are marked with an asterisk in the last column of Table 2. However, our photometric decomposition of their 2MASS images revealed the presence of a prominent disk component in all cases. In order to test our results within 1R e . Colours indicate the field (blue) and cluster (red) galaxies; the size of the symbol is proportional to the total mass indicated in Table 2. Circles identify galaxies from the MUSE sample.
against this possible contamination, we repeated the kinematic analysis after removing these 8 galaxies from the sample. We found that this test does not substantially change the results: the largest distinction between the kinematics is still observed when splitting the sample into field and cluster, although the significance moderately decreases due to the smaller sample statistics (see Figures 7 and 8). Regarding the Tully-Fisher relation, the removal of these 8 galaxies would leave in only 2 galaxies with a flat V circ profile in the field sample, so a sensible comparison with cluster galaxies is not possible. We also explored the extreme and unlikely case of 100% contamination, in which all the sample galaxies are either spirals or ellipticals. From the measured values of λ(< R e ) of spirals and ellipticals in Graham et al. (2018) and Emsellem et al. (2011), we created 1000 mock random catalogues with 12 randomly selected galaxies to represent the field population and 9 randomly selected galaxies to represent the cluster population. We made sure not to have duplicates in each mock catalogue, and that all the mock catalogues are unique. We selected only galaxies with λ(< R e ) < 0.7 as in our observations. We then computed the mean λ(< R e ) and its error for each of the field and cluster mock catalogues and check the significance in the difference between λ of field and cluster galaxies. We found that in only 10% of the simulations the difference between the average values of λ FIELD and λ CLUSTER is at least as significant as our results (i.e., at 1.43σ level, in the sense that the 1σ error bars associated to λ FIELD and λ CLUSTER touch each other if multiplied by 1.43). This fraction is small, but not insignificant; therefore a larger sample is desirable to confirm the results of this pilot study.

Other parameters
As evident from Figure 1, there is not a fully homogeneous mass distribution, luminosity, or bulge-to-total light ratio between field and cluster galaxies. In addition, galaxy properties may depend on inclination and distance from the cluster centre, which would require a much larger sample to marginalize out. Therefore, we might be seeing the effects of the contamination from other parameters rather than the environment.
For example, a mass dependency would be indeed consistent with the recent results of Fraser-McKelvie et al. (2018), in which they clearly show a bimodal dependency of stellar population properties with mass. Their results indicate that mass is the main driver for the formation of lenticular galaxies. According to their interpretation, massive lenticulars mostly form by morphological or inside-out quenching (stellar mass higher than 6 × 10 10 M ⊙ ), whereas less massive lenticulars mostly form by the fading of a progenitor spiral galaxy.
In order to evaluate possible selection biases in our results, we repeated the analysis of Sections 3.2 and 3.1 splitting the sample into massive versus less massive (with the division placed at 6 × 10 10 M⊙), into bulge dominated vs disk dominated (split at a disk-to-total luminosity ratio D/T = 0.5), and into luminous versus faint galaxies (either side of the median value of K-band magnitudes M k = −23.83). Results are shown in Figure 6. It is evident that the kinematic properties of the galaxies do not differ when split in these ways. Differences in the distribution of radial profiles and their values at 1R e are less significant than when the sample is divided by environment.
Thus, this kinematic analysis indicates that mass is not the main driver for different formation scenarios; this may seem at a first glance to be at odds with the recent findings by Fraser-McKelvie et al. (2018). However, one should note that we cover a significantly different range of masses and environments in this work. Our sample galaxies have masses higher than log(M/M⊙) = 10.2, so all belong to the upper mass range of Fraser-McKelvie et al. (2018), where mass-driven differences between their properties are expected to be less evident. Contrastingly, by construction our sample spans a much wider range of environments than Fraser-McKelvie et al. (2018), so should be more sensitive to any dependency there. Moreover, their study was mainly focused on stellar population, so did not explore the kinematic differences uncovered here.
Interestingly, if we consider the results of the stellar population analysis, we reach a somewhat different conclusion. The size of the symbols in Fig. 5 are proportional to the log of the total mass. There does seem to be at least some correlation with mass between points in Fig. 5, but not with environment. The massive galaxies tend to concentrate in the region of higher [Z/H] and age, and low t 50 , consistent with the findings of McDermid et al. (2015) and Fraser-McKelvie et al. (2018), and also as expected from the broader mass-metallicity relation in galaxies (e.g., Gallazzi et al. 2006).
Although our preliminary results suggest that environment shapes the kinematic properties of lenticular galaxies, whereas mass shapes the stellar populations, one would not expect the effects of mass and environment to be totally independent of each other. Indeed, merging and stripping processes can change the total mass of a galaxy and its morphology. Moreover, the time-scale of star formation (which is known to correlate with mass; e.g. Thomas et al. 2005;McDermid et al. 2015) or quenching mechanisms (that can be related to the environment; e.g., Davies et al. 2019) also play a role, and suggest that the reality must be a more complex interplay between mass and environment in determining the global properties of a galaxy.

Formation scenarios
A number of scenarios have been proposed to explain the formation of S0s. As mentioned in the introduction, one can group them into two main categories: the first category involves galaxy interaction (e.g., mergers, tidal interactions); the second involves changes in the gas content (e.g., rampressure stripping, starvation).
Interactions in general tend to destroy the stellar disk leading to less rotationally-supported systems or to an elliptical galaxy. We might therefore expect that S0s formed through this channel would have significantly lower V rot /σ and λ than their spiral counterparts of similar mass (e.g., Bournaud et al. 2005;Querejeta et al. 2015;Bekki & Couch 2011;Lagos et al. 2017). Also, merger remnants would be systematically shifted to lower V circ for the same mass right after their formation, but they would then move to the local Tully-Fisher relation after 4-7 Gyr of passive evolution. On the other hand, if a large amount of gas is involved in the merger, then new stars in rapidly-rotating structures can boost V rot /σ. However, simulations shows that such enhancement affects only the very central regions in the majority of cases (e.g., Lagos et al. 2018). In general, the modification of the gas content has less effect on the dynamical structure of the galaxy, so the resulting S0s would be expected to have the dynamical structure of their spiral progenitors.
Our results point in the direction that the field environment is dominated by the merger processes that decrease V rot /σ and λ. The fact that we do not observe systematic differences between the Tully-Fisher relations of S0s in different environments (at least within the scatter of our measurements) suggests that such mergers occurred at high redshift, and then the merger remnant has evolved passively and in isolation for 4-7 Gyrs (Tapia et al. 2017) moving back to the local Tully-Fisher S0s relation. On the other hand, dense environments seem to be dominated by the stripping processes that change V rot /σ and λ less dramatically. We clarify that here we consider a process to be "dominating" only with respect to the other mechanisms that contribute to form S0s, not with respect to all the mechanisms that Figure 6. Comparison between the values of profiles of V /σ (upper panel) and the proxy for specific angular momentum λ (lower panel) at 1 effective radius of the sample galaxies. Galaxies are divided into massive and less massive (left panels), bulge dominated and disk dominated (central panels), and luminous and faint (right panels). The threshold values used to separate the sample are: total mass = 6 × 10 10 M ⊙ , bulge-to-total ratio = 0.5, M k =-23.68 mag, i.e. the median of the M k magnitudes. Same results are obtained using W1= -20.48, i.e. the median of the W1 magnitudes.
can occur in a given environment and that can change, for example, the morphological type of an S0 galaxy.
From a qualitative point of view, the kinematic results are consistent with the idea that i) isolated S0 galaxies are the end product of past mergers resulting in a stellar component "dynamically hotter" than cluster S0s, and ii) cluster S0s are formed through processes that involve the rapid consumption or removal of gas resulting in a stellar component "dynamically colder" than field S0s. S0s in cluster could therefore be "faded" spirals, although a small contribution from minor mergers cannot be rule out (at least in some cases) to explain why the mean λ(R e ) of cluster S0 is not as high as those of spirals. Also the effect of the cluster environment is more directed to the gas content (e.g. ram pressure stripping), consistent with the fact that the majority of our cluster S0s have have even less gas in their central regions than S0s in the field (see Table 1).

Summary and future perspective
We have studied the properties of a sample of lenticular galaxies in the field and in clusters, with the aim of investigating the influence of the environment on their formation processes. Even with the relatively small number of galaxies in our sample, we find clear indications that galaxies in the field are generally more pressure supported than galaxies in the cluster, suggesting that galaxy interactions (mergers, tidal interactions) play a more major role in the field in shaping thee systems. This scenario is also reflected in the fact that the specific angular momentum of field S0s is closer to the range of values of elliptical galaxies, which are believed to form mainly via mergers. On the other hand, the kinematic properties of S0s in the cluster are closer to their spiral counterparts, indicating that processes involving the modification of the gas content (stripping, starvation, rapid star formation) are dominant in the cluster environment .
Mass or luminosity do not seem to be the main drivers of these kinematic properties; however in our study we have only considered galaxies with total mass log(M/M ⊙ ) > 10.2, so we do not cover the faint end of the S0s population. By contrast, the properties of the stellar populations (age, metallicity, t 50 ) are more correlated with mass than with environment, consistent with previous studies that note the correlation between, for example, mass and metallicity (e.g., Thomas et al. 2005).
One of the main cautionary aspects of this pilot study is the sample size, which is not entirely representative of the S0 galaxy population nor is it uniform in parameter distribution. A larger sample is therefore highly recommended to obtain unequivocal conclusions that either confirm or negate the findings and interpretation presented here. Notwithstanding this limitation, this work highlights that environment does seem to be an important parameter in shaping the kinematic properties of lenticular galaxies. While the previous study by Fraser-McKelvie et al. (2018) indicated stellar mass as the main driver for dictating the formation mechanism of S0s, a direct comparison between their work and this analysis is difficult to make at this stage because of the differences in mass range end environment explored, Figure 7. As Figs. 2 (upper panels) and 3 (lower panels) but after having removed the 8 galaxies marked with asterisks in Table  2, which have uncertain morphological classification. Green and pink symbols represent the mean values of spirals and ellipticals, respectively, as in Figs. 2 and 3. and because the previous work focused less on the kinematic signatures.
In summary, the results and limitations of this pilot study argue strongly for a much larger survey of S0s that spans a wide range of both masses and environments, and investigates the signatures of both kinematics and stellar populations in seeking to establish formation mechanisms. This approach will be the key to disentangling the contributions of mass and environment to dictating the channels by which lenticular galaxies form.  Table 2, which have uncertain morphological classification. Figure 9. Significance of the measured difference between λ(< R e ) for field and cluster galaxies, in the extreme assumption that all our sample is totally contaminated by spirals and ellipticals. The blue line at 1.43 represent the significance of our pilot project, in the sense that our measurements of λ(< R e ) for the field and cluster sample agree with each other if we multiply the error bars by 1.43. 7% of the simulations are above this line.   Figure A3. Ionised-gas velocity maps for the galaxies in our MUSE sample for which there is a significant detection of ionised gas.