Recovering the origins of the lenticular galaxy NGC 3115 using multiband imaging

A detailed study of the morphology of lenticular galaxies is an important way to understand how this type of galaxy is formed and evolves over time. Decomposing a galaxy into its components (disc, bulge, bar, ...) allows recovering the colour gradients present in each system, its star formation history, and its assembly history. We use GALFITM to perform a multiwavelength structural decomposition of the closest lenticular galaxy, NGC 3115, resulting in the description of its stellar light into several main components: a bulge, a thin disc, a thick disc, and also evidence of a bar. We report the ﬁnding of central bluer stellar populations in the bulge, as compared to the colour of the galaxy outskirts, indicating either the presence of an active galactic nucleus (AGN) and/or recent star formation activity. From the spectral energy distribution results, we show that the galaxy has a low luminosity AGN component, but even excluding the effect of the nuclear activity, the bulge is still bluer than the outer-regions of the galaxy, revealing a recent episode of star formation. Based on all of the derived properties, we propose a scenario for the formation of NGC 3115 consisting of an initial gas-rich merger, followed by accretions and feedback that quench the galaxy, until a recent encounter with the companion KK084 that reignited the star formation in the bulge, provoked a core displacement in NGC 3115 and generated spiral-like features. This result is consistent with the two-phase formation scenario, proposed in previous studies of this galaxy.

The formation history of NGC 3115 2147 for the larger number of high-mass galaxies in the local universe (Bernardi et al. 2010). However, the processes that lead to the formation of lenticular galaxies are still an open field of research. Different formation scenarios would leave different imprints on the galaxy kinematics and stellar populations. In lenticular galaxies that are the result of spirals undergoing ram-pressure stripping (Gunn & Gott 1972), we might expect the disc stellar population to be older than the one of the spheroid (Johnston, Aragón-Salamanca & Merrifield 2014), either because the gas truncation leads to an outside-in quenching of the star formation (Schaefer et al. 2017;Finn et al. 2018), or because the bulge stellar population could be rejuvenated by a last star formation episode (Bekki & Couch 2011), caused by the collapse of the stripped material into the central region. Lenticular galaxies that are formed through the merger and/or accretion of companion galaxies and satellites, might have younger discs than the bulge due to the accretion of gas-rich satellites on to the disc, which will enhance the star formation (Diaz et al. 2018). Finally, in lenticulars that are relics of high-redshift galaxies, such as primordial galaxies, which never passed through a spiral phase, the bulge, and the disc would have the same stellar population, as they are formed by the same material (Saha & Cortesi 2018).
The mechanisms that form lenticular galaxies in low-density environments would not be the same as the ones creating their highdensity counterparts, since mergers are more likely to occur in lowdensity environments (Tapia et al. 2017;Eliche-Moral et al. 2018). Tidal interactions are also typical of galaxy groups, while hot gas and a deep gravitational well are necessary to harass or strip a galaxy infalling a cluster of galaxies. Cluster lenticular galaxies might have been pre-processed in group or field environments, and it is therefore harder to isolate the mechanism responsible for their formation. On the contrary, isolated lenticular galaxies provide an ideal laboratory for studying the formation of this intriguing class of galaxies.
Using integral field unit (IFU) spectroscopy of a sample of 279 lenticular galaxies in the Mapping nearby galaxies at apache point observatory (MaNGA) survey, Fraser-McKelvie et al. (2018), have shown that lenticular galaxies have different stellar populations depending on their masses, such that low mass galaxies (M < 10 10 M ) are probably quenched spiral galaxies, whereas high mass lenticulars (M > 10 10 M ) are the result of merger events, either minor (Bekki 1998;Bournaud, Jog & Combes 2005) or major (Borlaff et al. 2014;Querejeta et al. 2015;Tapia et al. 2017;Eliche-Moral et al. 2018), independently of the environment. However, Coccato et al. (2020) showed that the kinematics of lenticular galaxies living in low and high density environments are different at 1.5σ confidence level at 1R e , suggesting that both mass and environment play a role in shaping the fate of the progenitor(s) of lenticular galaxies. It is still unclear, however, the exact role of these two factors in the formation of S0 galaxies.
In this work, we aim to answer the question of how the closest field lenticular galaxy, NGC 3115, was formed. This galaxy has been studied in a wealth of projects using different data sets (Guérou et al. 2016;Poci et al. 2019;Usher et al. 2019;Dolfi et al. 2020). NGC 3115 is the nearest galaxy hosting a billion solar mass black hole (Kormendy et al. 1996;Emsellem, Dejonghe & Bacon 1999) and several works have shown that it has a weak AGN activity (Lauer et al. 1995;Wrobel & Nyland 2012). Moreover, there is a general agreement that NGC 3115 formed through a two phase formation scenario, and it is quenched and isolated (Arnold et al. 2011;Guérou et al. 2016;Poci et al. 2019). Yet, several aspects of its formation are still unsolved, such as the presence of a low-luminosity AGN (Almeida et al. 2018), as well as of a large sample of UCDs and companion galaxies (Dolfi et al. 2020). Did they play a role in the quenching mechanism and how did they affect the galaxy kinematics and structure? For example, by recovering the galaxy kinematics using Planetary Nebulae, Cortesi et al. (2013) showed that the disc of NGC 3115 is rotationally supported, but that the value of V/σ is lower than what found for spiral galaxies, suggesting that this galaxy is not simply a faded spiral. Similar results are found by studying the globular cluster (GC) population (Zanatta et al. 2018) and integrated stellar light (Arnold et al. 2011).
Spectroscopic data allows for precise measurements of age and metallicity of the stars in a galaxy, but can generally only cover a small wavelength range. Spectroscopy is also limited in aperture, even IFU observations are usually restricted to the central parts of the galaxies, within a few effective radii. Furthermore, the large observing times required for IFU spectroscopy restrict this technique to small samples, whereas multiband observations may be more easily used to provide volume-limited, representative samples. A compromise can be achieved with narrow-band photometric surveys (e.g. Cenarro et al. 2019;Dupke et al. 2019;Mendes de Oliveira et al. 2019), which provide large number statistics and allow us to break some of the degeneracies involved in recovering the formation history of galaxies.
Spectral energy distributions (SEDs) obtained from broad-band photometry are far less restricted than IFU observations in spatial coverage, and the ability of this technique to gather information from different wavelengths, from ultraviolet to infrared, make up for the loss of detailed λ-by-λ constraints provided by spectroscopic data (Salim et al. 2007;da Cunha et al. 2012;Leja et al. 2017;Werle et al. 2019).
We provide in this work a study of a mimicked multiwavelength image, followed by a structural decomposition of NGC 3115. Previous works Vulcani et al. 2014;Dimauro et al. 2018Dimauro et al. , 2019Mosenkov et al. 2019;Psychogyios et al. 2020) have demonstrated the ability to decompose galaxies into their different structural components consistently and accurately, indicating that such decomposition may provide a more detailed picture of the formation history of the galaxy than by considering a galaxy as a whole. This pilot study of NGC 3115 was achieved by compiling ultraviolet, optical, near-infrared, and infrared data covering the galaxy out to 6.5 arcmin from its centre, thus studying the system out to ∼7 effective radii (Re = 58.1 arcsec; Cortesi et al. 2013), which provided a large-scale perspective of its main structural components. We then performed a detailed, consistent, multiwavelength decomposition of the main structural components of the galaxy. With that, we expect to enhance the understanding of the galaxy formation and dynamical history and the role of the environment on its evolution. This paper is structured as follows: In Section 2, we describe the data gathered for this study, as well as the process of the registration, calibration, and PSF homogenization. In Section 3, we explain the GALFITM setup and the details of the models used for the decomposition. In Section 4, we show all results obtained with GALFITM. In Section 5, we report our analysis. In Section 6, we discuss the results and propose a formation history scenario for NGC 3115.
Throughout this work, we adopt the distance to NGC 3115 of D = 9.4 Mpc (Brodie et al. 2014). We assume a standard CDM model, with H 0 = 70.5 km s −1 Mpc −1 (Komatsu et al. 2009).

Data sources and registration
To study NGC 3115, we use 11 images of the galaxy obtained with different telescopes, with the intent of covering a wide wavelength To create a consistent multiband image with images obtained with different instruments, we did a thorough treatment of the data before starting the fit. The images were all set to have the same pixel scale, world coordinate system (WCS), dimension and centre as the 2MASS H-band image (1 arcsec px -1 ). This was done using MONTAGE (Jacob et al. 2010), a tool for gathering fits images into custom mosaics using a reference image, in our case, the H band of 2MASS. We choose to use the H band as a reference, because the pixel scale is 1 arcsec, which is easier to deal with, and also this pixel scale is a compromise point between the small pixel scale of Subaru (0.2 arcsec) and large of WISE (2.75 arcsec). Differences in the gain, point spread function (PSF) and seeing are taken into account in the process of fitting with GALFITM. We created a synthetic PSF of each image using the Moffat distribution (Moffat 1969), which is based on a softened exponential profile to model the PSF. The parameters used to create such PSFs were obtained using the IMEXAMINE function in IRAF. To measure the sky background, we used the median value of small regions in the images, with little contamination from stars, followed by a sigma clipping, until convergence to the best value was reached (Hernandez-Jimenez et al. 2013. The best sky value for each wavelength band is defined and fixed in the configuration file of GALFITM to be considered by the routine when creating the best-fitting model. The zero points of these images were retrieved from each respective telescope, all in AB magnitudes. However, for the Subaru Suprime Cam images, these values were not available. To retrieve such information, we used the detailed study of the GC system around NGC 3115, which made use of the same images analysed in this work, presented in the paper of Jennings et al. (2014). Comparing the magnitudes of each object in this paper and our determinations of the instrumental magnitudes, we have estimated the offset, which corresponds to the zero points of these images. These are 31.40 ± 0.28, 30.5 ± 0.28, and 31.5 ± 0.76 mag for g, r, and i, respectively. We also recovered the zero points by using DECam images in the g, r, and i bands as a reference and we obtained the same values as from the GCs.

M E T H O D : S E T T I N G O F G A L F I T M
Historically, images of galaxies have been modelled in each passband individually, or aperture photometry was used to extract galaxy's luminosity into an area defined by the same radius in each band. This type of analysis might lead to inconsistent colours and photometric properties due to different image resolution and depth (see Vika et al. 2014). To overcome such a problem, MEGAMORPH Vika et al. 2013) was developed to perform simultaneous multiwavelength fitting of images. In particular, this technique allows to simultaneously fit bands that have lower signal-to-noise ratios (SNR), like ultraviolet and infrared data, with images with high SNR, as optical deep images, resulting in consistent and better-constrained parameters. This is obtained by expanding each fitted parameter in a wavelength-dependent series, using a smooth polynomial function (Chebyshev polynomial), where the maximum order of the series (hereafter orders) must be at most equal to the number of passbands. This method makes the fit possible and robust, even for low SNR bands. MEGAMORPH is a project that has produced a new version of GALFIT (Peng et al. 2002) and GALAPAGOS (Barden et al. 2012) routines, respectively GALFITM and  In this work, we make use of GALFITM to perform the analysis.

Tuning the initial configurations for GALFITM
The fundamental steps in recovering the galaxy components using GALFITM are: identifying the number of physically motivated components, identifying the order of the series associated with each parameter for every component, and defining the initial conditions for every component. On one hand, the highest the number of free parameters (and the number of components) and the highest the number of orders, the lowest will be the reduced χ 2 of the fit. On the other hand, the recovered components might be purely mathematical artefacts. Moreover, when including more than one component in the fit, we can assume that the colour variation inside each component is negligible when compared with the colour variation between different components. This approximation can be used to extract robust parameters from faint images, by anchoring the fit to the neighbouring deeper images. Again, to find a balance between freedom (i.e. large number of degrees of freedom) and no freedom at all (fixed aperture) is a fine-tuning process that requires trial and error. In the next subsections, we describe the two methods that we developed to find the best number of fitted components, best orders combination, and initial conditions of each component.

A frequentist approach: selecting the best-fitting model among 32 000 realizations
One important aspect in the GALFITM modelling is the order of the series allowed in the fitting process, which dictates how much each parameter may vary as a function of the wavelength. The higher the order of the series, the more flexible the function becomes, such that when the order is the same as the number of input images, that parameter can vary independently between bands as if the fit would be performed on every image separately. On the other hand, if the order of the series is zero, the best model would be exactly the one defined by the initial conditions in every band. Thus, to maximize the fit performance, we need to find the best combination of orders of each parameter.
To do so, we created a set of ∼32 000 GALFITM models using different combinations of orders that define the variation of each parameter with wavelength and an increasing number of fitted components [from a Single Sérsic (SS) model to a bulge-disc (BD) decomposition, to a bulge-disc-disc (BDD) decomposition]. The analysis of these 32 000 GALFITM models and the results obtained are shown in the Appendix A.
Throughout the process of finding the best model for NGC 3115, we inspected not only the best combination of orders for the parameters but also the best set of components that would better capture the galaxy's morphology. We investigated from the most extreme case of three Sérsic profiles until the other extreme of three exponential discs. The best-fitting model was the result of fitting a Sérsic profile and two exponential disc profiles. (A careful analysis of the results of the simulations of different combinations of orders is also shown in Appendix A.) By analysing the residual images and the BIC and AIC parameters, we obtained the following results (i) Three components are needed to recover the morphology of NGC 3115 using this data set.
(ii) The best fit is obtained using low orders for each parameter.
(iii) The best initial conditions are recovered.
Yet, the large number of realizations of the fit makes unfeasible a human evaluation of the residuals, which remains the most reliable source of judgement of the goodness of this type of fits. Therefore, we refine the fitting procedure as explained in the following section. We keep in mind, though, that the best-fitting model must comprise the three results obtained above.

An intuitive approach: using physical assumptions and visual inspection to find the best-fitting model
An alternative approach, rather than generating a high number of models to statistically encounter the best one, is to make physically based assumptions to define the fitting parameters. For example, when including more than one component in the fit, we can assume that the colour variation inside each component is negligible when compared with the colour variation between different components. Such assumption translates into using a low number of orders for the fitted parameters, as also occurring in the case of the best-fitting model obtained from the simulations. Moreover, the deepest and with highest SNR images, the g, r, and i images (Subaru), can be used to fix the shape of the galaxy components' at all wavelengths. The fit of the Subaru images started with a SS model, evolving to a BD decomposition, and finally to the BDD model -as found as the best model in the previous section. For the SS model, we retrieved the initial magnitudes from the Nasa Extragalactic Data base (NED) and the other initial parameters were set by the best SS model out of the 32 000. Then, we used the outcome of the SS model to set the initial conditions of the bulge-to-disc light decomposition and the results of the latter to define the initial conditions for the BDD model (see Fig. 1). Fig. 2 shows the evolution of the residuals for the three fits (SS, BD, BDD). We then adopt the best-fitting model obtained for the three deep images of Subaru to fix the shape of the galaxy's components at all wavelengths, returning to use the images homogenized with MONTAGE. In other words, we use the outcome of the fit of the Subaru images to define the initial conditions for the fit of the 11 bands, allowing the parameters to vary with 1 order (i.e. they must converge to the same value in all bands), except for the magnitude and disc-scale length, which we allow to vary with orders 9 and 2 to define the series, respectively, for all components. In fact, an order of 9 for the magnitudes is large enough to ensure an independent evaluation of the parameter, but still low enough to decrease the error bars and recover a smooth profile. While the choice of 2 for the order of the disc-scale length was made to improve the residuals.
The parameters found for this best model, as well as the allowed order of each parameter are summarized in Table 2, in the results section (Section 4). A summary of the steps taken to create the model that best describes NGC 3115 in all 11 bands is shown in Fig. 1.

A bulge, thick disc, and thin disc representation
As explained in the previous section, to find the best model, we passed from a SS model, to a BD model, to finally a BDD model, which provided the smallest residuals, while maintaining the physical reliability of the model, i.e. with physically meaningful components.
The shape and size of the components of a model can bring great insight to the actual morphology of a galaxy. For example, a classical bulge usually is modelled by a de Vancouleurs profile [i.e. Sérsic indec (n) 4], while a pseudo-bulge would have values of n 2 (Kormendy & Kennicutt 2004;Gadotti 2009). A thick disc would appear more oblate and extended than a thin disc, resulting in a higher value of b/a. Based on these considerations, we interpret our best-fitting model of NGC 3115, composed of three components, as a representation of a galaxy containing a classical bulge, a thick disc, and a thin disc.
Looking at the model images in Fig. 3 and given the high value of the thick disc axial ratio (b/a = 0.45; see Table 2), this thick disc could be interpreted as an oblate spheroid as proposed by Arnold et al. (2011) and Poci et al. (2019). From Table 2, we see that we obtain a Sérsic index of n = 3.5 for the bulge, which is also in agreement with the black hole mass of NGC 3115, retrieved by Kormendy & Kennicutt (2004; 2 × 10 9 M ), according to Savorgnan et al. (2013) and Graham & Driver (2007). Finally, in Table 3, we compare the total magnitudes obtained in our GALFITM model with available measurements in NED. We apply a conservative error of 10 per cent to our magnitudes, since this result is the combination of the magnitudes obtained for each component. We can see that, within errors, our GALFITM model reproduces well the expected magnitudes. The only exception is the total magnitude obtained for the FUV band using a BDD model, which is almost one magnitude brighter than the magnitudes available in the literature. The FUV image is too faint to fit three components and the model is oversubtracting the data, creating negative residuals, as possible to see in Fig. 3. In the rest of the work, we will use the SS value for the FUV magnitude, when possible. In the other cases, we will test the difference between including or not the FUV value in the fit.
In Fig. 3, we show the SS, BD, and BDD decomposition of NGC 3115. In it, we compare the accuracy of the fit when using one, two, and three components. First, we show in the right-hand side of the plot the quadratic sum of the standard deviation of the residuals over all bands, which consistently shows how the model with three components significantly improves the fit of the galaxy, achieving the smallest residual. We also show, for each band, the standard deviation of the residuals after the BDD decomposition, revealing that the bands that have the largest residuals are the ones in the optical, as it would be expected, since these are the deepest images, and as previously discussed, only in these bands we are able to identify and even fit extra components, such as a central bar (see Fig. 9) and spiral arms (see Fig. 10). The other bands are so faint that these extra components are not identified and so their residuals are smoother, since the main components of the galaxy are well fit.   Table 2, as well as the order of the series of each parameter.
In Fig. 4, we show the surface brightness profile of NGC 3115, highlighting where each component dominates the galaxy in two wavelength regimes: UV (NUV) and optical (r), as well as the profile recovered calculating the light encompassed in ellipses of fixed axial ratio obtained using MORFOMETRYKA (Ferrari, de Carvalho & Trevisan 2015). One important comment on Fig. 4 is that even if the modelled bulge is not dominant at large radii, it is a component that contributes from the innermost region of NGC 3115 until its outskirts, since by construction, the GALFITM fit extends to infinite, becoming negligible where the galaxy light reaches the background value. As a consequence, the galaxy components also extend to infinite. The relative contribution of each component to the entire flux varies with radius and from galaxy to galaxy. Allen et al. (2006), modelling the light distribution in 10095 galaxies from the Millennium Galaxy Catalogue (MGC), identified six types of surface brightness profiles. In the type 3 case, the Sérsic profile accounts for the bulge light and the outer part light, being the disc embedded in the spheroidal component. The question whether the outer spheroid is simply an extension of the inner bulge or not, although interesting, is out of the scope of this paper. This method is different from other structural decomposition methods (e.g. Poci et al. 2019) that cut the domain of the bulge to a defined radius, and, therefore, are not directly comparable with this work.
The comparison with the light profile recovered with MORFOME-TRYKA shows that the total (bulge + thin disc + thick disc) GALFITM model is recovering relatively well the galaxy light both in NUV and r bands, shown as an example. The displacement between GALFITM and MORFOMETRYKA at small radii can be explained by a Sérsic core displacement (Graham & Driver 2007;Dullo & Graham 2014), this is only apparent with MORFOMETRYKA, since this method fits isophotes to recover the galaxy light, while GALFITM fits an extrapolation of the bulge Sérsic profile. However, an important difference between the two methods is that in the models recovered by GALFITM, the axial ratio varies for different components, and the axial ratio of the total model is the sum of the contribution of the bulge, thin disc, and thick disc, where the latter might vary with radius as well, depending on which component dominates the light profile. On the other hand, MORFOMETRYKA calculates the light embedded in ellipses of fixed ellipticity. This difference does not seem to cause discrepancies, as the light profiles recovered with both methods are consistent within errors at all radii in the two wavelength regimes, probably as a consequence of the fact that the thick disc is the predominant component, in this radial range, and its axial ratio is the same obtained with MORFOMETRYKA. With this test, we show that the total parametric model of GALFITM does not miss any light when compared with the flux enclosed in ellipses following the isophota profile.

The residuals
The residuals of NGC 3115, i.e. the image obtained after subtracting the best-fitting model from the input images, are very peculiar and show the same patterns independently of the code used to fit the image (see Arnold et al. 2011;Guérou et al. 2016;Poci et al. 2019), raising the question of whether there are other components in the galaxy.
To clearly understand the residuals, we provide in Fig. 2 the evolution of the residual images from the simplest GALFITM model with one component to the one with three components for the r band. By looking at Fig. 2, we begin to understand the necessity of a model with at least three components to describe NGC 3115. The left-hand panels show the residual images resulting from fitting a single component: only the outer parts of the galaxy are fitted, leaving behind two semicircular symmetric regions in the major axis and a central equatorial component. In the two-components model (central panels of Fig. 2), the outer residuals around the major axis disappear, but the central equatorial component is still not fitted. Note that this equatorial component extends until nearly 1R e . Finally, adding a third component (right-hand panels of Fig. 2), which we interpret as a thin disc, the model is capable of fitting a significant part of this equatorial component seen in both previous models. The fit of this third component reveals the existence of other two components on NGC 3115: two structures distributed along the major axis of the galaxy that can be either remnant of spiral arms or an edge-on ring; and a central hourglass shape, oriented along the minor axis of the galaxy, resembling an end-on bar. These findings are consistent with Capaccioli, Held & Nieto (1987) and Guérou et al. (2016), which proposed that the outskirts of the galaxy exhibit remnant spiral arms, which can be seen as very shallow symmetric structures in almost every residual image of the galaxy (see Fig. 2). Guérou et al. (2016) also proposes that the galaxy holds a bar and a nuclear ring in its centre.
As mentioned, analysing the residual image, we also observe a central component similar to a bar. Adding a fourth component to our GALFITM fit (a Sérsic model with Sérsic index lower than 1, typical of a bar), we were able to fit such structure. However, we could only add this fourth component in the model of the Subaru images, because of the higher SNR (see Section 6.1 for more details).
It is hard to decouple the presence of the thin disc from a more general 'equatorial' component, which would include not only the thin disc but also a nuclear ring, an inner bar, spiral-features, i.e. all the components lying along the central part of the major axis of the galaxy (Arnold et al. 2011;Guérou et al. 2016). In an attempt to quantify the contribution of the equatorial component, we obtained the flux of this inner region using the package ISOPHOTE from the PHOTUTILS (Bradley et al. 2016) library in PYTHON, a tool to fit elliptical isophotes to a galaxy image, i.e. we fitted the total model -bulge -thick disc to quantify the flux of this inner region, and we compare it with the flux of the thin disc, as recovered by GALFITM. The difference between the two is lower than 0.23, 0.16, and 0.27 mag in the g, r, and i bands, respectively (which is consistent with the errors in the GALFITM model), so we conclude that the thin disc as extracted with GALFITM is a fair representation of the central region of the galaxy. In Section 6, we discuss in detail the result of the ellipse fitting for the different sub-components.

Colours
Our first step to understanding the stellar populations present across the galaxy is deriving 1D colour gradients of our models. In Fig. 5, we show different colours of the galaxy retrieved from the GALFITM models and input images, i.e. FUV-NUV, .
To create this figure, we draw a rectangle of 50 arcsecs height along the galaxy's major axis, from the innermost pixel of NGC 3115 until its outskirts, covering out to 3 arcmin, and we obtain the median value of the colour at each radius (in bins of 1 arcsec). In this figure, the red vertical line represents the end of the bulge domain in the r band (21.69 arcsec) as derived from the bulge model of GALFITM (see Table 2), and the yellow vertical line shows the bulge-to-disc transition (i.e. where the bulge light becomes subdominant in the fit) derived from Fig. 4.
To create this figure, we use both the GALFITM models (blue line) and the input images (green line), for comparison.
In Fig. 5, we show the colour gradients obtained for the two cases and we see that we recover, within errors, the same trends using the input image and modelled image, indicating that we are retrieving reliable colour gradients present in the galaxy. From Fig. 5, we can observe, moreover, the power of GALFITM in delivering a smoother colour profile, and enhancing regions of transition and inversion of the colour gradient.
It is worth noticing in Fig. 5 that the region that demarcates the transition from the bulge to the disc domain (extracted from Fig. 4) coincides with the inflexion in the colour gradient in the optical bands, suggesting the presence of different stellar populations between the components. Moreover, it is interesting to see that the galaxy becomes bluer with radius in the optical, indicating that we may be observing the effect of several low-mass accretions in NGC 3115 that would account for the bluer colours in the outskirts of the galaxy (Mapelli M. et al. 2015).
Moreover, from the GALFITM model gradients in Fig. 5, it is clear that in the ultraviolet range, the inner region is bluer than the outer one. Similar behaviour is present in the optical colour, but it is inverted at large radii, being the transition point the radius where the bulge light gets subdominant. The near-IR colour gradient is opposite to the FUV-NUV and NUV-r ones, being the inner region the reddest.
According to Kaviraj et al. (2007), a recovered NUV-r 3.3 colour would imply that the bulge of this galaxy hosts a young stellar population, as discussed in details in Section 6.3. Finally, the optical colours g−r and r−i are consistent with the results obtained by Guérou et al. (2016), Poci et al. (2019), and Dolfi et al. (2020): from the outskirts inwards, we detect a bluer outer envelope, followed by a redder disc-like region and by a mid-colour inner region.
The bluer colour of the bulge in the ultraviolet regime could be attributed to a young stellar population, the presence of an AGN (Wong et al. 2011) or the emission of white dwarf stars (Lisker, Grebel & Binggeli 2008). However, the latter is less likely, because white dwarfs would yield bluer FUV-NUV, instead of NUV-r, and the contribution of white-dwarfs emission would not be enough to explain alone the retrieved colour gradients in the FUV-NUV and NUV-r images. Moreover, this galaxy presents extended strong emission in the Spitzer 8 and 24 micron bands, covering especially the central part of the galaxy. The presence of the mid-infrared emission in passive early-type galaxies (ETGs) might be of stellar . Annotated in red in the seventh row are the standard deviation of the residuals for each band in the BDD model. Eighth, ninth, and tenth rows correspond to the bulge, thick disc, and thin disc models present in the BDD model, respectively. The columns represent each wavelength band. Annotated in red in the right-hand side of the image are the quadratic sum of the standard deviation of the residuals over all bands for the SS, BD, and BDD models. 10.0 ± 1.0 r 9.0 ± 0.9 9.05 ± 0.03 i 9.0 ± 0.9 8.72 ± 0.03 z 8.7 ± 0.9 -J 6.7 ± 0.7 6.80 ± 0.02 H 7.4 ± 0.7 7.9 ± 0.02 K 7.6 ± 0.8 7.6 ± 0.02 3.4 7.6 ± 0.8 -4.6 8.3 ± 0.8 -Note. a this magnitude was retrieved with the SS model. origin, coming from the dusty envelopes around the asymptotic giant branch (AGB) stars (Clemens et al. 2011) or could, again, indicate the presence of an AGN. In fact, there is evidence for the presence of a billion mass supermassive black hole (2 × 10 9 M ) in this galaxy and hot flows in the centre, studied in detail by Lauer et al. (1995), Kormendy et al. (1996), Emsellem et al. (1999), Almeida et al. (2018), and Jones et al. (2019). By looking again at Fig. 5, we see that using optical and near-IR colours, we retrieve the same gradient as the one found by Guérou et al. (2016), i.e. a bulge with an older/more metallic stellar population than the disc. However, when we include ultraviolet colours, we can identify these blue colours in the central region of the galaxy. Several pieces of evidence strengthen this hypothesis: (1) The very central region of the galaxy (identified as blue in the NUV−r panel) is the only one that shows X-ray emission.
(2) This same region shows high dust temperature when observed with Spitzer (Amblard et al. 2014), indicating that the UV light was likely absorbed and re-emitted by the dust, confirming the presence of strong UV light in this region (strong evidence of young populations), and (3) when focusing only on optical and near-IR colours, our results agree well with previous works.
To understand the origin of the blue colour in the central region of the galaxy, i.e. if it is due to the presence of an AGN or recent star formation activity, we perform a fit of the SED of the galaxy and its components, as described in the next section.

SED fitting
In this section, we describe the SED fitting process implied to derive the physical properties of the galaxy and its components. For such task, we used the code CIGALE (Burgarella, Buat & Iglesias-Páramo 2005;Noll et al. 2009).

CIGALE
CIGALE is a robust SED fitting code, which includes the contribution of dust emission and attenuation models, AGN models, as well as the simple stellar population models from Bruzual & Charlot (2003), providing accurate estimates of the galaxy properties.
Since we aim to understand if the blue colour observed in the bulge is the result of recent star formation events in this region or emission of the AGN (or a composite of the two), we use CIGALE to perform . Surface density profiles of NGC 3115 and its components for two wavelength regimes: ultraviolet (NUV) and optical (r). Dotted red, cyan, and yellow lines represent the extrapolation of the bulge (Sérsic), thick disc (exponential disc), and thin disc (exponential disc) fit performed with GALFITM to the stellar surface brightness data. The surface brightness profile as recovered with MORFOMETRYKA (MFMTK) fitting (black dotted line) and with GALFITM models with two (BD) and three (BDD) components (blue and magenta lines, respectively) overlap until the faint outskirts in NUV and from 60 arcsec in the optical. the fitting accounting for models that include both the contribution of stellar and AGN emission.
We fit a total of 11 free parameters to the SED of the galaxy and its components: four parameters for the star formation history (SFH), one parameter for the dust attenuation, one parameter for the simple stellar population, and five parameters for the AGN emission templates, described in Table 4.
In the SED fitting process, we include the AGN emission models of Fritz, Franceschini & Hatziminaoglou (2006) only in the bulge component and we use the stellar population models of Bruzual & Charlot (2003). We assume the dust attenuation model described by Calzetti et al. (2000), a Chabrier initial mass function (Chabrier 2003), and a double exponential SFH, which, according to Burgarella et al. (2005), provides a good estimate of the SFR of ETGs and galaxies that had a recent episode of star formation, but is not good for recently quenched or starburst galaxies. The concept behind this star formation history is to allow for two decaying exponentials to fit the star formation rate of the system, one for the long term SFR and one for the more recent events. We point out, however, that it is not because we are using a double exponential SFH that both peaks will have significant participation in the fit. In fact, if the data requires, one of the peaks of the SFH has the freedom to be zero.
It is especially important to include AGN emission in the SED fitting of the bulge component, as it can have a large impact on the derived properties of the component and bring light to the role of the AGN in the evolution of lenticular galaxies. For example, in the UV and optical wavelength range, a large portion of the emitted light of a galaxy hosting an AGN can be attributed to the accretion disc, and could thus explain the blue colours that are observed in the bulge of NGC 3115.
In Table 4, the fitted parameters and the values that each parameter can assume in a CIGALE SED modelling are shown. The total number of models, shown in the bottom line of Table 4, states the total number of SED models that are created by varying and combining these parameters. The parameters described in Table 4 are: τ main , the e-folding time of the main stellar population model, f burst is the mass fraction of the late burst population, age burst is the age of the late burst, and τ burst is the e-folding time of the late starburst population model. While for the dust attenuation, the only fitted parameter is A v , the V-band attenuation of the young population. For the AGN, τ stand for the optical depth at 9.7 microns, while β is linked to the radial dust distribution in the torus and γ is linked to the angular dust distribution in the torus. The opening angle regards the full angle of the dust torus (see Fig. 1 of Fritz et al. 2006), and the AGN fraction relates to the percentage of the AGN mass concerning the total mass of the galaxy.
To find the best-fitting SED, we performed millions of tests using CIGALE, using both the delayed SFH and double exponential SFH. What we see is that both models fit well the data, however, the delayed SFH fails to reproduce simultaneously the UV (young stellar population) with the NIR and IR (old stellar populations), delivering sometimes too high SFR for very old stellar populations. As discussed in Ciesla et al. (2016) and Ciesla, Elbaz & Fensch (2017), it is very hard to recover both age and SFR reliably using SED fitting techniques. Moreover, there is a general agreement on the difficulty of constraining the age of galaxies from broad-band SED fitting (e.g. Maraston et al. 2010;Pforr, Maraston & Tonini 2012;Buat et al. 2014;Ciesla et al. 2015). What we find, however, is that using a double exponential SFH, the delivered ages are closer to the ones expected from other works while maintaining the physical reliability of the SFR.
Therefore, we decided to use the double exponential SFH for all SEDs shown in this work.
The resulting SEDs of the galaxy and its components are shown in Fig. 6, while the SED of the bulge including AGN models is shown in Fig. 7. In these figures, we also indicate the contribution of dust attenuation and stellar emission.
Furthermore, to better constrain the SED results and be certain, we are obtaining fiducial results, we took two complementary roads: (1) we added to our fitting procedure two literature data in the infrared regime (Spitzer 8μ, IRAS 12μ; Amblard et al. 2014;Knapp et al. 1989).
(2) we compare the results obtained with GALFITM with those obtained with ELLIPSE in IRAF. The results of this comparison, with the infrared complementary data included in the fitting, are shown in Fig. 8.
The results presented in Figs 6, 7, and 8 are summarized in Table 5, showing the main physical properties obtained for the galaxy and its  components, including and not including AGN models in the fitting process. In the two bottom rows, we show the results when including literature data to constrain the SED fitting and compare it with the fit obtained with ELLIPSE in IRAF. From Figs 6 and 7 and Table 5, we begin to understand the influence of the AGN in the emission of the galaxy bulge at different wavelength ranges. As other works have shown before (Guérou et al. 2016;Dolfi et al. 2020), the galaxy is overall old and with low amounts of gas in all components. However, using photometric data in the UV, we detect blue colours in the bulge of the galaxy (5.1), maybe revealing the presence of younger/less metallic populations in its central region. The emergence of this blue colour could be explained by a recent star formation event in the centre of the galaxy, or it could be due to AGN emission. What we see, however, in Fig. 7, is that the AGN emission has very low luminosity (Almeida et al. 2018;Jones et al. 2019) and it is always subdominant in the galaxy, especially in the UV regime, where this emission seems not to play any role. This indicates that the bluer colours observed in the bulge can not be due to the presence of the AGN. Moreover, looking at Table 5, we see that the derived properties obtained for the bulge including and not including AGN models are similar, within errors, reinforcing that the presence of this AGN is not interfering in the physical properties of the bulge. Note that the reduced χ 2 is smaller in the SED fitting including the AGN model.
From Table 5, we see that the mass of the galaxy is mainly in the bulge, where the supermassive black hole lies, in contrast with what was found by Poci et al. (2019) that the majority of the galaxy's mass is in the spheroidal/thick disc component. This discrepancy might be related to the different methods used to recover the galaxy's components. In fact, in Poci et al. (2019), an orbit-superposition dynamical modelling is used, while GALFITM fits the galaxy light distribution. This implies that the Sérsic model used to recover the bulge component would be dominant in the innermost region, but would keep a contribution until very large radii, resulting in the high mass recovered for the bulge component. On the other hand, Poci et al. (2019) associate all the outer part of the galaxy to the spheroidal/thick disc component. To investigate the validity of the decomposition carried out in this work, we calculate the mass-to-light (M/L) ratio for the bulge and thick disc components using the 3.4 μm band. The results obtained are comparable to what found by the DiskMass Survey (DMS; Bershady et al. 2010;Lelli et al. 2016), being M/L bulge = 0.26 for the bulge and M/L disc = 0.16 for the disc. Moreover, we find that M/L bulge = 1.6 M/L disc , nearly as suggested by SPS models (Schombert & McGaugh 2014), granting the feasibility of the recovered parameters. Thus, we conclude that our decomposition is valid. Interestingly, Lelli et al. (2016), for Spitzer data at 3.6 μm, find that M/L bulge = 0.7, M/L disc = 0.5. Keeping the light fixed, these values of M/L would imply almost double the mass recovered with SED fitting in this work, so nearly the same mass obtained by Guérou et al. (2016) and Poci et al. (2019).
In fact, when comparing the total stellar mass of the galaxy with the results obtained using MUSE data (Guérou et al. 2016;Poci et al. 2019), we can see that our measurement is lower, even using a bigger field-of-view. In this respect, it has been shown in several works (e.g. Ciesla et al. 2015 and references therein) that SED fitting techniques deliver stellar masses systematically lower than their spectroscopic counterparts. Particularly, using CIGALE, Buat et al. (2014) and Ciesla et al. (2015) found that the output stellar mass, indepen-dently of the SFH chosen, is approximately 7 per cent lower than expected.
When comparing the results obtained with GALFITM for the total model of NGC 3115 (bulge+thick+thin disc) with results obtained with ELLIPSE in IRAF, we see that the results we are obtaining with GALFITM do reflect realistic features of the galaxy. One important reinforcement to make in this comparison is that ELLIPSE integrates the light of the galaxy until the outermost detectable (with respect to the background light) isophote, whilst GALFITM integrates until infinite, so the magnitudes obtained are expected to be slightly fainter for ELLIPSE, and the SED results might reflect these differences. Moreover, while the fitting with ELLIPSE integrates the flux inside several ellipses, GALFITM recovers the flux of the parametric models that better recover, when combining all the components, the light distribution of the galaxy. Thus, this comparison is only valid in order of magnitude, as a consistency check. Looking at Table 5, it is possible to see that the total physical parameters obtained with the two models are in agreement, inside error bars. The results of stellar mass obtained using both GALFITM and ELLIPSE are consistent with the ones obtained by Forbes et al. (2017) using Spitzer imaging.
Therefore, excluding the mass, the stellar population properties recovered in this work are in agreement with what found by Guérou  2016) using MUSE data. One should keep in mind, none the less, that the two methods are not directly comparable 1 ; for instance, it has been shown that SED fitting techniques tend to underestimate the recovered ages (Ciesla et al. 2018). Nevertheless, both works find that the stellar population of the galaxy is on average around 11 yr, with slight variation along the major axis, being the innermost part the oldest and more metallic, followed by a region with a mixed population (young and star-forming/old and quenched) and a middleaged, low metallicity outer region.
Finally, using the recovered stellar mass and SFR, we obtained the specific star formation rate (sSFR) for each component. We find that the thin disc has sSFR of around 2.7 × 10 −10 yr −1 and it is higher than the sSFR of the thick disc and the bulge, see Table 5. Following the definition of Bruce et al. (2012; sSFR <10 −10 yr −1 is quiescent), the thick disc and the bulge are quenched, while the thin disc is moderately active, again in agreement with what found by Guérou et al. (2016). Interestingly, even if quenched, the bulge presents a recent event of star formation, as we shall discuss in detail in the next section.

The inner bar
As mentioned in Section 4.2, the model in the optical bands is significantly improved with the addition of a bar, however, the inclusion of this bar to other bands with lower SNR caused a strong over subtraction and so we decided not to include this component in the final model. Looking for alternative ways to prove the hypothesis of the presence of this bar, we compared NGC 3115 with two simulated galaxies: both with the same fitted components as NGC 3115 (bulge, thick,and thin disc), but one with a bar and one without a bar to see which had more similarities (and showed similar residuals to those) of the model. The unbarred simulation uses the star-forming simulations described in Clarke et al. (2019) and Beraldo e Silva et al. (2020; who refer to it as FB10), which forms clumps at early times producing a thick disc in the system. The barred galaxy, on the other hand, is a pure N-body (no gas or star formation) simulation initially comprised of a thin disc+thick disc. This model is similar to the thin disc+thick disc simulations described in Debattista et al. (2017), but with the initial thick disc having a scale-height of 900 pc. After forming a bar, we view the system edge-on with the bar end-on.
For the comparison, we have created unsharp masks for the r band of NGC 3115, the unbarred galaxy and the barred galaxy. The unsharp mask is an image sharpening technique, which consists of digitally blurring or 'smoothing' the original image. This operation suppresses smooth features (i.e. which have a structure on large scales) in favour of sharp features (those with a structure on small scale), resulting in a net enhancement of the contrast of fine structure in the image. The unsharp masks used in this work were created using a 25-by-25 pixels median box to create the median subtracted image and a circular Gaussian smoothing with σ = 5.
The top panel of Fig. 9 shows the unsharp mask of the simulated unbarred galaxy, where it is possible to see a stripe of stars in the centre of the galaxy. The simulated barred galaxy is shown in the central panel of Fig. 9 and the unsharp mask of the r band of NGC 3115 in the bottom panel of Fig. 9. It is clear that the hourglass shape present in the unsharp-mask of the r-band image is also present in the simulation of the barred galaxy, suggesting that this feature in NGC 3115 is probably caused by the presence of a bar. Guérou et al. (2016) also reported the possible existence of an end-on bar, based on the correlation between the mean and skewness of the measured line-ofsight velocity distribution (V and h3), usual of barred-galaxies.

Recent dynamical history: encounter with a dwarf companion
NGC 3115 is known to have a very rich GC system (Cantiello et al. 2018) and several companion galaxies, which are mostly ultracompact dwarfs (Dolfi et al. 2020). KK 084 is one of the dwarf spheroidal companions of NGC 3115, located at about 5.5 arcmin to the SE, (Harris & van den Bergh 1981;Puzia & Sharina 2008). KK084 is very faint and, from all data sets available to us, it could only be detected in the deep Subaru images. Using these images, we fitted KK 084 with ELLIPSE in the PHOTUTILS library in PYTHON, to recover the main physical properties of this dwarf galaxy. Also with ELLIPSE, we fitted the observed spiral-like features in the residuals of the models obtained with GALFITM. In total, we fitted three components with ELLIPSE, the companion and the NE and SW spiral arms. The results of this fit are shown in Fig. 10 and in Table 6.
Analysing Fig. 10, it is intriguing how the colours of the two hypothetical remnants (or reborn) spiral arms differ from each other. This difference may be due to the very low SNR of the images of both structures, or, given that one arm is much redder than the other, it could be due to the presence of dust. However, only new deep multicolour imaging (μ r 28) could provide a robust understanding of the stellar populations present in these sub-components. Now, looking at Table 6, we notice that the g−r and r−i colours of the companion galaxy differ approximately 0.3 dex from the colours of the components of NGC 3115 (see Fig. 5), indicating that they probably do not have similar stellar populations. We could speculate that the companion galaxy suffered a recent star formation episode caused by the interaction with NGC 3115, given its blue colours. Unfortunately, this galaxy is too faint to be observed in the other Figure 8. SED of the NGC 3115 photometry obtained using GALFITM and ELLIPSE, respectively. The first eleven points in our fit represent the data used throughout this work, while the final two points are literature data from Spitzer 8μ and IRAS 12μ. In each panel, we show the retrieved physical properties. The magenta squares stand for the input fluxes, the cyan dots are the model fluxes, the black line is the modelled spectrum, the orange line is the stellar attenuation, the blue dashed line is the stellar emission unattenuated and green line is the AGN emission. Red triangles are the literature data. Table 5. SED fitting results obtained from CIGALE. The properties obtained with CIGALE are divided in three: containing AGN models in the fitting process, without AGN models, and including literature infrared data.
bands or to perform SED fitting to get the SFR and age, so to confirm or discard this possibility. Despite the need for deeper data to understand the stellar populations of both the companion KK084 and spiral arms and their relation to a possible encounter with NGC 3115, we can still probe the hypothesis of this encounter by analysing if there exists an offset in the core of NGC 3115, typical of recent interactions. In fact, it has been shown that the external parts of bulges and elliptical galaxies can have a crossing time longer than the time of an interaction (Combes et al. 1995;Binney & Tremaine 2008). Therefore, these regions suffer a tidal impulse in response to an encounter (Aguilar & White 1985;Binney & Tremaine 2008) rather than a differential tidal one, in which the internal dynamical time is shorter than the encounter time. Then, the external parts, due to an impulse answer, expand, while the core is displaced, as a tidal typical answer. This offset can be up to 20 per cent of the observable radius of the galaxy (Lauer 1986(Lauer , 1988Davoust & Prugniel 1988;Combes et al. 1995;González-Serrano & Carballo 2000;Mora et al. 2019). Such offsets are excellent indicators of recent galactic encounters (Combes et al. 1995;Mora et al. 2019). To search for signs of recent interactions, we performed a photometric analysis following Mora et al. (2019), who analysed the distortions of the Penguin system (Arp 142). We modelled the i-band image of NGC 3115 with the ELLIPSE task in IRAF, letting all geometrical parameters, such as position angle, ellipticity, and centre of the ellipses unconstrained. In order to quantify the displacement, we used δ/a ratio (Lauer 1988;Mora et al. 2019), where δ is the separation of the isophotal centre to the nucleus and a is the length of the semimajor axis (SMA). In Fig. 11, we show the radial profile of δ/a and the outermost isophotes of NGC 3115 and respective ellipses Figure 10. ELLIPSE fitting of the spiral-like features observed in NGC 3115 and its dwarf companion KK084, located approximately at (5',-2'). The rows represent the g, r, and i bands, respectively. The columns show the input image (residual from the fit with GALFITM), the model of the spiral features and companion KK084 and residuals of the fit, respectively. Table 6. Colour of the components of NGC 3115 and the dwarf companion KK084 recovered using GALFITM and ELLIPSE. Component Total Model (GALFITM) 0 . 9 9 ± 0.12 1.05 ± 0.14 0.05 ± 0.01 Bulge (GALFITM) 0 . 9 5± 0.13 1.03 ± 0.13 0.08 ± 0.02 Thick disc (GALFITM) 0 . 9 5± 0.16 0.99 ± 0.14 0.05 ± 0.01 Thin disc (GALFITM) 1 . 1 1± 0.14 1.17 ± 0.19 0.06 ± 0.02 KK084 (ELLIPSE) 0 . 6 3± 0.07 0.68 ± 0.09 0.05 ± 0.01 NE spiral feature (ELLIPSE) 0.93 ± 0.14 1.92 ± 0.34 0.77 ± 0.11 SW spiral feature (ELLIPSE) 0.2 ± 0.01 − 0.16 ± 0.03 − 0.36 ± 0.07 fitted, and their centres. We can see, in the bottom panel of Fig. 11, that from about 160 arcsec from the centre of the SMA, the centre of the isophote starts to show an offset towards the South, revealing a maximum displacement of about 5 per cent (≈800 pc). Actually, the direction of the offset could serve as a constraint to the location of the pericentral passage, e.g. Combes et al. (1995) and Mora et al. (2019). This is a strong indication that NGC 3115 had a pericentre passage with a companion galaxy just a few hundred thousand years ago (Combes et al. 1995;Mora et al. 2019). A natural candidate as the disturbing agent would be KK084. From the magnitudes measured by Sharina et al. (2008), we have estimated a mass for KK084 of M = 1.5 × 10 8 M . Such a low mass would imply a high-speed encounter with a small impact parameter for the interaction between NGC 3115 and KK084 (Aguilar & White 1986). Moreover, KK084 might have lost part of its mass during previous passages around NGC 3115 (Aguilar & White 1986;Mayer et al. 2001).

The formation history of NGC 3115
NGC 3115 is a unique and extremely interesting object. Studying this galaxy, we found several extra components to the system, even when other lenticular galaxies usually present only a bulge, a disc, and a bar (or a lens; Johnston et al. 2021). The existence of these extra components in NGC 3115 (nuclear ring, end-on bar, spirallike features) reinforce a complex formation scenario for this galaxy. Including data from an extended wavelength range in the analysis brings us important information regarding the stellar populations of each of these components. As previously discussed, NGC 3115 was largely studied in the past, using photometric, spectroscopic, discrete tracers and IFU data. Cortesi et al. (2013), using PNe as tracers of the overall stellar kinematics, were able to rule out the scenario in which NGC 3115 is simply a stripped spiral galaxy, given the high value of V/σ 3.3 in the thick disc of the galaxy. The same result was obtained by studying the GCs kinematics (Zanatta et al. 2018). Reinforcing this point, Falcón-Barroso et al. (2019) show, using CALIFA data of 300 galaxies across the Hubble sequence, that spiral galaxies hardly fade into S0s, unless specific conditions of light concentration, angular momentum, and V/σ are satisfied (as for Sa galaxies). In the other cases, pre-processing might have taken place in the formation of the lenticular, either through tidal interactions or mergers.
Various works proposed the two-phase formation scenario as the most plausible to explain the formation of this galaxy (Arnold et al. 2011(Arnold et al. , 2014Guérou et al. 2016;Dolfi et al. 2020). Arnold et al. (2011) described this formation scenario as the central region being formed in a violent and dissipative event, such as a gas-rich major merger, followed by a sequence of dry minor mergers and accretion events, that would have been responsible for the assembly of the outer spheroid and halo of the galaxy. Such a scenario is supported by the declining rotation and metallicity profiles measured for this galaxy, in their work, using optical imaging and multislit spectra. Guérou et al. (2016), using IFU data, recovered the same kinematics and metallicity profiles as Arnold et al. (2011;. In their paper, Guérou et al. (2016) also hinted at the possibility of new components to this galaxy. They suggest the presence of an end-on bar, a nuclear disc (suggested by the correlation between the mean and skewness of the measured line-of-sight velocity distribution), and of remnant spiral arms in the faint end of the thin disc of the galaxy. Poci et al. (2019) suggest that the bulge and subsequently the thin disc of NGC 3115 could have been formed via a 'Compaction' (Dekel, Sari & Ceverino 2009;Zolotov et al. 2015) scenario at early times, where a dissipative contraction would bring gas on to the central region of the galaxy via cold streams, which could trigger star formation events in its centre, being rapidly quenched in such massive galaxies as NGC 3115 (Zolotov et al. 2015). They also confirmed that the dynamics and stellar population of the thick disc is the result of the accretion of dwarf companions and satellites, which may increase the thick disc dispersion velocity and flatten the metallicity gradient. In a more recent study, Dolfi et al. (2020), using a large amount of data on NGC 3115 and faint companions, covering up to 6.5 effective radii and implementing an innovative kinematic method, confirmed the assembly history proposed by Arnold et al. (2011), finding an older pressure-supported bulge, followed by a transition in the kinematic profile at 0.2 Re to a rotationally supported disc, that goes out to 2-2.5 effective radii and then changes again to an outer pressure-supported spheroid. Such transitions agree well with the two-phase formation scenario.
Moreover, the stellar halo mass of NGC 3115 has been estimated to be nearly 14 per cent of the total stellar mass of the galaxy (e.g. Peacock et al. 2015). This could be strong evidence that a merger had a key role in the formation of this galaxy. The formation of the thick disc could be attributed to this initial merger, although other mechanisms could also be responsible for that, e.g. clump instabilities in the early disc -see Bournaud, Elmegreen & Martig (2009), Clarke et al. (2019), and Beraldo e Silva et al. (2020. The results of this work, from the analysis of the SED fitting of the galaxy's components and the galaxy colour gradients, agree with the previously proposed scenario for the formation of NGC 3115. We see that all of the components of the galaxy are old, as it is retrieved by Guérou et al. (2016) and Dolfi et al. (2020), and we also recover the declining metallicity profile from the bulge out to the thick disc (outer spheroid), as recovered by Arnold et al. (2011) and Poci et al. (2019).
One difference, however, between our work and the previously cited works, is that the inclusion of images in the ultraviolet has suggested the presence of blue stellar populations in the bulge of the galaxy, that we propose to be the result of star-formation Figure 12. Formation scenario proposed for NGC 3115. The first panel shows two high-redshift galaxies, one with a central black hole. The second panel shows the moment after the encounter between the two galaxies, forming the bulge (red ellipse), thick disc (light blue ellipse) and halo (light green circle), and, afterwards, the formation of the thin disc (blue equatorial component). The third panel shows the formation of the bar (mustard equatorial rectangle), igniting the AGN activity (yellow central ring). Companion dwarf galaxies (yellow circles) are accreted, increasing the thick disc and the halo of the galaxy. The fourth panel shows the effect of AGN feedback (yellow lines coming out of the bulge) and the quenched thin disc (light pink equatorial component). The fifth panel shows the nearly passive galaxy with low AGN activity. The sixth panel shows the passage of the dwarf companion KK 084 (left blue ellipse), the direction of displacement of the core due to the passage (black arrow), its possible trajectory (red arrow), and the effect of the passage: bulge becoming bluer in its centre (central blue ellipse) and formation of spiral-like features (blue semicircles in the middle of the thin disc). enhancement by an encounter with the companion KK084 that was capable, on one hand, to increase the star formation in the bulge and, on the other hand, to cause instabilities in the disc, generating spiral-like/ ring-like features in the faint end of the thin disc. This interaction would have caused a core displacement in NGC 3115, as discussed in Section 6.2.
The process of rejuvenating S0 galaxies has been largely studied in the past few years (Bresolin 2013;Mapelli M. et al. 2015;González Delgado et al. 2017;Barway & Saha 2020). In fact, Mapelli M. et al. (2015) showed, using simulations, that recent encounters of S0 galaxies with gas-rich satellites usually form gas-rings and can reignite star formation in the central region of galaxies, rejuvenating the S0. They show that minor mergers may trigger episodes of star formation in the S0s that would last for ∼10 Gyr, which could explain both the presence of the ring-like features and the recent star formation in the bulge of NGC 3115.
Incorporating these results into the previous scenarios proposed for its formation, we propose a novel scenario for the formation of NGC 3115: initially, a violent, dissipative gas-rich major merger created the bulge, halo, and thick disc of the galaxy, then the remaining gas from the merger was able to generate the observed thin disc. Instabilities in the disc would have created the central bar, and such bar would be responsible for dragging gas from the disc to the centre of the galaxy, feeding the SMBH and igniting the AGN activity. With time, the galaxy starts accreting companions, which, in turn, enrich and increase the size of the now-massive thick disc, resulting in a disc with a mixed population of accreted stars and stars formed in-situ. The continuous accretion and increasing mass of the galaxy cause feedback events, ceasing the star formation in the galaxy. The now gas-less thin disc is not able to feed the AGN any further, diminishing the AGN activity (observed to be low nowadays; Almeida et al. 2018;Jones et al. 2019) and the galaxy enters a phase of passivity. A recent encounter with the companion KK084 could have been capable of creating a starburst event in the centre of the galaxy, explaining the blue population and high SFR observed, and have also caused instabilities in the disc that would force the stars to rearrange in the spiral-like features observed today. Furthermore, the passage of the companion KK084 would have caused a core displacement to the south in NGC 3115 just a few hundred thousand years ago, and the direction of this displacement may set a constrain in the trajectory of the companion galaxy. This formation scenario is summarized in Fig. 12.

C O N C L U S I O N
In this work, we study the field lenticular galaxy closest to the Milky Way, NGC 3115. For this purpose, we use 11 images from UV to IR wavelengths, observed with different telescopes, to perform a multiwavelength broad-band fitting. Using GALFITM, we model this galaxy with a set of initial structural parameters to find the ones that best fit the galaxy and its components. The best-fitting model of NGC 3115 has three components: a bulge with Sérsic profile described by n = 3.5 and Re = 21.69 arcsec (0.97 kpc), a thick disc and a thin disc, whose disc scale length increase linearly with wavelength. Given the high value of the axial ratio of the thick disc (b/a = 0.44), it might be interpreted as an oblate spheroid. We also find evidence of the presence of a bar and spiral-like features from the residuals of these fits. We find that the age of the galaxy as a whole and its components are around 11 Gyr, in agreement with previous works.
In particular, we recover the same age/metallicity gradients as found by Guérou et al. (2016) and Poci et al. (2019) using IFU optical data. The results of this work are also consistent with the NUV-r results from Kaviraj et al. (2007) and the Sérsic index (n) versus BH mass relation found by Dullo & Graham (2014) and Graham & Driver (2007).
In the present analysis, by recovering the star formation history of the galaxy components, after performing SED fitting using CIGALE, we show, for the first time, that the galaxy may have had a last burst at a few Myr ago, during a pericentre passage of the companion KK084, resulting also in a core displacement of the lenticular to the South, where the direction of this displacement can be used to constrain the trajectory of the companion around the lenticular.
We propose a formation scenario for NGC 3115 that incorporates results from previous works (Arnold et al. 2011(Arnold et al. , 2014Cortesi et al. 2013;Peacock et al. 2015;Guérou et al. 2016;Zanatta et al. 2018;Poci et al. 2019;Dolfi et al. 2020) with new findings from this work, as explained at the end of the Section 6.3 and summarized in Fig. 12.
We conclude that multiband fitting of galaxy images followed by multicomponent SED fitting are powerful tools to recover galaxies' formation histories and physical properties, consistent with spectroscopic studies. In this sense, the novel multiband surveys S-PLUS, J-PLUS, and J-PAS provide a new window to the understanding of galaxy evolution, especially when complemented with archive allsky surveys (as GALEX and WISE).

AC K N OW L E D G E M E N T S
We deeply thank the referee and Boris Haussler for their insightful comments and suggestions on the paper.
MLB, AC and CMdeO acknowledge the financial support provided by the São Paulo Research Foundation (FAPESP; grants 2018/09165-6 and 2019/23388-0). AC acknowledge the financial support provided by CAPES.

DATA AVA I L A B I L I T Y
The data used in this work is available via the GALEX, Subaru Suprime Cam, DECam, 2MASS, and WISE mission archives.

A P P E N D I X A : G A L F I T M M O D E L S
The fundamental steps in fitting a galaxy and correctly decompose its light into its main components using GALFITM are: identifying the number of physically motivated components, identifying the order of the series associated with each parameter for every component, and defining the initial conditions of each parameter. The components of a galaxy are usually defined by some key parameters, such as Sérsic index (n), effective radius (Re), disc scale-length (Rs), position angle (PA), axial ratio (b/a), and diskyness/skewness (C0). Each of these parameters in GALFITM is expanded in a wavelength-dependent series, and the variation of these parameters with the wavelength is defined by the order of this series (orders). As previously explained, this order can vary from the most restrict case, i.e. orders = 0, where all bands are tied to the initial condition, to the freest case, i.e. orders = number of fitted images (11, in our case), where this parameter can vary freely between bands. The number of fitted parameters will increase with the number of components, and so determining the best combination of the order of the series of each parameter can be a hard task, such as manually combining these orders and evaluating each output can take a lot of time. In order to find this best combination, we created several GALFITM models using different combinations of orders of these parameters. This was done for models with one, two, and three components.
For the models with one component, since we are not decomposing the galaxy in subcomponents, and, therefore, not separating its different colours, one might expect that the series would need to vary with a larger order, to account for differences in the structural parameters of the galaxy. Therefore, for the models with one component, we let five parameters vary with the order spanning from 0 to 7, creating models with each combination of these orders. These parameters were n, Re, PA, b/a, and C0. This led to a total of 32 768 models.
For the models with two components, we were more restrictive, since we are now decomposing the galaxy in at least two colours, we expect that the variation would be smaller between bands, and, therefore, the orders were only allowed to span from 0 to 3. Also, we can assume that there is no significant change among bands in the position angle or axial ratio of these components, therefore, we fix the order of these parameters to be 1, i.e. constant between bands, for both components. Moreover, analysing the models with one component, we realized that the C0 parameter was not affecting our results, so we excluded this parameter from the further modelling. We now have three parameters: n, Re, and Rs, and their order can vary from 0 to 3, leading to 64 models.
Finally, for the model with three components, applying the same cuts as we did for the model with two components, we have four parameters that we allow to vary their order spanning from 0 to 3, these are n, Re, Rs-disc1 (hereafter Rs1) and Rs-disc2 (hereafter Rs2), which led to 256 models.
The models created and the orders of the series that each parameter can vary with are summarized in Table A1.
These models were first analysed taking into consideration three statistics that we find representative of the goodness of the fit, these are: (i) Model/Input: this value is obtained by dividing the model images by the input images and taking the median of the resulting image. The concept behind this is to understand how well the model is mimicking the input image. Consequently, the best models are those who have Model/Input closest to one.
(ii) Residual/Input: this is obtained by dividing the residual image by the input image in each band and taking the median of the result. With these statistics, we can understand how significant the residual Table A1. Order of the series that each parameter can vary to find the best-fitting model with GALFITM.
fluxes are with respect to the data. The best models are expected to have values of Residual/Input close to zero, where, if this is smaller than zero, means that the galaxy is being overfitted, and if higher, underfitted.
(iii) Residual/Sigma: these final statistics are obtained by dividing the residual image by the error of the input data. The idea is to understand how significant the residuals are concerning the error. Therefore, if sigma >> residuals, the residuals are well within the error and this is a good fit. If residual > sigma, the residuals are bigger than the errors and this is not considered such a good fit. Therefore, the best models will have the smallest Residual/Sigma.
We also analyse our models in terms of the Bayesian Information Criteria and Akaike information criteria. However, these methods are taken less into account, since they are strongly biased by the best χ 2 , which not necessary resembles the best-fitting model, since we are looking for models with physically meaningful components, and sometimes the models with the lowest χ 2 are found by creating mathematical artefacts. In Fig. A1, we show the distribution of the models with one, two, and three components in the 3D plane, shaped by the Model/Input, Residual/Input, and Residual/Sigma statistics.
From the left-hand side of Fig. A1, we can directly understand that models with one component are not sufficient to recover the total light of NGC 3115, as they fall substantially lower than the models with two and three components in the Model/Input axis, meaning that the model found are comprising about half of the light of the input images. In the right-hand side of Fig. A1, we focus on the models with two and three components. As we can see, the majority of these models share a locus in the plot. Both models with two and three components behave quite well in the three statistics. We can see, however, that there are about five models with three components that outstand this locus, presenting values close -or exactly the expected -for the three statistics, i.e. Model/Input ≈ 1, Residual/Input ≈ 0, and Residual/Sigma ≈ 0. To further investigate the difference between the models with two and three components, we analyse in Figure A1. GALFITM models created for NGC 3115 using one, two, and three components, with different combinations of orders that define how each parameter can vary with the wavelength. Models are evaluated according to the statistics: Residual/Input (how the residual compares with the input image: the smallest the better), Residual/Sigma (how the residual compares to the error of the input image, the smallest the better), and Model/Input (how well the model is mimicking the input image: the closest to one the better). The bottom line shows the 2D projection of each pair of axes of the 3D plot for the models with one, two, and three components. Fig. A2, the BIC and AIC of these models. Confirming what we found previously, what we can see in Fig. A2 is that the models with three components (Bulge + Thin disc + Thick disc) deliver the minimal value of BIC and AIC, representing the best fits for the data. Now, to better understand the variation of these three component models in the plane, we show in Fig. A3, the variation of the models according to the order of the series of each parameter (n, Re, Rs1, and Rs2).
There are insightful conclusions that we can take from Fig. A3. First, analysing the first panel, we can see that the five best-fitting models are the ones where the Sérsic index can either be fixed to the initial condition (orders = 0) or constant among bands (orders = 1), showing that, when we carefully separate the galaxy into components, the parameters that define each component do not vary significantly between bands, because we are most likely not dealing with mixed colours/structural parameters anymore. The radius, however, as seen from the other panels in Fig. A3, can still have some small variation (i.e. orders = 1, constant or orders = 2, linear) among bands, as it would be expected, since we know from observing the galaxy, that its radius severely changes from UV to IR. Most importantly, from all panels, what we see is that large freedom of the parameters seems to worsen the models, reinforcing once again the power of GALFITM to tie the images together and deliver reliable and consistent physical parameters.
From this analysis, we can take several important information that the best model of NGC 3115 must comprise: (i) Morphology: three components (a bulge and two exponential discs).
(ii) Orders of the parameters: n, PA, and b/a should either be fixed to the input value or constant among bands. Re, Rs1, and Rs2 can either be constant or vary linearly, but not more than that.
(iii) The parameters recovered in the best-fitting model in this simulation can be used as reliable initial conditions to future models.
Although the best-fitting model found by these simulations is not the one chosen as the best for this work, with this analysis, we recovered the three fundamental steps, as described in the first paragraph of this appendix, to create the most consistent and robust GALFTIM model, described in the main body of this paper and used for our analysis. Figure A2. BIC and AIC of the models with two and three components, where the different models are created by using different combinations of orders of the wavelength-dependent series that define the variation of the parameters among bands. Figure A3. GALFITM models created for NGC 3115 using three components and different combinations of orders of the parameters. Models are evaluated according to the statistics: Residual/Input (how the residual compares with the input image: the closest to zero the better), Residual/Sigma (how the residual compares to the error of the input image, the closest to zero the better), and Model/Input (how well the model is mimicking the input image: the closest to one the better). Upper left-hand panel: distribution of models coloured with the orders that define the variation of the Sérsic index among bands. Upper right-hand panel: distribution of models coloured with the orders that define the variation of the effective radius among bands. Middle left-hand panel: distribution of models coloured with the orders that define the variation of the disc scale-length of the first disc among bands. Middle right-hand panel: distribution of models coloured with the orders that define the variation of the disc scale-length of the second disc among bands. Bottom line: the 2D projection of each pair of axes of the 3D plot. These panels are coloured with the orders of variation of the effective radius, such as the top left-hand panel. This paper has been typeset from a T E X/L A T E X file prepared by the author.