Combining Stellar Populations with Orbit-Superposition Dynamical Modelling - the Formation History of the Lenticular Galaxy NGC 3115

TL;DR

This paper introduces a novel method combining dynamical modeling and stellar population data to unravel the formation history of the galaxy NGC 3115, revealing distinct structural components and their evolutionary paths.

Contribution

The study develops an innovative approach integrating orbit-superposition models with stellar population maps, enabling detailed analysis of galaxy formation processes.

Findings

Identification of a metal-rich, dynamically-cold disk component.

Detection of an old, diffuse stellar halo with varying kinematics.

Correlation between dynamical support and stellar age similar to the Milky Way.

Abstract

We present a combination of the Schwarzschild orbit-superposition dynamical modelling technique with the spatially-resolved mean stellar age and metallicity maps to uncover the formation history of galaxies. We apply this new approach to a remarkable 5-pointing mosaic of VLT/MUSE observations obtained by Gu\'erou et al. (2016) extending to a maximum galactocentric distance of $120''$ ($\sim 5.6\ {\rm kpc}$) along the major axis, corresponding to $\sim 2.5\ {\rm R}_e$. Our method first identifies `families' of orbits from the dynamical model that represent dynamically-distinct structures of the galaxy. Individual ages and metallicities of these components are then fit for using the stellar-population information. Our results highlight components of the galaxy that are distinct in the combined stellar dynamics/populations space, which implies distinct formation paths. We find evidence for a dynamically-cold, metal-rich disk, consistent with a gradual in-situ formation. This disk is embedded in a generally-old population of stars, with kinematics ranging from dispersion-dominated in the centre to an old, diffuse, metal-poor stellar halo at the extremities. We find also a direct correlation between the dominant dynamical support of these components, and their associated age, akin to the relation observed in the Milky Way. This approach not only provides a powerful model for inferring the formation history of external galaxies, but also paves the way to a complete population-dynamical model.

Figures

Equations

$$M_{\mathrm{DM}}(r)$$

$$\Phi (r)=-\frac{G{M}_{\bullet }}{\sqrt{r^2+a^2}}$$

$${\chi }_r^2$$

$${\chi }_{\mathrm{kin}}^2$$

$$N_{\mathrm{kin}}$$

$$=6\times 4881$$

$$=29286$$

$${\lambda }_z$$

$$\overline{L_z}$$

$$\overline{r}$$

$${\overline{V_2}}^2$$

$${\Phi }^{\mathrm{SFH}}$$

$${\tilde{\omega }}^j$$

$$\left(\begin{array}{cccc}{\tilde{\omega }}_1^1& {\tilde{\omega }}_1^2& \cdots & {\tilde{\omega }}_1^{N_{comp.}}\\ {\tilde{\omega }}_2^1& {\tilde{\omega }}_2^2& \cdots & {\tilde{\omega }}_2^{N_{comp.}}\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{\omega }}_{N_{\mathrm{aper}}}^1& {\tilde{\omega }}_{N_{\mathrm{aper}}}^2& \cdots & {\tilde{\omega }}_{N_{\mathrm{aper}}}^{N_{comp.}}\end{array}\right)\cdot \left(\begin{array}{c}{\varphi }^1\\ {\varphi }^2\\ ⋮\\ {\varphi }^{N_{comp.}}\end{array}\right)$$

$$⟨V_k^2⟩$$

$${\sigma }_k$$

$$\int {\left[{\vec{\varphi }}^{\prime \prime }(x)\right]}^2\mathrm{d}x$$

Full text

††pubyear: 2019††pagerange: Combining Stellar Populations with Orbit-Superposition Dynamical Modelling - the Formation History of the Lenticular Galaxy NGC 3115–C

1 Introduction

The present-day observed state of a galaxy is the result of the integration over its entire formation history, including external accretion/mergers, in-situ star-formation, and passive stellar evolution. To determine how and when a galaxy has built up its mass, it is necessary to disentangle its present-day state into spatially- and chemically-distinct events. Typically, the studies of stellar populations and dynamical properties have remained independent, however it is the union of these two aspects that is necessary to be able to investigate the origin of a galaxy’s mass over its formation history. These ideas have been investigated for some time for the Milky Way, where a wealth of chemical and dynamical information can be directly obtained for individual stars. We endeavour here to extend these analyses to external, unresolved galaxies.

Historically, dynamical models of galaxies have been utilised for a wide range of applications, from individual galaxy analyses to large statistical samples of the galaxy population. These efforts have attempted to place constraints on, among other properties, the mass of the super-massive black holes (SMBH) at the centres of external galaxies (see Kormendy & Ho, 2013 and references in Table 1 therein, and more recently Seth et al., 2014; Krajnović et al., 2018); the IMF mass normalisation (for instance, see Thomas et al., 2011; Cappellari et al., 2012; Lyubenova et al., 2016; Davis & McDermid, 2017; Li et al., 2017; Oldham & Auger, 2018); the dark-matter content/distribution in galaxies using stellar kinematics (for example, see Cappellari et al., 2013; Tortora et al., 2019) all the way down in mass to dwarf spheroidals (Jardel & Gebhardt, 2012), and using gas kinematics (for example, see Corbelli & Salucci, 2000; Gentile et al., 2004; Di Teodoro et al., 2019). Dynamical models have also been used to uncover new relationships between galaxy structural parameters, including the widely-used stellar spin-ellipticity $(\lambda_{r}-\epsilon)$ correlation (introduced in Emsellem et al., 2007), as well as other observed correlations with dynamical properties (Cappellari, 2016). Finally, dynamical models of individual galaxies have been used to probe internal dynamical structures in great detail (for instance, see Krajnović et al., 2005; van den Bosch et al., 2008; Krajnović et al., 2015; Zhu et al., 2016a), placing strong constraints on the formation mechanisms of those specific galaxies.

Independently, stellar-population models of galaxies have been utilised for their own range of applications, covering similar scales in sample size. These models have been applied to constrain the IMF shape, normalisation, and low-mass cut-off (for instance, see Conroy & van Dokkum, 2012; Smith et al., 2012; Martín-Navarro et al., 2015; Alton et al., 2017; Rosani et al., 2018; Dries et al., 2018; Vaughan et al., 2018); the measurement and interpretation of $\alpha$-enhancement (as in Thomas et al., 2005; Conroy et al., 2014; Greene et al., 2015; McDermid et al., 2015); the measurement of spatial gradients in stellar properties (for instance, see Mehlert et al., 2003; Sánchez-Blázquez et al., 2007; Kuntschner et al., 2010; Li et al., 2018); assembly timescales for galaxy formation (Martin-Navarro et al., 2018). Similarly to dynamical models, stellar population models have been used to uncover new relationships between galaxy structural parameters through a variety of scaling relations (see Gallazzi et al., 2008; Graves et al., 2009a, b; McDermid et al., 2015; Walcher et al., 2015, for some such relations). Again, stellar population models of individual (or handfuls of) galaxies in great detail offer insight into the specific formation path of these galaxies, as well as the presence of sub-structures with distinct stellar populations (for instance, see Barbosa et al., 2016; Mentz et al., 2016; Streich et al., 2016).

Driven by the influx of spatially-resolved observations coming from integral-field units (IFU), more recent investigations have focused on attempting to infer structural, dynamical and/or chemical properties for localised regions of galaxies, by decomposing them into physically-motivated components. In fact, this concept pre-dates the large-scale use of IFU with works dealing with, for instance, photometric disk/bulge decompositions of surface brightness profiles (for example, see Kent, 1985; Cinzano & van der Marel, 1993; Scorza & Bender, 1995; Moriondo et al., 1998; Krajnović et al., 2013), including using multiple filters (Dimauro et al., 2018). Aside from surface brightness, radial profiles of other parameters have also been subject to analogous decompositions. For instance, the decomposition of gas (usually $\nuclide{H}\textsc{i}$) circular velocity profiles into contributions from different galaxy sub-components is a well-established practise (as in van Albada et al., 1985; Carignan et al., 1988; Battaglia et al., 2006; Noordermeer et al., 2007; Swaters et al., 2012; Sofue, 2017; Aniyan et al., 2016; Aniyan et al., 2018), while decompositions of mass profiles have attempted to infer the contributions from dark matter and baryons (stars and globular clusters, gas, et cetera; Poci et al., 2017; Annunziatella et al., 2017; Bellstedt et al., 2018). These concepts have been extended to two dimensions, including multi-band photometric disk/bulge decompositions of images, rather than profiles (for instance, see Scorza et al., 1998; de Souza et al., 2004; Norris et al., 2006; Simard et al., 2011; Méndez-Abreu et al., 2017; Dalla Bontà et al., 2018). Moreover, there have been recent efforts to conduct the decomposition directly on an observed spectrum (Johnston et al., 2012; Coccato et al., 2015; Tabor et al., 2017; Coccato et al., 2018), to similarly determine the contributions to various spectral features coming from ‘distinct’ galaxy subcomponents. This type of component-based approach attempts to isolate the distinct contributions to observed galaxies from regions which may or may not have had different origins and/or formation paths, however they have thus far dealt with the problem from only one perspective - dynamics or stellar populations.

These works, and others, have motivated a need for combining the aforementioned dynamical and chemical models in order to investigate the full formation history of a galaxy. This combination is in fact necessary to be able to determine the origin of the different components of a galaxy. There have been only a small number of works in which different chemical and dynamical populations are simultaneously treated, such as the models of Zhu et al. (2016b, a), which use globular clusters (GC) as discrete tracers of the kinematics, while simultaneously fitting for two (chemical) populations of GC (based on the discrete modelling prescription of Watkins et al., 2013). This showed that there are observable distinctions to be made in the combined populations/dynamics space (and indeed physical space) between different components of galaxies.

In this paper, we describe a new approach which aims to decompose a galaxy into dynamical and chemical components in order to infer its formation history. We present also an application of this method to a real galaxy; the nearby lenticular galaxy NGC 3115. Section 2 describes the observational data used in this work. Sections 3 and 4 describe in detail the multiple aspects of our chemical and dynamical analyses, respectively. The results of our show case are presented in Section 5. Section 6 presents the formation history of NGC 3115 as determined from our method, and the connections between dynamical and chemical properties.

2 Observational Data and Kinematic Extraction

The observations used in this work were obtained and reduced by Guérou et al. (2016). NGC 3115 is the nearest lenticular (${\rm S}0$) galaxy to the Milky Way, with an orientation very close to edge-on. Tonry et al. (2001) used surface brightness fluctuations to measure a distance modulus for NGC 3115 of $(m-M)=29.93\pm 0.09$, placing it at a physical distance of $9.7~{}$\mathrm{M}\mathrm{p}\mathrm{c}. It has an effective (half-light) radius of $R_{e}=47.32$\mathrm{\SIUnitSymbolArcsecond}$\big{/}2.23~{}$\mathrm{k}\mathrm{p}\mathrm{c} (Emsellem et al., 1999). Our data extends to a galactocentric radius of $\sim$120\text{,}\mathrm{\SIUnitSymbolArcsecond}\left(5.6\text{,}\mathrm{k}\mathrm{p}\mathrm{c}\right) along the major axis, and so we have coverage out to $\sim 2.5R_{e}$.

The data set has a pixel-scale of $$0.2\text{,}\mathrm{\SIUnitSymbolArcsecond}\ \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}^{-1}$\ ($9.4\text{,}\mathrm{p}\mathrm{c}\ \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}^{-1}$)$, and consists of over $360,000$ individual spectra. We refer the interested reader to Guérou et al. (2016) for further details regarding the observational procedure and data-reduction techniques.

For all subsequent analysis, we consider only the spectral range $\lambda\in\left[4760,6400\right]~{}$\mathrm{\SIUnitSymbolAngstrom}$$. This is to reduce the impact of residual sky emission lines on the extracted stellar kinematic and population properties. We spatially bin the reduced datacube using a Python implementation111Available at http://www-astro.physics.ox.ac.uk/~mxc/software/ of the Voronoi tessellation algorithm (Cappellari & Copin, 2003) to a target signal-to-noise ratio $(S/N)$ of 90 per bin. This relatively high threshold is set to ensure an accurate recovery of (the moments of) the line-of-sight velocity distribution (LOSVD), as well as the subsequent stellar population analyses.

Although kinematics were extracted by Guérou et al. (2016), we re-derive the kinematics222These data products are also available at https://datacentral.org.au/teamdata/NGC3115/public/ here with the new $S/N$ threshold for the spatial binning. This allows us to extract the first six Gauss-Hermite coefficients of the parametrised LOSVD, which provide important additional constraints on the dynamical model.

For the kinematic extraction, we use a Python implementation1 of the parametric Penalised PiXel-Fitting code ppxf (Cappellari & Emsellem, 2004; Cappellari, 2016), with 985 stellar templates from the empirical MILES library (Sánchez-Blázquez et al., 2006; Falcón-Barroso et al., 2011). ppxf finds the linear combination of the provided templates that, when convolved with an LOSVD parametrised by the 6 Gauss-Hermite coefficients, best matches the observed spectrum. Computing such a fit for every (Voronoi-binned) spectrum in the cube provides the best-fitting LOSVD (or parametrisation thereof) at every spatial location. Since the spectral resolutions of the MILES models and MUSE data are comparable, both sets of spectra are kept at their native resolution. Every bin is fit with the freedom of the full template library, and with an additive polynomial of order 16. This combination of additive polynomial and stellar templates is used to achieve the best possible fit to the spectrum, in order to most accurately recover the LOSVD, without being tied to any particular stellar population model. During the fitting, spurious spectral pixels / artefacts in the data are iteratively clipped. Errors for all 6 Gauss-Hermite moments are computed by Monte Carlo simulations on every spectrum individually for which kinematics are extracted. We derive a mean uncertainty (and apply a floor on these uncertainties during the dynamical model described in Section 4), for $V$ and $\sigma$ of $10.30\ (7.00)$ and $3.17\ (2.00)\$\mathrm{km}\text{,}{\mathrm{s}}^{-1}$$, respectively, and for moments $3$ through $6$ of $0.0016\ (0.0016)$, $0.0020\ (0.0019)$, $0.0016\ (0.0016)$, $0.0021\ (0.0020)$. The floor on the uncertainties prevents the dynamical model from preferentially fitting the central pixels (which have the smallest errors).

3 Spatially-Resolved Star-Formation History

Our analysis of NGC 3115 requires the computation of a spatially-resolved star-formation history (SFH). To do this, we employ full-spectral-fitting techniques, with the aim of investigating chemically-distinct components in our dataset. To fit the spectra, we again utilise ppxf, but use the MILES-IndoUS-CaT (MIUSCAT) Single Stellar Population (SSP) model library (Vazdekis et al., 2010) as templates. While we aim to compare our results with those from Guérou et al. (2016), we adopt different SSP templates for our SFH. From an initial run of the dynamical model, we inferred a dynamical (enclosed) mass that was lower than the stellar mass derived assuming a Salpeter (1955) IMF (which was the assumption in Guérou et al., 2016). Such seemingly non-physical discrepancies have been found previously (Lyubenova et al., 2016), and is interpreted as evidence to exclude single-power-law IMF shapes. We therefore assume a Kroupa (2002) IMF for all stellar-population models in this work. The SSP library used here is based on the Padova isochrones (Girardi et al., 2000), with ages $0.1\leq t\leq 14.1\$\mathrm{G}\mathrm{y}\mathrm{r} and metallicities $-1.71\leq Z\leq 0.22\ $\mathrm{d}\mathrm{e}\mathrm{x}333Determined from the “safe ranges” taken from the MILES website over the MUSE wavelength range. These models do not consider any $\alpha$ abundances explicitly, and since they are based on empirical stars, therefore share the same $\alpha$-enhancement characteristics as the solar neighbourhood. During the ppxf fit, we employ a small fractional regularisation to reduce the intrinsic degeneracy of the fit. This regularisation prefers a smooth solution that would otherwise be degenerate with a ‘spiky’ SFH. We do not use any additive polynomial, as this would change the relative strengths of the spectral features, which would in turn impact the derived stellar population properties. We do, however, employ a multiplicative polynomial of order 16 in order to account for the continuum. The LOSVD are left free during the stellar-population fits. This is done to minimise any possible systematics and template mismatch, while maintaining good fits across the FOV. We confirmed that the resulting SFH2 is consistent within the errors when the kinematics were fixed to those extracted in Section 2.

Fig. 1 illustrates this fitting concept on the highest and lowest $S/N$ spectra in the dataset, as well as the associated SSP weights for the ages, $t$, and metallicities, $Z$.

We compute similar fits for every binned spectrum in the datacube to investigate the spatial behaviour of the stellar populations. For reasons discussed in Section 4.2, we compute luminosity-weighted SSP properties by removing the relative normalisation of all the template spectra prior to fitting.

4 Dynamical Model

In order to identify intrinsically-distinct components, we employ the Schwarzschild orbit-superposition technique (Schwarzschild, 1979). In this work we use a triaxial extension of the axisymmetric method originally presented in van der Marel et al. (1998) and Cretton et al. (1999), and developed further by Cappellari et al. (2006). This triaxial implementation is detailed extensively in van de Ven et al. (2008); van den Bosch et al. (2008) and developed further by Zhu et al. (2018a, b). We present here a brief summary of the relevant aspects of this implementation, and refer the reader to the above references for further details.

4.1 Mass Model

The Schwarzschild orbit-superposition method generates galaxy models within a given stationary gravitational potential. Fitting Schwarzschild models to observations therefore requires finding the best input gravitational potential that reproduces the observable constraints. Since the gravitational potential can not be measured directly, we construct each input mass model as the combined contributions from a stellar mass model, a dark matter parametrisation, and a point-source central super-massive black-hole (SMBH) component.

4.1.1 Stellar-Mass Model

The stellar mass is a derived quantity, dependent on a number of assumptions relating the directly-observed stellar light to an implied stellar mass. This in turn requires a model of the observed surface brightness from which a mass can be derived. In this work, we use the surface-brightness model of NGC 3115 presented in Emsellem et al. (1999), which used the multi-Gaussian Expansion (MGE) technique (Monnet et al., 1992; Emsellem et al., 1994). MGE models fit a series of $2{\rm D}$ Gaussians to the observed photometric isophotes. One advantage of using the MGE approach is that for any inclination that is trialled by the dynamical model, the Gaussians have an analytic deprojection into an intrinsic $3{\rm D}$ model. This results in a fast (though not necessarily unique) description of the mass for that inclination, forming the framework within which the dynamical model is computed. Emsellem et al. (1999) used a combination of high-resolution HST and ground-based photometry to compute a photometric fit out to $\sim 300$\mathrm{\SIUnitSymbolArcsecond}$$, which is necessary in order to ensure that the stellar mass model comfortably encloses the extent of such a nearby galaxy. The collapsed data cube and MGE surface brightness model are shown in Fig. 2.

Typical Schwarzschild model implementations, and many dynamical models in general, assume that the conversion from light to mass can be done with a single global scale for a given galaxy; that is, a spatially-constant stellar mass-to-light ratio ($M_{\star}/L$). This implies that the observed light originates from only a single population of stars. In our work, however, we consider the interplay between dynamical and stellar-population structures, and so the spatial variations in the $M_{\star}/L$ are of particular importance. Since we have already characterised the presence of multiple stellar populations within this galaxy in Section 3, it is therefore possible to incorporate the spatial structures in the stellar populations to obtain a more accurate conversion from light to mass. It is imperative that we attempt to construct an accurate input mass model that takes into account this information in order to maintain self-consistency when we analyse the outputs. There have already been recent implementations of this for other purely-dynamical modelling techniques by parametrising the derived $M_{\star}/L$ map into a $1{\rm D}$ radial profile, which is then used to scale the surface-brightness model accordingly (Poci et al., 2017; Mitzkus et al., 2017). We compute a $V$-band $M_{\star}/L_{V}$ map for NGC 3115 based on the ppxf fits from Section 3, shown in Fig. 3. $V$-band is consistent with both the spectral range used in the MUSE observations and the original photometry used by Emsellem et al. (1999).

It is immediately clear, however, that the $M_{\star}/L_{V}$ map of NGC 3115 can not be well-approximated by a radial profile due to the ‘lobe’ features along the major axis created by a relatively young stellar disk. In order to derive an accurate stellar-mass model, we take a different approach here in order to maintain the $2{\rm D}$ information from the SFH. We scale the surface-brightness MGE (MGEμ) by the MUSE $M_{\star}/L_{V}$ map directly in order to obtain a mass ‘image’, to which we fit a new mass density MGE (MGEΣ).

One issue when comparing photometric and IFU observations is the difference in the size of the FOV. To overcome this, we first evaluate MGEμ on an image grid that is sufficiently large (comparable to the FOV used in Emsellem et al., 1999). We then cast the $M_{\star}/L_{V}$ map onto the same image scale. To populate the pixels that lie outside of the MUSE FOV, we assume that the $M_{\star}/L_{V}$ is constant at large radius. This already appears to be the case in the MUSE observations for $R\gtrsim 80$\mathrm{\SIUnitSymbolArcsecond}$$, and so the exterior pixels are fixed to the average value of the outermost bins of Fig. 3. Fig. 4 presents MGEμ and MGEΣ on the approximate scale of the stellar disk, which clearly highlights the impact of considering the $2{\rm D}$ $M_{\star}/L_{V}$ map.

The MGEΣ has rounder iso-mass contours (compared to MGEμ) in the region where relatively young stars contribute a lot of luminosity but not much mass. Assuming a constant $M_{\star}/L$ would have attributed too much mass to the disk region (along the major axis), which in turn would have biased the dynamical model and the inferred stellar populations. MGEΣ is tabulated in Appendix A2.

4.1.2 Dark-Mass Model

To include contributions from non-luminous mass, we include a parametrisation of the dark matter (DM) halo by assuming it has the form of a generalised spherical Navarro, Frenk, and White (NFW; Navarro et al., 1996) halo, as described in Eq. (3) of Zhu et al. (2018b), but included here for completeness in Eq. 1.

[TABLE]

where $g\left(C_{\rm DM}\right)=\left[\ln\left(1+C_{\rm DM}\right)-\frac{C_{\rm DM}}{1+C_{\rm DM}}\right]^{-1}$, and $M_{200}=200\times\frac{4}{3}\pi\rho_{c}r_{200}^{3}$. Here, $\rho_{c}=\frac{3H^{2}}{8\pi G}$ is the critical density of the Universe with Hubble constant $H$ and gravitational constant $G$. While we do not test different dark matter prescriptions in this work, we note that the broader method described in this work is not tied to the specific parametrisation of the dark matter halo.

The central SMBH similarly contributes non-luminous mass to the mass model, and therefore affects which orbits exist within the gravitational potential. Our models include a spatially-localised non-luminous Plummer potential (Dejonghe, 1987) to represent the mass of the SMBH, $M_{\bullet}$, defined as

[TABLE]

where $a$ is the Plummer core radius (effectively the SMBH softening length). In the region $r\leq a$, the gravitational potential of the model is dominated by the Plummer potential. We fix $a=$0.008\text{,}\mathrm{\SIUnitSymbolArcsecond}$$, which is below the MUSE resolution, in order to avoid numerical issues when integrating orbits close to the SMBH.

4.2 Orbit Solution

The Schwarzschild code generates a library of orbits that are physically permitted for the given mass model. In our model, we sample the orbits over 42 logarithmically-spaced radial starting locations, which conserve the first integral of motion; the energy. There are also 15 linearly-spaced locations in $I_{2}$, the second conserved integral of motion. Finally, we sample the third (non-classical) conserved integral of motion, $I_{3}$, from 12 linearly-spaced starting locations. We may cheaply (without re-integrating) increase the size of the orbit library, by simply mirroring orbits about their integrals of motion to allow for counter-rotation. In the general triaxial code, this occurs only for those orbits which are sampled on the $(x,z)$ meridional plane, doubling the number of those orbits. Moreover, to boost the number of box orbits accessible to the model, we launch a new set of box-only orbits, adding another factor to the orbit library (van den Bosch et al., 2008). Therefore, the full number of orbits in the model is $N_{\rm orb.}=3\times N_{E}\times N_{I_{2}}\times N_{I_{3}}=3\times 42\times 15\times 12=22680$. In order to minimise the discreteness of the orbits, each starting position is ‘dithered’ (see Cappellari et al., 2006) in the three-dimensional integral space by a factor of 5, thereby sub-sampling each $(E,I_{2},I_{3})$ location into a grid of $5^{3}$ adjacent starting positions. This high sampling ensures that the small physical scale of the observation is probed sufficiently, while also covering a large enough volume to describe the full extent of the galaxy. Compared to the extensive work done on CALIFA data (Zhu et al., 2018b), our data has $5$ times higher spatial resolution, and our orbit library is slightly over $5$ times larger. Specifically, we can estimate the approximate density of orbits across our FOV. By taking the conservative approach of considering the orbital starting locations, we find an average of $\sim 102\$\mathrm{o}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{2}$$. From this density of orbits, we conclude the orbit library contains sufficient resolution to accommodate our data. Each of these orbits is integrated for 200 complete periods.

With a complete orbit library, the Schwarzschild code then solves for the weighted linear combination of these orbits that best reproduces the (projected) observations. The orbital weights are computed by a Non-Negative Least-Squares (NNLS) fit to all observed constraints simultaneously, producing a weight per orbit, per spatial bin. For a given combination of orbits, the combined model is projected and convolved with the PSF to mimic the observations, binned onto the same Voronoi bins, and finally, the intrinsic LOSVD of the model are compared to the measured kinematic moments in each spatial bin. It is clear here that each kinematic moment extracted from the observed spectra therefore provides additional constraining power, in order to discriminate between models with similar lower-order kinematic moments.

While the orbit integration is carried out in the gravitational potential (mass model), the orbit solution is fit to the observed luminosity-weighted properties - namely, the surface brightness MGE, and all 6 kinematic moments. The linear combination of orbits is therefore luminosity-weighted also, consistent with the observations and the SFH described in Section 3.

4.3 Model Free Parameters

The process described in Section 4.2 is for the single fixed gravitational potential within which the orbits were integrated. While the projected mass model is fixed by MGEΣ, the intrinsic galaxy shape that gives rise to this projection remains unknown. It is therefore necessary to optimise the free parameters of the intrinsic gravitational potential in order to obtain not only a best-fitting set of orbits for a given potential, but also the best-fitting potential itself. The parameters of the gravitational potential are described here.

For a triaxial system, the baryonic component can be described by three intrinsic shape parameters, $(q,p,u)$. These intrinsic shapes are translated into projected viewing angles according to Eqs. (7)$-$(9) of van den Bosch et al. (2008). While we fit for $q$ and $p$, we fix the third intrinsic shape parameter to $u=0.99999$. This imposes the assumption that the shape of the baryonic component of the galaxy does not depend on the azimuthal angle at which we are observing it (that is, invariant under rotations about the rotation axis). This is a reasonable assumption in our case, where we are modelling a nearly-edge-on, approximately axisymmetric galaxy.

There are also free parameters of the dark mass contributions to the gravitational potential. From Eq. 1, our model gains two free parameters; the concentration of dark matter $C_{\rm DM}$, and the fraction of dark-to-stellar mass within $r_{200}$, f_{\rm DM}\left(r_{200}\right)=M_{200}\big{/}M_{\star}. We also vary $M_{\bullet}$ in order to accurately probe the orbits in the central region of the galaxy.

The final free parameter of the Schwarzschild model is the $M_{\rm dyn.}/L$ scaling, which we refer to as $\Upsilon$ to avoid confusion with the SFH-derived $M_{\star}/L_{V}$. Although we have incorporated the shape of the stellar component into the gravitational potential using the MUSE $M_{\star}/L_{V}$ map, this final free parameter is a global (spatially-constant) scale of the total mass. It allows the model to produce a deeper or more shallow potential, as required by the dynamics, due to possible systematics in the derived $M_{\star}/L_{V}$ including the assumption about the IMF, as well as the accuracy of the dark matter parametrisation. Since $\Upsilon$ is a global scaling, the shape of the gravitational potential does not change. Therefore, this parameter can be efficiently optimised by scaling all of the orbital velocities by $\sqrt{\Upsilon}$, and refitting the model. A value of $\Upsilon=1$ implies that the input gravitational potential, in terms of the the absolute stellar $M_{\star}/L$, assumed $M_{\bullet}$, and assumed DM parametrisation, contains the correct enclosed mass to reproduce the observed kinematics.

We are thus left with 6 free parameters for the Schwarzschild model, described in Table 1.

To find the global best-fitting model, we employ the following $\chi^{2}$ grid search. From a starting guess for each parameter, we consider those locations which are $\pm 5\Delta$ (the step size defined in Table 1) away. This is to ensure that a sufficient volume of the parameter-space is probed to avoid local minima. Including the central position, there are then 3 initial trials for each free parameter. In total, for $N$ gravitational free parameters (excluding $\Upsilon$), the model begins with $3^{N}$ trials. Each of these $3^{5}=243$ initial gravitational-potential models is evaluated for 3 values of $\Upsilon$. Once these have completed, the best-fitting location becomes the new ‘centre’, and each free parameter is explored in increments of $\pm\Delta$ in a similar manner. This is repeated until a location is found such that all surrounding models in the $N{\rm D}$ space produce worse fits to the data, which is followed by one final run with $\pm 0.5\Delta$ in order to accurately characterise the best-fitting region.

The metric for model comparison used in this work follows the formalism of Zhu et al. (2018b). That work introduced a normalised $\chi^{2}$ metric in order to account for spatial pixels which are not truly independent. It is defined as

[TABLE]

for the number of Gauss-Hermite coefficients fit in the Schwarzschild model, $N_{\rm GH}$, and the number of spatial Voronoi apertures, $N_{\rm aper}$. In this way, a standard deviation of this $\chi_{r}^{2}$ distribution is given by $\sqrt{2N_{\rm kin}}$, and the grid search for NGC 3115 is illustrated in Fig. 5.

The corresponding best-fitting Schwarzschild model is shown in Fig. 6 for all seven projected constraints.

4.4 Intrinsic Angular Momentum

Following the optimisation of the gravitational potential parameters, the best-fitting Schwarzschild model provides intrinsic information about the galaxy, instead of being typically limited to projected quantities. The orbits of the model have intrinsic $3{\rm{D}}$ angular momentum and orbital anisotropy, and the model has true LOSVD rather than parametrisations of these distributions. These intrinsic properties are a key requirement for the application of our methodology.

One such property that we utilise directly here is the intrinsic angular momentum. Each orbit has a corresponding angular momentum vector $\vec{L}=(L_{x},L_{y},L_{z})$. We consider a scaled version of this intrinsic angular momentum, the ‘circularity’, which was introduced by Zhu et al. (2018b) and is defined as

[TABLE]

We use the bar notation to denote that these values are averaged over all of the $(x,y,z,V_{x},V_{y},V_{z})$ points along a given orbit’s integrated path. The probability distribution of the $\lambda_{z}-r$ plane for the best-fitting Schwarzschild model is shown in Fig. 7 (see Section 4.5.1). It shows only the region constrained by the kinematics, even though the model orbits can extend beyond this in order to fit the full stellar MGE at large radius.

Circularity is a projection of the full $3{\rm D}$ orbital distribution, and ‘clumps’ in this projection identify families of orbits with similar orbital characteristics in the meridional plane. This definition separates circular and box/radial orbits, having circularity values of $\left|\lambda_{z}\right|=1$ and $\lambda_{z}=0$, respectively. It is analogous to the ‘stellar-spin’ parameter (defined in Emsellem et al., 2007) in that it is a metric to discriminate between rotation- and dispersion-supported systems (or orbits).

4.5 Applying Stellar Populations to the Orbital Structures

At this point, we have a dynamical model, fitted to the observed kinematics, that provides information on the orbital distribution function by way of orbital luminosity weights. We now wish to couple this information to the observed stellar population properties, in terms of the (luminosity-weighted) mean age and metallicity maps. We do this by assigning individual ages and metallicities to the fitted orbital components, noting that a ‘component’ can be composed of a collection of orbits from our dynamical model. We first consider how such components may be defined from the model, and then describe how we associate them with stellar population information to fit the observations.

4.5.1 Dynamical Selection

The Schwarzschild model provides a full description of the orbital structure of our galaxy, being a combination of many hundreds of distinct orbital families. However, in terms of coherent galaxy sub-structures, we may expect far fewer components. In our own Milky Way, for instance, we broadly recognise a ‘thin disk’, ‘thick disk’, ‘bulge’, and ‘halo’. Indeed, as described in the introduction, such broad component definitions are often applied when considering the decomposition of external galaxies. By combining both stellar kinematics and populations, we here seek to test such ‘component’ concepts.

As described above, circularity phase-space gives a simple projection of the full orbital phase-space, from which broadly-distinct orbits can be tracked as a function of radius in the galaxy. We take advantage of this as a way of potentially reducing the number of distinct orbital families by instead grouping them to form ‘components’. This brings computational advantages by reducing the dimensionality of the problem, as well as giving a clear definition of what the (dynamically, chemically, and chronologically) distinct components are. We trial a number of conceptually-different approaches here, and compare their results.

**Sample I: **Conventional Decomposition

We refer to Sample I as a ‘conventional’ decomposition because we define bins in circularity that closely approximate classic galactic components; namely, a thin disk, thick disk, and bulge. We use the ranges defined by Zhu et al. (2018b); $\lambda_{z}>0.8$, $0.25<\lambda_{z}\leq 0.8$, and $\lambda_{z}\leq 0.25$, representing the thin disk, thick disk, and bulge, respectively.

**Sample II: **Small Sampling in Circularity

In this instance, we group the orbits into small bins of $\lambda_{z}$. We do not assign physical meaning to the resulting ‘components’, but merely investigate how the extra freedom affects the stellar-population fitting.

**Sample III: **Small Sampling in Circularity and Radius

This approach samples the $2{\rm D}$ circularity space in small bins of circularity and radius. Once again, we avoid attributing physical meaning to these components, but investigate specifically if radial changes (compared to $\lambda_{z}$ sampling only) can improve the fit significantly. For consistency, the $\lambda_{z}$ bins are the same as Sample II. Also, to align with the sampling of the underlying orbits from the Schwarzschild model, the sampling of the circularity phase-space in radius is logarithmic in the inner region (but linear outside the FOV to avoid overly-large bins).

**Sample IV: **Orbit-Based Decomposition

Finally, in the limit of increasingly smaller bins in circularity, each individual orbit may be considered as a separate component. Dealing directly with orbits is conceptually a more robust implementation, as it avoids issues with the ‘projection’ of the intrinsic $3{\rm D}$ orbital phase-space into the circularity domain. It is, however, computationally expensive, and can result in an ill-conditioned problem for large numbers of orbits.

For each criterion, we consider only those components/orbits that are given non-zero weight in the best-fitting Schwarzschild model, as all others do not contribute any mass to the model.

We refer to Sample III as the illustrative case for the remainder of this work (see Section 5.1), and Fig. 7 shows the corresponding components in circularity phase-space. The gaps in the grid correspond to ‘components’ that have zero weight.

By the nature of the Schwarzschild model, we can investigate these dynamical components in greater detail, in order to gauge their spatial extent and physical properties. Fig. 8 shows the surface brightness and kinematics of a subset of characteristic dynamical components, selected to contrast the different regions in the circularity phase-space.

These properties are derived as statistical moments of the intrinsic LOSVD of each component - the integral, mean, and dispersion for the surface brightness, mean velocity and velocity dispersion, respectively. The dynamical decomposition effectively divides the LOSVD in each aperture amongst the resulting components. Each component’s LOSVD, therefore, need not be Gaussian. This is illustrated in Fig. 9, using the Sample I decomposition for clarity, and is why we restrict our visualisation of the components to the low-order kinematics.

Nevertheless, it is clear from Fig. 8 that the lowest $\lambda_{z}$ components are truly pressure-supported, while the highest $\lambda_{z}$ components are rotation-supported.

4.5.2 Dynamics/Stellar Populations Associations

The Schwarzschild model provides the luminosity-weights per aperture (Voronoi bin), for each component. The SFH provides the mean luminosity-weighted stellar-population properties per aperture. We can therefore assign such properties to the dynamical components, by equating their luminosity-weighted average with the SFH values. For a single aperture, this can be expressed as

[TABLE]

where $\Phi^{\rm SFH}$ is the luminosity-weighted stellar-population value ($t$ or $Z$) from the SFH, $\omega^{j}$ is the orbital weight of the $j$-th dynamical component, and $\phi^{j}$ is the unknown stellar-population value ($t$ or $Z$) of the $j$-th dynamical component. If we normalise the orbital weights beforehand, such that

[TABLE]

we can then express the entire problem, over all apertures, as a matrix equation

[TABLE]

where one row in the matrix (superscripts) denotes the weights for all dynamical components in a single aperture - equivalent to Eq. 6 - and one column (subscripts) corresponds to the weights of a single dynamical component in all apertures. We now simply require a linear (computationally-efficient) matrix inversion in order to find the $\vec{\phi}$, which thereby associates the dynamically-identified components with individual stellar populations. For $\phi$ and $\Phi$ in this work, we investigate both age and metallicity. Unlike other similar methods recently proposed in the literature (Long & Mao, 2018), our method does explicitly allow consideration of the $3{\rm D}$ distribution of stellar populations in projection, as the age and metallicity parametrisation explicitly adds linearly (unlike the line strength indices considered in that work). The matrix inversion is done in the framework of a Bounded-Value Least-Squares (BVLS) fit, which is akin to NNLS problems, but rather than a positivity constraint, the solution bounds are specified explicitly. Here we impose the boundary values of the SSP library in order to maintain consistency with our SFH. This is not necessary in general - for instance, if $\vec{\Phi}$ is generated without the use of an SSP library - and we emphasise that the best-fitting $\vec{\phi}$ is continuous within the bounds and not tied to the specific sampling of any SSP library. While there exists mild covariances between the stellar age and metallicity, our spectral fits from Section 3 employ a small linear regularisation. This regularisation assumes a smooth SFH in both age and metallicity in order to break the degeneracy. Moreover, in the subsequent fitting described in this section, our model needs to reproduce the spatial structure in the SFH maps, not just the age and metallicity values of an individual bin. This combination of spatial and temporal coherence allows our model to mitigate effects due to covariances. We therefore treat age and metallicity as orthogonal parameters, and simply fit them independently.

Note that the different decomposition criteria described in Section 4.5.1 are merely re-distributions of the weights $\tilde{\omega}^{j}$ into different numbers of $N_{\rm comp.}$, while Eq. 8 and the method itself remain unchanged.

5 Results

Here we present the main results of this work. We first present the various combinations of dynamical components and stellar population properties, and then how these combinations lead to inferences of the formation history of this galaxy.

5.1 Constraining the Required Number of Dynamical Components

Fig. 10 shows the fits to the measured luminosity-weighted mean stellar age and metallicity for all of the sample criteria tested in this work (see Section 4.5.1).

It is immediately clear that the conventional few-component decompositions are completely unable to reproduce the structure in the stellar populations. This implies a rather dramatic disconnect between the photometric, kinematic, and chemical structures seen in the galaxy, at least from the point of view of a handful of components. However, we note that while the surface brightness of the components from Sample I follow typical Sérsic-law expectations, the assumption that they have a single monolithic formation is strongly inconsistent with the data. This suggests that studies assuming the conventional Sérsic approach and fitting a few components should include stellar population gradients, rather than a single population per component as done here, even for evolved objects like NGC 3115.

Similarly, it appears that regardless of the degree of freedom in $\lambda_{z}$ as in Sample II, the fits are unable to reproduce the SFH, in particular the metallicity. This reaffirms the existence of radial gradients in stellar populations, with these gradients being stronger in metallicity, as has been seen by many studies previously both in general (for example, see Sarzi et al., 2018; Parikh et al., 2018; Poetrodjojo et al., 2018) and specifically for NGC 3115 (for example, see Strom et al., 1976; Pastorello et al., 2014). However, we can additionally conclude here that these gradients must exist within each component even at fixed $\lambda_{z}$.

Unsurprisingly, as the number of degrees of freedom (dynamical components) increases, so too does the quality of fit. Surprisingly, though, once the circularity phase-space is well-sampled in both dimensions as in Sample III, there already appears to be sufficient freedom in the model to reproduce the stellar-population maps, without having to consider every individual orbit. This is not necessarily expected a priori, but allows us to dramatically reduce the size of the parameter-space for subsequent analyses by using Sample III as the default dynamical decomposition moving forward. This is because, while there is a factor of 10 increase in $N_{\rm comp.}$ between Sample III and IV, their fits to the stellar-population maps are almost indistinguishable.

The fit described in Section 4.5.2 is done on each bin of the observations. However, due to the fact that orbits from the Schwarzschild model overlap in projection, multiple orbits will contribute to the mean age and metallicity of a single bin. There are necessarily different weighted mixtures of ages and metallicities that give the same average values, producing a degeneracy when associating stellar-population properties to dynamical components. One option would be to reduce the number of components until the fit is adversely affected, however this requires arbitrary choices of which components to remove. A more natural approach is to apply linear regularisation. The regularisation scheme used here prefers solutions where adjacent dynamical components/orbits have similar stellar population properties given an otherwise degenerate alternative. It is described in detail in Appendix B.

5.2 Recovering the Dynamical SFH

Following these fits and adopting the results from Sample III, we now have a set of components for NGC 3115 with both fitted kinematics and stellar populations. We have already been able to investigate the structures in components of the galaxy based on a dynamical selection criterion (see Fig. 8). Now, however, we can investigate NGC 3115 from an orthogonal perspective: the nature of spatial structures in components selected on their stellar populations. This produces spatially-resolved maps for components of the galaxy in bins of age and metallicity, which is effectively a conventional SFH, except that it originates from the spatial distribution of dynamical components. A conventional full-spectral-fitting SFH was computed in Guérou et al. (2016). For comparison, we construct the same bins in age and metallicity as that work, which are shown in Figs. 11 and 12, respectively.

By considering for now just the left columns of Figs. 11 and 12, it is clear that the bulk of the stellar mass is old, and in a spheroidal structure concentrated at the centre but extending across the FOV. We see a small portion of relatively young stars, that exist in a much more flattened structure along the plane of the galaxy. There is remarkable agreement between the purely spectral analysis in Guérou et al. (2016) and our combined approach, despite the fundamental differences between the two methodologies. Moreover, from the metallicity panels, we find a central elongated metal-rich component. This is surrounded by a near-solar diffuse component, which is slightly overdense in the centre. We find a portion of the mass in an even more diffuse metal-poor halo-like component. The absolute masses presented in each panel of Figs. 11 and 12 are computed as the fractional mass (determined by the Schwarzschild model) of MGEΣ, integrated within the MUSE FOV. This ensures that these stellar mass measurements include the information from the stellar $M_{\star}/L_{V}$ map, and are consistent with the dynamical model and the spectral SFH.

In addition to the conventional SFH analysis possible with this method, we can extract intrinsic kinematics for each of these age and metallicity bins. This is possible at present uniquely from our combined methodology. It allows us to investigate the chemistry and kinematics simultaneously and self-consistently, in a spatially-resolved manner. These kinematics are shown in the middle and right columns of Figs. 11 and 12. As mentioned in Section 4.5.1, the decomposition of the full Schwarzschild model by either dynamical (that section) or population (this section) properties necessarily divides up the LOSVD in each aperture. This is why the kinematics in Figs. 11 and 12 may appear non-physical and discontinuous. These components do not exist in isolation, and the kinematics of each component are physically meaningful only in the context of the full dynamical model. Moreover, the bins in age and metallicity are discrete, further contributing to the discontinuities in the kinematic maps. However, visualising the results in this way is useful because it allows us to investigate the relative dynamical properties of the different stellar populations in a completely spatially-resolved manner. It is in fact quite clear from these figures that, for instance, the youngest ages are heavily rotationally-supported, and that the most metal-poor components are highly pressure-supported - findings that we expand on in the following sections.

5.3 Intrinsic Stellar Velocity Dispersion and Stellar Populations

We may delve further into the model by investigating its intrinsic properties, and any direct correlations between the dynamical and stellar-population properties. To this end, we study how the properties of the intrinsic velocity ellipsoid and the assigned stellar populations are related, focusing on the vertical velocity dispersion in order to isolate features in the disk plane. We present in this section the results of this investigation, based on the intrinsic properties derived from the Schwarzschild model and the subsequent population fits.

We consider the $z$-component of the velocity dispersion, $\sigma_{z}$, in order to probe the structure and properties of the disk region of NGC 3115. As we have the complete phase space information for every orbit in our model, we can compute the intrinsic moments of the orbital motions. This includes the first ‘true’ moment of the velocity distribution along the three principle axes $(V_{x},V_{y},V_{z})$, the second moments, and all the cross-terms. The second moment of the velocity distribution can be approximated as

[TABLE]

for first and second moments $\langle V_{k}\rangle$ and $\langle V_{k}^{2}\rangle$, respectively, and velocity dispersion $\sigma_{k}$, where $k\in[x,y,z]$ is the Cartesian axis. Therefore, to compute the velocity dispersion, it follows trivially that

[TABLE]

The outputs of the Schwarzschild model are first transformed from Cartesian to cylindrical coordinates. To construct radial profiles, we transform the logarithmic energy sampling from the original Schwarzschild model from spherical coordinates to cylindrical coordinates. This produces a cylindrical radius equivalent to the original sampling in $E$. For clarity, we combine all the orbits within one PSF of the galaxy centre and treat this as a single position, as we cannot directly resolve this region. Radial profiles of the first and second moments - $\left\langle V_{z}\right\rangle(R)$ and $\left\langle V_{z}^{2}\right\rangle(R)$, respectively - are constructed by computing the luminosity-weighted sum of the orbits’ moments that lie within radial annuli. These are then converted into $\sigma_{z}(R)$ with Eq. 10. Finally, all components belonging to a particular age and metallicity bin are averaged to produce stellar-population-selected velocity-dispersion profiles. These are shown in Fig. 13.

The radial extent of each profile is determined by two factors; whether there is mass at a given radius from a given age/metallicity bin (that is, the spatial extent of the corresponding surface-brightness distribution), and whether Eq. 10 is numerically defined at a given location in a given bin, which depends on its particular combined kinematics from the dynamical model. We present the marginalised plots (over age and metallicity) in Appendix C

5.4 Modelling Systematics

Before interpreting the wealth of information in Fig. 13, we briefly discuss here the possible sources of systematic uncertainties arising from this method, and how reliable the subsequent formation history is.

In all of the dynamical samples presented in Section 4.5.1, the components are defined by hard rectangular grids in the orbital probability distribution. It is possible that this is responsible for some of the substructure that can be seen in some panels of Fig. 8. Some components appear to be a superposition of underlying dynamical structures, indicating that those components could somehow be further divided until unique dynamical features are isolated. It may be possible to mitigate these effects to some extent by introducing a more complex method for grouping dynamically-similar orbits to define the components, but such a method would have to remain contiguous in the circularity phase-space in order to conserve the mass of the dynamical model. Testing and implementing a diverse range of methods for the dynamical selection is beyond the scope of this work.

By construction, each row of the matrix in Eq. 8 - each spatial bin of the observations - is assigned the same set of $\vec{\phi}$. The consequence of this is that each component is mono-age and mono-metallicity. For the triaxial Schwarzschild model, each component is also at least point-symmetric. Therefore, the resulting fit in Fig. 10 can not reproduce any asymmetry in the SFH maps, whether that asymmetry is physical or otherwise. There is mild (non-physical) asymmetry in the SFH maps for this particular data-set. The effect of this would be a slightly larger discrepancy between the data and model maps in Fig. 10 due to the implicit averaging between the $\pm x$ sides. While we mitigate the most significant deviations by masking the right-most bins during the fit, there is still a gradual sky residual gradient across the FOV. This can be further seen in the youngest ages of the disk on the left and right sides, which would contribute to the uncertainties from that fit. However, such an issue is specific to the data-set in use here, and general applications of our method would not be subject to these uncertainties.

In Sections 5.1 and B, we address the possibility that the solutions to Eq. 8 may be unstable/degenerate. To estimate what impact this degeneracy may have specifically on our interpretation of Fig. 13, we run 100 Monte Carlo simulations by perturbing the measured SFH maps within their uncertainties, and re-fitting the dynamical components as described in Section 4.5.2. Running through the entire analysis, we generate a new set of $\sigma_{z}(R)$ curves for each trial, and these are plotted within each panel of Fig. 13. It is clear from the dense clustering of curves that the solutions are reasonably stable. It is possible for some components to traverse the bins in stellar populations, but this is only because their best-fit age and metallicity are close to the boundaries of the bins, while the absolute change in the individual ages and metallicities remains small.

More broadly, in order to accurately estimate the possible systematic uncertainties associated with a method of this nature, it is necessary to run our analysis on mock data derived from high-resolution hydrodynamical galaxy simulations. This effort is considerable in and of itself, and has been in preparation in parallel with our application here to real data (Zhu et al., in prep.). That work will test the accuracy and reliability of the recovery of stellar population properties for dynamically-selected components. However, in this work specifically, there is such clear structure in the stellar populations and kinematics due to the proximity of the galaxy and the quality of the data, that this method can be applied with confidence in its ability to recover the underlying properties of the galaxy.

6 Discussion

6.1 Galactic Components of NGC 3115

The many facets of our comprehensive model can be seen concisely in Fig. 13, in which we study the galaxy across stellar age and metallicity, and simultaneously with a metric for the spatially-resolved, intrinsic dynamical properties. We characterise these galactic features here.

Old, Metal-Rich, Compact, Hot Spheroidal Bulge

Firstly, we find a bulge-like component at the oldest age, with a mild spread in metallicity. This is clear from the central, compact spheroidal peak in surface brightness in the three most metal-rich panels at $14.2\$\mathrm{G}\mathrm{y}\mathrm{r}; panels $(f)$, $(l)$, and $(r)$. This is accompanied by a corresponding central $(R\lesssim 15$\mathrm{\SIUnitSymbolArcsecond}$)$ peak in $\sigma_{z}$ in these bins $(\sigma_{z}\gtrsim 100\ $\mathrm{k}\mathrm{m}\ \mathrm{s}^{-1}$)$, indicating that this region is pressure-supported. This component has only undergone mild chemical enrichment, likely due to secular stellar-evolution processes, with no other significant evolution - neither in its kinematics nor populations. By making the appropriate cuts in age, metallicity, circularity, and radius in order to isolate only the bulge contributions to the relevant panels, we estimate that this component contains $\sim 3.9\times 10^{9}\ $\operatorname{$\mathrm{M}_{\odot}$}.

Metal-Rich, Extended, Cold Disk

We also find evidence for a disk-like structure that is enriched, and present at all ages. This is clear from the surface brightness in panels $(a)-(f)$, and also panel $(h)$, which all show elongated structures in the plane of the galaxy. Again, this is corroborated by the dynamics, which show that this component has $\sigma_{z}\lesssim 50\$\mathrm{k}\mathrm{m}\ \mathrm{s}^{-1} at intermediate radii where the disk dominates, implying that it is supported primarily by coherent rotation (cold orbits). In particular, even panel $(f)$ has a signature of the disk in both the surface brightness and the velocity dispersion. This is an interesting finding, as a galactic component that is simultaneously cold, old, and enriched would be sensitive to any dynamical perturbations over the galaxy’s history. Therefore, the mere existence of such an old population in a disk configuration immediately implies a fairly quiescent history for NGC 3115. We see that star formation in the disk continued smoothly from the earliest times, gradually declining until it ceased $\sim 4\ $\mathrm{G}\mathrm{y}\mathrm{r} ago. The resulting stars are naturally enriched, and are progressively more dynamically cold - a progression that we quantify further in Section 6.4. In our model, the disk component is $\sim 1.6\times 10^{10}\$\operatorname{$\mathrm{M}_{\odot}$}$$.

Old, Metal-Poor, Diffuse, Hot Stellar Halo/Thick Disk

Finally, we find a diffuse halo/thick-disk-like component that is dynamically hot ($\sigma_{z}\gtrsim 100\$\mathrm{k}\mathrm{m}\ \mathrm{s}^{-1}, including at large radius), metal-poor $(\nuclide{[Z/H]}\lesssim-0.50\ $\mathrm{d}\mathrm{e}\mathrm{x}$)$, and mostly old $(t\gtrsim 10\ $\mathrm{G}\mathrm{y}\mathrm{r}$)$. This can be seen in panels $(p)-(x)$, which all have a similar featureless spatial distribution, and relatively high $\sigma_{z}$. Panels $(j)$ and $(l)$ appear to have contributions from this component as well, but superimposed with contributions from the disk and bulge, respectively. This stellar halo/thick-disk component contains $\sim 4.2\times 10^{10}\ $\operatorname{$\mathrm{M}_{\odot}$}.

6.2 Implications for the Formation History of NGC 3115

The current paradigm for galaxy formation, in particular of ETG, details the hierarchical assembly of massive galaxies via the merging of smaller systems, which result in expectedly pressure-supported systems due to the nature of the merging process (for instance, see Somerville & Davé, 2015, and references therein for a general review). Merging during the lifetime of a galaxy is a key aspect, and often the focus, of many cosmological and/or hydrodynamical simulations (for example, see Bird et al., 2013; Athanassoula et al., 2016; Eliche-Moral et al., 2018). Moreover, merging is often invoked to explain many phenomena seen in observational studies (for example, see Arnold et al., 2011; Lidman et al., 2013; Guérou et al., 2016). However, these mergers are typically assumed to have been of two (or more) fully-formed progenitor galaxies, implying that sufficient cosmological time had passed for such galaxies to first form, then merge. As discussed above, we find a very old and dynamically cold population in the disk of NGC 3115, implying that any destructive major merger had to have occurred prior to the formation of this disk. Laurikainen et al. (2013) arrived at the same conclusion when they detected cold ‘lenses’ (an elliptical structure with a sharp inner edge in surface brightness; Kormendy, 1979) within ${\rm S}0$ galaxies that were very old - also $\sim 14\$\mathrm{G}\mathrm{y}\mathrm{r}$$. The Universe had, at such early times, considerably different physical conditions to the environments in which the simulated progenitor galaxies were formed. For instance, the progenitor galaxies in Athanassoula et al. (2016) were explicitly modelled after “nearby galaxies”. As a result, such scenarios may not appropriate for our inferences here, and we therefore look for a possible alternate formation mechanism to explain the results we see in Fig. 13.

One such mechanism is ‘compaction’ (Dekel et al., 2009; Zolotov et al., 2015) at early times. Compaction involves the rapid dissipative contraction of highly-perturbed gas into the central regions of galaxies in the form of cold streams. These streams can trigger star-formation in this region, which is reportedly quenched rapidly for massive galaxies (Zolotov et al., 2015). Therefore, compaction could also explain the enrichment gradient in the central bulge component, which is all at early times. Moreover, since the gas cooling is dissipative, and the cold streams are not as disruptive as major mergers, this method of formation could also explain the persistence of the old cold disk structure (Dekel et al., 2009), even if the streams entered after the old disk formed. The remaining gas which had been accreted via the cold streams would go on to form the main stellar disk at progressively younger ages until the gas reservoir is exhausted.

Irrespective of the formation mechanism for the central bulge and main disk, the outer stellar halo/thick-disk component is strongly consistent with an accreted origin, mainly from dry minor mergers. This is due to the combination of old ages, low metallicities, pressure-supported kinematics, and featureless yet extended surface brightness. This component forms a significant portion $(\sim 0.68M_{\star})$ of the stellar mass of NGC 3115. Interestingly, the cosmological simulations of Oser et al. (2010) predict that for ‘intermediate mass’ galaxies such as NGC 3115, the accreted material constitutes $\sim 65\%$ of the stellar mass at $z=0$, which is remarkably consistent with our model. Given such a significant contribution to the stellar mass budget, persistent minor mergers could also explain how NGC 3115 transformed into an ${\rm S}0$ galaxy, by building up the ‘thick’ disk and diluting any spiral arms that may have been present in the progenitor object. While previous works have claimed that either environmental perturbations (Bekki & Couch, 2011), internal disk instabilities (Saha & Cortesi, 2018), or major mergers (Querejeta et al., 2015; Tapia et al., 2017; Diaz et al., 2018; Fraser-McKelvie et al., 2018) are the likely formation paths for ${\rm S}0$ morphologies, our model is inconsistent with these mechanisms. NGC 3115 is a field ${\rm S}0$, making environmental perturbations unlikely, and both internal disk instabilities and major mergers would likely destroy the old disk structure that we find in our model.

6.3 Comparison to Previous Inferences

Due to its proximity, NGC 3115 has been studied widely in the literature, using many techniques that produce orthogonal constraints on the inferred formation history. We investigate here how our results compare with these works.

Arnold et al. (2011) use Supreme-Cam imaging and long-slit GC observations to study the kinematic and metallicity properties of NGC 3115. They found that an early violent major merger would explain the relatively high rotation and flattening of the central bulge component. Guérou et al. (2016) proposed that a small number of progenitors with \log_{10}(M_{\star}/\operatorname{\mathrm{M}_{\odot}})\sim 10 is consistent with both the surviving angular momentum of this central component, as well as its enrichment at such early times.

Guérou et al. (2016) suggested that some portion of the gas reservoir that was present during the major merger survived the interaction, went on to cool and eventually form the dynamically cold, younger disk stars. Our results indicate that the cooling of the gas and the formation of new stars happens ‘immediately’ (within the oldest age bin), and stars continued to form on progressively colder orbits through to the $4\$\mathrm{G}\mathrm{y}\mathrm{r}$$ bin, until this gas was consumed. Their interpretation of the excess gas does indeed agree with our results, however our model does not need to invoke a major merger as the source, and an alternate origin for the gas could be the cold streams associated with an early compaction phase.

Arnold et al. (2011) claim that steeper metallicity gradients can also be the result of passive accretion of low-mass, low-metallicity satellites which lower the average metallicity in the outskirts, thereby steepening the gradient. This accretion is also consistent with the observed metallicity map for NGC 3115. Brodie et al. (2012) and Cantiello et al. (2014) find strong evidence that the observed colour bimodality in the GC population in NGC 3115 is driven by an underlying metallicity bimodality. This supports not just an accreted origin for the blue population of GC, but specifically accretion from dry minor-mergers, locking in the low metallicity of the in-falling objects. Our metal-poor, dynamically-hot, diffuse components also strongly favour passive accretion. This is because low-mass objects would be relatively metal-poor, and many such in-falling objects would impact at arbitrary angles imparting no net angular momentum, while increasing the $\sigma_{z}$.

6.4 Resolved Studies and Galaxy Formation Simulations

By exploiting the intrinsic, spatially-resolved properties of our model, we can begin to compare directly to results from resolved studies of the Milky-Way and Local Group, as well as cosmological hydrodynamical galaxy formation simulations. This allows us to leverage the look-back capabilities of the simulations as well as the resolving power of the local observations in order to strengthen the interpretation of our results. Moreover, by extending the intrinsic correlations seen in the Milky Way to external galaxies, we begin to gauge how unique the Milky Way is.

Age-velocity dispersion relations have been found in many resolved (Nordström et al., 2004; Grieves et al., 2018) and simulated (Martig et al., 2014) studies. They primarily conclude that stars born earlier have a higher velocity dispersion compared to those formed later. The underlying physical cause of this relation, however, has evaded a general consensus in the literature. Some scenarios broadly claim that the increase in velocity dispersion is a dynamical effect of internal interactions that build up over time, with a number of different specific mechanisms being proposed as the culprit. Since the older stars have been experiencing these interactions for a longer period of time, they should therefore show the largest increase in velocity dispersion. Evidence in favour of such internal mechanisms has come from both observations (Yu & Liu, 2018) and simulations (Saha et al., 2010; Grand et al., 2016; Aumer et al., 2016). An alternative explanation is that in the early Universe, conditions were generally more chaotic (Wisnioski et al., 2015), so that any stars born at that time were more likely to have higher velocity dispersion. As conditions gradually settled over cosmic time, stars were being born in progressively lower-dispersion conditions. A number of studies have identified the conditions at birth as the dominant effect in determining a population’s present-day velocity dispersion, again both from resolved observations (Leaman et al., 2017) and simulations (Bird et al., 2013; Ma et al., 2017). Finally, the comparison between observations and simulations in Pinna et al. (2018) has identified the underlying complexity and inherent degeneracy in discriminating between these scenarios. They claim that many of the effects described above likely play a role to a varying degree, and that the imprint of some mechanisms fade over time, further complicating any attempt to constrain the physical cause of disk heating.

Our model allows us to estimate the intrinsic age-velocity dispersion relation, even though NGC 3115 is unresolved. It can be seen from Fig. 13 that, at fixed metallicity (the enriched components in particular), there is a mild increase in velocity dispersion with age, notably at all radii. This is quantified in Fig. 14, for the disk region of our model. To isolate this region, we make a conservative cut in circularity of $\lambda_{z}>0.5$ in order to exclude the most pressure-supported orbits, as well as select the highest metallicity bins.

Owing to how distinct our methodology is for the derivation of these data, comparison to other works is complicated. Specifically, observations of the stellar velocity dispersion across redshift are difficult to obtain. We therefore look for qualitative comparisons to other works. For instance, Wisnioski et al. (2015) measured, and compiled literature observations of, the gas velocity dispersion as a function of redshift. These data are taken from a range of surveys containing galaxies with \log_{10}(M_{\star}/\operatorname{\mathrm{M}_{\odot}})\in[10.1,11.0]. We include these data in Fig. 14, which are specifically from HERACLES (Leroy et al., 2009), DYNAMO (Green et al., 2014), GHASP (Epinat et al., 2010), PHIBBS (Tacconi et al., 2013), MASSIV (Epinat et al., 2012), OSIRIS (Law et al., 2009), AMAZE-LSD (Gnerucci et al., 2011), SINS (Schreiber et al., 2009) and zC-SINF (Schreiber et al., 2014), and KMOS3D (Wisnioski et al., 2015). Pillepich et al. (2019) investigate both the gaseous and stellar velocity dispersion as a function of redshift in the IllustrisTNG cosmological simulations for star-forming galaxies. They find the same general redshift evolution as we see here for both gas and stars in their simulations. We emphasise here, however, that the data types, galaxy types, and techniques for measuring the velocity dispersion are not directly comparable between these three works. Both Wisnioski et al. (2015) and Pillepich et al. (2019) study ensemble population properties of gas and/or stellar kinematics in star-forming galaxies, while we consider only the stellar kinematics of a single, relatively-quiescent galaxy. Nevertheless, that the shape of this redshift evolution is so consistent between these works suggests that it is robust.

Interestingly, we also see in Fig. 13 an increase in velocity dispersion towards more metal-poor components, even at fixed age and again at all radii. This implies that metallicity has a significant contribution to the increase in $\sigma_{z}$, which is obscured when considering age at fixed metallicity as in Fig. 14, or when marginalising over all metallicities as in Fig. 16.

We can thus further constrain the physical mechanism driving the age-velocity dispersion relation by leveraging the formation scenario we’ve established in the previous sections. Following the reasoning in Section 6.2, the presence of a relatively cold, yet very old disk implies an upper limit on how much internal ‘heating’ could have occurred in this galaxy, since such heating would gradually inflate the scale height of the disk’s old stellar population over the galaxy’s lifetime. Furthermore, any internal dynamical interactions that drive an increase in $\sigma_{z}$ over time should be completely impartial to the metallicity of those stars. This is inconsistent with the rather significant change in $\sigma_{z}$ with metallicity at fixed age in our model. Finally, if we interpret at face value the absolute agreement between our stellar velocity dispersion evolution and that of gas measurements in Fig. 14, this suggests that the stars inherited their dynamical properties from the gas and there was little-to-no subsequent evolution. Once again, this is suggestive of minimal, if any, amounts of disk heating throughout the entire evolutionary history of NGC 3115. Therefore, while we can not make strong conclusions about the true physical interpretation of the age-velocity dispersion relation, nor definitively exclude any disk heating over the lifetime of the galaxy, our model indicates that the conditions at birth have the most dominant impact on a population of star’s present day velocity dispersion.

6.5 Non-Axisymmetric Structures

The possibility of a bar-like structure and associated resonances within NGC 3115 has been discussed in the literature (Kormendy & Richstone, 1992, 1995; Guérou et al., 2016). These works all conclude that there is tentative evidence for such structures, but they can neither confirm nor refute their definitive existence. Guérou et al. (2016) finds many coincident kinematic features that would indicate the presence of a bar. They show, however, that a simple $N$-body model without a bar is able to reproduce these features. We can investigate our dynamical model for evidence of kinematic signatures indicative of a bar. We find no such evidence in the orbital composition, which is dominated by short-axis tube orbits everywhere, with only a small gradual increase in the contribution from box orbits towards the centre. Our model reproduces the positive correlation between $V$ and $h3$ in the central off-axis region that is usually associated with a bar, despite being built in a static gravitational potential. We conclude therefore, as Guérou et al. (2016), that a bar is not strictly necessary to form this feature, but whose existence can not be conclusively ruled out by our model.

In general, the box/long-axis-tube orbits are a unique feature of a triaxial model, and the presence of a small fraction of these orbits towards the centre of our best-fit model implies a degree of triaxiality - and may also imply a bar. We measure a very small constant oblate-triaxiality over the full model of $T=(1-p^{2})/(1-q^{2})=0.042$, with $p=b/a$ and $q=c/a$, for the major, intermediate, and minor axis lengths, $a$, $b$, and $c$, respectively. Despite the slight increase in non-axisymmetric orbits towards the centre, there is no corresponding increase in triaxiality. However, we note here that the potential for strong triaxiality is limited in our model, since we have fixed the intrinsic shape $u=0.999$, and are using a projected mass model which has only a single position angle. More accurate constraints on the triaxiality would thus require relaxing both of these assumptions.

We therefore find that while our static triaxial model is able to fully reproduce the kinematics, conclusive evidence of a bar (or otherwise) would require a more direct modelling approach including a tumbling time-variable gravitational potential, such as the NMAGIC made-to-measure models (de Lorenzi et al., 2007; Morganti & Gerhard, 2012). More importantly for the general method presented here, Zhu et al. (2018a) showed that even for an intrinsically barred galaxy, this triaxial Schwarzschild code is able to accurately recover the underlying orbital distribution for the non-bar/resonance regions.

7 Conclusions

In this work, we have presented the application of a conceptually-new approach for the deduction of formation histories from IFU observations of external unresolved galaxies. We extracted and presented stellar kinematics and stellar-population properties across a $\sim 240$\mathrm{\SIUnitSymbolArcsecond}$$ FOV of a nearby ${\rm S}0$ galaxy, NGC 3115. By fitting detailed Schwarzschild orbit-superposition dynamical models to the kinematics, we defined components within the galaxy that are dynamically distinct. These components were then assigned a mean stellar age and metallicity in order to reproduce the observed stellar population properties. This combination of spatial resolution, kinematics, and stellar populations allows us establish a complete history of the formation of NGC 3115:

• •

We find that in the early gas-rich Universe, cold streams funnelled into the core of the progenitor of NGC 3115 early in its formation. These streams caused the compaction of the bulge and its mild metallicity gradient

• •

The remaining gas from this event cooled and formed stars, which began shortly after the compaction phase. Star formation continued (though declining) through to the youngest stars on the coldest orbits and with the highest metallicity

• •

Meanwhile, many low-mass satellites were being accreted, fleshing out the halo/thick-disk with low-metallicity material, gradually converting NGC 3115 into its present-day early-type lenticular morphology

More generally, we have combined the stellar dynamics and populations in a comprehensive and self-consistent manner. This has allowed us to empirically conclude that conventional galactic decomposition techniques - with few components - are unable to simultaneously fit a galaxy’s shape, kinematics, and stellar populations, unless gradients are considered. We have also determined that in the case of NGC 3115, the conditions in which a population of stars forms has the dominant effect on their observed present-day kinematics. This approach amounts to a significant step towards a completely simultaneous population-dynamical model, that will drive progress in the field of galaxy formation with remarkable detail and accuracy. In future work, we will incorporate other stellar-population parameters (namely, spatially-resolved measurements of the IMF and abundances) to further improve the accuracy of the combined model presented here.

This method is currently being tested on mock data in order to estimate the reliability of such an approach. Leveraging the simulations will also inform this methodology on which physical properties are the most important for discriminating between the different formation paths that built up the galaxy. Finally, the successful application of this methodology to a sample of galaxies in different environments will be able to uncover the dominating mechanism(s) of formation during these galaxies’ histories.

Acknowledgements

We thank Adrien Guerou and Eric Emsellem for providing the reduced and calibrated MUSE data cube and the photometric MGE model, as well as for the stimulating discussion. We thank Emily Wisnioski for providing the literature measurements for Fig. 14. RMcD is the recipient of an Australian Research Council Future Fellowship (project number FT150100333). LZ acknowledges support from Shanghai Astronomical Observatory, Chinese Academy of Sciences under grant NO.Y895201009. GvdV acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 724857 (Consolidator Grant ArcheoDyn). We acknowledge support from the Australia-Germany Joint Research Cooperation Scheme funded by Universities Australia and the German Academic Exchange Service (DAAD). This work makes use of the queenebee compute cluster at Max-Planck-Institut für Astronomie, and the OzStar supercomputer at Swinbourne University. The software used for data analysis include matplotlib (Hunter, 2007), AstroPy (Astropy Collaboration et al., 2013), and the SciPy collection (Jones et al., 2001). We finally thank the anonymous referee, whose feedback greatly improved the clarity of this work.

Appendix A Mass Surface Density MGE

We present in Table 2 the MGE fit to the mass surface-density ‘image’ described in Section 4.1.1, which requires $16$ Gaussians in order to accurately describe the mass distribution.

Appendix B Regularisation in the Dynamical Star-Formation History

By fitting the spatially-resolved luminosity-weighted maps of stellar-population properties (described in Section 4.5.2) with a given number of dynamical components, there is inevitably some level of degeneracy between solutions to Eq. 8, which increases with increasing $N_{\rm comp.}$. In order to reduce the impact of this degeneracy on the analyses that followed, we implemented a linear regularisation into the BVLS fit. This regularisation, as is widely used in astrophysical problems, penalises solutions which vary sharply in the parameter-space under consideration. For instance, regularisation in the context of Schwarzschild dynamical models would penalise solutions which have significantly different contributions from neighbouring orbits that have similar physical properties. The motivation for this is that most physical systems should undergo changes in their properties smoothly across physical parameters, rather than discretely.

In our application of regularisation to Eq. 8, we wish to impose a smoothness in the distribution of ages and metallicities that are assigned to the dynamically-selected components. In this way, we prefer solutions (which would otherwise be degenerate) that assign similar stellar-population properties to dynamical components that have similar physical properties.

For the regularisation, we minimise a linear approximation (as it is a linear problem) to the integral of the second-order Laplacian. Given in §19.5 of Press et al. (2007), this can be expressed as the following:

[TABLE]

This formalism assumes that sequential $\mu$ are adjacent in the physical parameters of interest - in the case of our dynamically-selected components, this means neighbouring in $\lambda_{z}$ and $R$. However, as seen in Fig. 7, this is not the case for some components; where sequential labels wrap around to the next column, they have significantly different physical properties. In order to retain the information about which components are truly adjacent in physical properties, and avoid regularising over non-neighbouring components, we instead construct a $2{\rm D}$ and $3{\rm D}$ regularisation matrix for the dynamical components and orbits, respectively (which are defined by 2 and 3 parameters, respectively). This matrix preserves the second-order Laplacian given in Eq. 11 for each dimension, and is unique for every set of $[\mu,\mu+2]$. In order to apply the matrix to the BVLS fit, it is flattened into a single row in the same way in which the components and orbits are reduced to a single dimension along the $N_{\rm comp.}$ axis of Eq. 8. In this way, the location of the regularisation constraints preserves the memory of which components are neighbouring in physical parameters. Since Eq. 11 affects only 3 components, this operation is repeated for each set of $[\mu,\mu+2]$, adding a row to Eq. 8 each time. In practise it is implemented in an analogous fashion to what is described in Cappellari (2016). The effect of this is illustrated in Fig. 15, which compares a completely free, unregularised fit with a regularised fit which was used for the results presented in this work.

This figure illustrates that the regularisation does indeed act in the desired way, by producing ‘smoother’ variations between the grid cells, compared to the unregularised fit. Importantly, there is a statistically-insignificant difference between the $\chi^{2}$ values of the two fits, implying the regularisation is indeed acting only to break the degeneracy.

Appendix C Marginalised Intrinsic Velocity Dispersion Profiles

In Figs. 16 and 17, we present the radial profiles of the intrinsic vertical velocity dispersion for the dynamically-selected components of our model, binned only in age, and only in metallicity, respectively.