Ongoing deep integralfield spectroscopy IFS observations of disc galaxies provide opportunities for comparison with the Milky Way MW to understand galaxy evolution. However such comparisons are marred by many challenges such as selection effects differences in observations and methodology and proper validation of fullspectrum fitting methods. In this study we present a novel code GALCRAFT to address these challenges by generating mock IFS data cubes of the MW using simple stellar population models and a mock MW stellar catalogue derived from EGALAXIA. We use the widely adopted fullspectrum fitting code PPXF to investigate the ability to recover kinematics and stellar populations for an edgeon mock MW IFS observation. We confirm that differences in kinematics mean age inlineformulatexmath idTM0001 notationLaTeXmathrmMHtexmathinlineformula and inlineformulatexmath idTM0002 notationLaTeXmathrmalpha Fetexmathinlineformula between thin and thick discs can be distinguished. However the age distribution is overestimated in the ranges between inlineformulatexmath idTM0003 notationLaTeX24texmathinlineformula and inlineformulatexmath idTM0004 notationLaTeX1214texmathinlineformula Gyr compared to the expected values. This is likely due to the age spacing and degeneracy of SSP templates. We find systematic offsets in the recovered kinematics due to insufficient spectral resolution and the variation of lineofsight velocity distribution with age and inlineformulatexmath idTM0005 notationLaTeXmathrmMHtexmathinlineformula. With future higher resolution and multiinlineformulatexmath idTM0006 notationLaTeXmathrmalpha Fetexmathinlineformula simple stellar population templates GALCRAFT will be useful to validate key signatures such as inlineformulatexmath idTM0007 notationLaTeXmathrmalpha Fetexmathinlineformulainlineformulatexmath idTM0008 notationLaTeXmathrmMHtexmathinlineformula distribution at different R and inlineformulatexmath idTM0010 notationLaTeXztexmathinlineformula and potentially infer radial migration and kinematic heating efficiency to study detailed chemodynamical evolution of MWlike galaxies.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Validating fullspectrum fitting with a synthetic integralfield spectroscopic observation of the Milky Way
{'Validating full-spectrum fitting with a synthetic integral-field spectroscopic observation of the Milky Way': 'Zixian Wang ( 王 梓 先 ), 1 , 2 , 3 ★ Sanjib Sharma, 1 , 2 , 4 † Michael R. Hayden, 1 , 2 , 5 , 6 † Jesse van de Sande, 1 , 2 , 6  \nJoss Bland-Hawthorn, 1 , 2 , 7 Sam Vaughan, 1 , 2 , 8 , 9 Marie Martig, 10 Francesca Pinna 11 , 12 , 13  \n- 1 Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, NSW, 2006, Australia\n- 2 ARC Centre of Excellence for All Sky Astrophysics in Three Dimensions (ASTRO-3D)\n- 3 Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA\n- 4 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 5 Homer L. Dodge Department of Physics & Astronomy, University of Oklahoma, 440 W. Brooks St., Norman, OK 73019, USA\n- 6 School of Physics, University of New South Wales, NSW, 2052, Australia\n- 7 Miller Professor, Miller Institute, UC Berkeley, Berkeley CA 94720\n- 8 Astronomy, Astrophysics and Astrophotonics Research Centre, Macquarie University, Sydney, NSW 2109, Australia\n- 9 Centre for Astrophysics and Supercomputing, School of Science, Swinburne University of Technology, Hawthorn, VIC 3122, Australia\n- 10 Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK\n- 11 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany\n- 12 Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, E-38205 La Laguna, Tenerife, Spain\n- 13 Departamento de Astrofísica, Universidad de La Laguna, Av. del Astrofísico Francisco Sánchez s/n, E-38206, La Laguna, Tenerife, Spain  \nAccepted 2024 September 12. Received 2024 August 26; in original form 2023 October 27', 'ABSTRACT': 'Ongoing deep IFS observations of disk galaxies provide opportunities for comparison with the Milky Way (MW) to understand galaxy evolution. However, such comparisons are marred by many challenges such as selection effects, differences in observations and methodology, and proper validation of full-spectrum fitting methods. In this study, we present a novel code G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc to address these challenges by generating mock IFS data cubes of the MW using simple stellar population models and a mock MW stellar catalog derived from E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc. We use the widely adopted full-spectrum fitting code /p.pcPXF to investigate the ability to recover kinematics and stellar populations for an edge-on mock MW IFS observation. We confirm that differences in kinematics, mean age, [ M / H ] , and [ 𝛼 / Fe ] between thin and thick disks can be distinguished. However, the age distribution is overestimated in the ranges between 2 -4 and 12 -14 Gyr compared to the expected values. This is likely due to the age spacing and degeneracy of SSP templates. We find systematic offsets in the recovered kinematics due to insufficient spectral resolution and the variation of line-of-sight velocity distribution with age and [ M / H ] . With future higher resolution and multi-[ 𝛼 / Fe ] SSP templates, G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc will be useful to validate key signatures such as [ 𝛼 / Fe ] -[ M / H ] distribution at different 𝑅 and | 𝑧 | and potentially infer radial migration and kinematic heating efficiency to study detailed chemodynamical evolution of MW-like galaxies.  \nKey words: Galaxy: stellar content - galaxies: stellar content - galaxies: kinematics and dynamics - galaxies: star formation methods: numerical - techniques: spectroscopic', '1 INTRODUCTION': "Howgalaxies form and evolve with time remains one of the outstanding questions facing astrophysics today. Over the last decade, integralfield spectroscopy (IFS) instruments have enabled us to obtain integrated light observations of galaxies across different regions of the same object, which have much better spatial coverage than singleor multi-slit spectroscopy. Several IFS surveys have already been carried out such as SAMI (Croom et al. 2012), CALIFA (Sánchez et al. 2012), and MaNGA (Bundy et al. 2015) to measure properties of thousands of galaxies to understand their evolution history.  \n- ★ E-mail: wang.zixian.astro@gmail.com\n- † These authors contributed equally to this work.  \nIn extragalactic analysis, an integrated spectrum is assumed to be a weighted superposition of many simple stellar population (SSP) spectra and the spectral features are reshaped by line-of-sight velocity distributions, dust attenuation and gas emission (Conroy 2013). To decompose these complex components and measure the kinematics and stellar population properties, many stellar population synthesis software have been developed such as /p.pcPXF (Cappellari & Emsellem 2004; Cappellari 2017, 2023), STARLIGHT (Cid Fernandes et al. 2005), STECMAP(Ocvirket al. 2006), P/i.pc/p.pc/e.pc3D (Sánchez et al. 2016), and F/i.pc/r.pc/e.pc/f.pc/l.pc/y.pc (Wilkinson et al. 2017) to perform full-spectral fitting to the observational spectra. They have been broadly used to provide spatially resolved galaxy properties in data releases for a wide variety of scientific topics.  \nDue to our unique vantage point, the Milky Way (MW) is by far  \nthe best-studied galaxy in the Universe. It provides an ideal test bed for our understanding of galaxy formation and evolution, and is often used to compare with extragalactic studies. Precise astrometric data from Gaia (Gaia Collaboration et al. 2023) and accurate chemical abundances of individual stars from large spectroscopic surveys such as LAMOST (Zhao et al. 2012), GALAH (Buder et al. 2021) and APOGEE (Majewski et al. 2017) have been analysed in the last decade for Galactic archaeologists to reveal the detailed chemodynamical picture of the Milky Way (Freeman & Bland-Hawthorn 2002; Bland-Hawthorn & Gerhard 2016). This includes the [ 𝛼 / Fe ] -[ Fe / H ] bimodal distribution (e.g., Nidever et al. 2014; Hayden et al. 2015; Sharma et al. 2021b), radial migration (e.g., Schönrich & Binney 2009), accretion history and the interplay of chemical and dynamical processes (e.g., Xiang & Rix 2022). In the past decades, many Galactic chemical evolution (GCE) models have been developed to reproduce the [ 𝛼 / Fe ] bimodality and uncover its origin (e.g., Haywood et al. 2018; Spitoni et al. 2019; Sanders & Binney 2015; Sharma et al. 2021b; Chen et al. 2023). However, because all of these models can reproduce consistent chemodynamical distributions with observations, the origin of [ 𝛼 / Fe ] -bimodality is still under debate. Moreover, the MW is only one galaxy, whether the formation theories apply to other disk galaxies or whether the MW is unique in the Universe is still an open question. Therefore, it is essential to combine the MW with other Milky Way analogous (MWAs) to address all these questions (van de Sande & Scott 2021).  \nSeveral studies compared face-on MWAs from extragalactic IFS surveys with the MW: Boardman et al. (2020) selected 62 MWAs from MaNGA with the criteria of stellar masses and bulge-to-total ratios. They found most of these galaxies have flatter stellar and gas-phase metallicity gradients due to different disc scale lengths. A greater consistency can be found when scaling gradients by these lengths; Zhou et al. (2023) revisited 138 MaNGA galaxies by fitting the spectra with a semi-analytic chemical evolution model (Zhou et al. 2022) and measured their star formation and chemical enrichment histories. They detected similar [ 𝛼 / Fe ] bimodality as the Galactic APOGEEobservations (Abdurro'uf et al. 2022) in many of the galaxies. Compared to face-on galaxies, edge-on MWAs are more useful in providing elemental abundance distributions at different 𝑅 and | 𝑧 | . Several studies of nearby edge-on MWAs and lenticular galaxies using MUSE found similar kinematics and [ 𝛼 / Fe ] differences between the thin and thick disks (e.g., Pinna et al. 2019a,b; Scott et al. 2021; Martig et al. 2021). In particular, Scott et al. (2021) demonstrated that UGC 10738 has similar [ 𝛼 / Fe ] -[ Fe / H ] distributions at different projected 𝑅 and | 𝑧 | with the MW observation in Hayden et al. (2015) and model predictions in Sharma et al. (2021b), which supports the commonness of MW's chemical distributions.  \nDespite the efforts above, a direct comparison of the MW with its analogs in kinematics and chemistry is still challenging because the observables and methods used for studying the MW are different from those for external galaxies, i.e., utilizing properties of individual stars as opposed to integrated quantities from stellar populations with projection effects. Therefore, the comparisons may be impacted by systematics or biases (Fraser-McKelvie et al. 2019; Boardman et al. 2020). In addition, some key processes such as radial migration and kinematic heating have not been extensively explored like the MW on external galaxies, which are also essential to identify whether the MW's formation and evolution is distinct from the general population. Therefore, to take the MW as an ideal laboratory and test its evolution theories on other MWAs, one needs to remove these observational and methodological biases, transfer the observables of MWandexternal systems into the same definition, and apply models that consider both chemical enrichment and kinematic processes.  \nIn addition, even though the full-spectrum fitting methods have been useful in providing kinematics and stellar population properties of different components of external galaxies, due to the high degeneracy of all the parameters, it is still difficult to verify how close the measured values are to the 'actual' galaxy properties and whether any potential bias is introduced by the measuring methods. This could cause different fitting methods or strategies to obtain different results and lead to different conclusions. Therefore, every individual study has to validate the reliability of their measurements (e.g., Kacharov et al. 2018) and the conclusions are often only made qualitatively. Several studies examined the robustness of different full-spectrum fitting software and investigated effects due to different spectral models (e.g., Ge et al. 2018, 2019; Woo et al. 2024). However, they only tested mock spectra generated using one or two stellar populations, or only tested on global stellar population properties without spatial distribution. Until now, there has not been a comprehensive exploration to perform full-spectrum fitting on a realistic mock galaxy with chemical enrichment and kinematic processes and investigate the recovery ability of spatially resolved properties in detail. All the above considerations lead to the development of tools and studies presented in this work.  \nIn this study, to address the challenges of comparing the MW to other galaxies, we present a novel code G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc to generate mock integral-field spectrograph (IFS) data cubes of the MW with integrated spectra using simple stellar population (SSP) models and mock stellar catalog from E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc (Sharma et al. in perp.), which is based on the chemodynamical model of Sharma et al. (2021b) (hereafter S21) that was verified to be consistent to the current Galactic observations in both kinematics and chemistry across a range of 𝑅 and | 𝑧 | . This mock data cube is in the same format as extragalactic IFS observations. Therefore, extragalactic data analysis methods (e.g., the GIST pipeline Bittner et al. 2019 with /p.pcPXF Cappellari & Emsellem 2004; Cappellari 2017, 2023) can be applied to the mock data cube to measure directly comparable parameters in (age, [ M / H ] , [ 𝛼 / Fe ] , 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 , light/mass fraction distributions). To address the reliability of full-spectrum fitting methods in measuring galaxy properties, we apply /p.pcPXF on G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc-generated mock edge-on IFS cubes and compare the measured kinematics and stellar population properties to the input values according to the GCE model S21. This paper provides a comprehensive test of full-spectrum fitting methods on a mock galaxy that considers physical chemodynamical processes at different galaxy components. The results will be useful as a reference for future extragalactic data analysis using similar spectrum fitting algorithms.  \nWe describe the ingredients used in G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc and detailed procedures of this code in Section 2. In Section 3, we test the ability of broadly used spectral fitting tool /p.pcPXF to recover kinematics, stellar population parameters, and light/mass fraction distributions. In Section 4, we discuss the causes of deviations between measured and true (input) values and address potential reasons and future improvements on current algorithms for obtaining better results. We explore the effect of different fitting strategies and provide some suggestions when using such spectral fitting tools for future studies. We also give some caveats about using mock cubes generated from G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code. In Section 5, we talk about future plans on the improvements of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc along with spectral fitting methods and prospect some scientific topics that can be done by using G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc. A summary is presented in Section 6.", '2 DATA CUBE GENERATION': 'The purpose of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc is to take the mock stellar catalog of the MW obtained from E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc (Sharma et al. in prep.), and transform it into a mock data cube in three dimensions (two in spatial and one in spectral) as observed by integral-field spectrographs (IFS) such as PMAS/PPak (Roth et al. 2005; Kelz et al. 2006), AAT/SAMI(Croomet al. 2012), Sloan/MaNGA (Bundy et al. 2015), and VLT/MUSE (Bacon et al. 2010). The G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code takes as input the following set of user-specified input parameters: galaxy distance ( 𝑑 ), inclination ( 𝑖 ), extinction, SSP templates, instrumental properties, as well as a few additional parameters (see the full list in Table 1). The produced mock data cube can be processed in the same way as real IFS observations data by many methods, particularly codes like Voronoi binning (Cappellari & Copin 2003), Penalized Pixel-Fitting (/p.pcPXF, Cappellari & Emsellem 2004; Cappellari 2017, 2023), line-strength indices (e.g., Worthey 1994; Schiavon 2007; Thomas et al. 2011; Martín-Navarro et al. 2018), or a combination of them (e.g., the GIST pipeline, Bittner et al. 2019). The results can be compared directly with those from IFS observations of MWAs in terms of mass- or light-weighted parameter maps, line-of-sight velocity distribution (LOSVD), and mass fraction distributions. The ingredients, flexibility, procedures, and computational expenses of this code are described in detail in the following sub-sections.', '2.1.1 Galactic chemical evolution (GCE) model': 'We apply the analytical chemodynamical model of the MW developed by Sharma et al. 2021b (hereafter S21) which can predict the joint distribution of position ( x ) , velocity ( v ) , age ( 𝜏 ) , extinction 𝐸 ( 𝐵 -𝑉 ) , photometric magnitude in several bands and chemical abundances ( [ Fe / H ] , [ 𝛼 / Fe ] ) of stars in the MW. Compared with previous models (e.g., Chiappini et al. 1997; Haywood et al. 2019), this model included a new prescription for the evolution of [ 𝛼 / Fe ] with age and [ Fe / H ] and a new set of relations describing the velocity dispersion of stars (Sharma et al. 2021a). This model shows for the first time that it can reproduce the [ 𝛼 / Fe ] -[ Fe / H ] distribution of APOGEE observed stars (Hayden et al. 2015) at different radius 𝑅 and height | 𝑧 | across the Galaxy. In this model, the origin of two [ 𝛼 / Fe ] sequences is due to two key processes: the delay between the first star formation and sequential occurrence of SN Ia causes the sharp transition from high-[ 𝛼 / Fe ] to low-[ 𝛼 / Fe ] at around 10.5 Gyr ago, which is likely to create a valley between the two sequences; the radial migration or so-called churning is responsible for the large spread of the low-[ 𝛼 / Fe ] sequence along the [ Fe / H ] axis. This model also showed that without churning it is not sufficient to reproduce the two sequences (see their Fig. 6). Because this chemical evolution model is purely analytical, it is easy to be inserted into the forward-modeling tool E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc to generate mock stellar catalogs for further analysis. In addition, several free parameters such as radial migration and kinematic heating efficiency are adjustable, which will be useful to implement similar analysis on external galaxies.', '2.1.2 E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc': 'To mock the MW IFS data cube, we need a comprehensive stellar catalog from observations with well-measured parameters such as position ( x ) , velocity ( v ) , age ( 𝜏 ) and chemical abundances. However, this is impractical as the MW has hundreds of billions of stars being unexplored and the dust in the disk obscures distant light. An  \nalternative way is to employ N-body/hydrodynamical simulations (e.g., EAGLE Schaye et al. 2015, FIRE-2 Hopkins et al. 2018), but most of the simulations only contain ∼ 10 6 stellar particles, which is not enough to generate integrated spectra because each spatial bin would only contain less than a hundred particles. This in turn would increase the sampling noise of the integrated spectrum and make observables derived from spectra noisy. Even though oversampling can solve this problem, it is still challenging to find a simulation that is observably identical to the MW in all aspects quantitatively, especially the [ 𝛼 / Fe ] bimodal trends.  \nTherefore, here we generate a mock catalog to represent the MW. We use E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc (Sharma et al. in prep.), which is a tool in accordance with the chemical evolution model of S21 and can create a catalog with the user-defined number of stellar particles ( 𝑁 𝑝 ) containing parameters including position ( x , y , z ) , velocity ( v x , v y , v z ) , age ( 𝜏 ) , metallicity ( [ M / H ] ) and [ 𝛼 / Fe ] , and the parameter distributions are chemodynamically consistent to the MW observations. In E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc, each stellar particle is equivalent to a single-stellar population (SSP) with an initial birth mass of 1 M ⊙ , so one can directly assign spectra using SSP templates under the default unit without transferring from the mass of initial cloud to remaining stars. Other codes like TRILEGAL (Girardi et al. 2005), BESANCON (Girardi et al. 2005), and G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc (Sharma et al. 2011) can also create mock catalogs. However, compared to E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc, the underlying models of these codes lack the information of [ Fe / H ] and [ 𝛼 / Fe ] , and do not have the processes of radial migration and kinematic heating. Furthermore, the ability of E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc to control the observed properties by regulating the underlying physical process is important for future analysis of external galaxies, whose formation history and dynamical processes are expected to be different from the MW. One caveat is that SSP particles in the current version of E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc are distributed only in the disk and there is no bulge/bar, halo, or nuclear disk structure. This is because the chemical and kinematic distribution functions predicted by S21 are extrapolated from observations in the solar neighborhood. Nevertheless, particle parameter distributions in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc are still consistent with APOGEE observations in the range of 3 < 𝑅 gc < 15 kpc (Sharma et al. 2021b).', '2.1.3 Spectral libraries': 'To build a mock data cube, one important procedure is to turn particles in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc into stellar spectra based on their properties, so a spectra library is needed. Given each particle is an SSP in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc and the integrated spectrum is a sum of many stellar populations We will employ SSP spectral libraries in G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc. The SSP spectrum describes the spectral energy distribution (SED) of a stellar population with a single age, metallicity, and chemical abundance patterns. An initial mass function (IMF) is assumed, and the stellar population is evolved using a given isochrone (Conroy 2013). G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc currently supports MILES (Vazdekis et al. 2010, 2015) and PEGASE-HR (Le Borgne et al. 2004), both of which will be used in this study.  \nThe MILES SSP library (FWHM = 2 . 51 Å , 3500 Å < 𝜆 < 7500 Å ) is based on the model of Vazdekis (1999). It uses Padova 2000 (Girardi et al. 2000) or BaSTI (Pietrinferni et al. 2004, 2006) isochrones and IMF in Unimodal/Bimodal (Vazdekis et al. 1996), Kroupa Universal/Revised (Kroupa 2001) and Chabrier (Chabrier 2003) shapes. For Padova 2000 isochrones, the template grids cover 7 metallicity bins between (-2.32, 0.22) dex, 50 age bins between (0.063, 17.78) Gyr and one scaled-solar [ 𝛼 / Fe ] bin (Vazdekis et al. 2010). As for BaSTI isochrones, the template grids cover 12 metallicity bins between (-2.27, 0.40) dex, 53 age bins between (0.03, 14.00) Gyr  \nand two [ 𝛼 / Fe ] bins in 0.0 and 0.4 dex (Vazdekis et al. 2015). Since [ 𝛼 / Fe ] enhancement is essential in this study, for most of the cases we will use the 𝛼 -variable templates.  \nThe PEGASE-HR library (FWHM = 0 . 55 Å , 3900 Å < 𝜆 < 6800 Å ) is based on the Elodie 3.1 stellar library (Prugniel & Soubiran 2001; Prugniel et al. 2007). The SSP models are computed using Padova 1994 (Bertelli et al. 1994) isochrones with a Salpeter (Salpeter 1955), Kroupa, or top-heavy (Elmegreen & Shadmehri 2003) IMF. The templates grid consists of 7 metallicity bins between (-2.30, 0.70) dex, 68 age bins between (0.001, 20) Gyr, and one scaled-solar [ 𝛼 / Fe ] bin. In this work, the PEGASE-HR templates are mainly used to explore the effect of spectral resolution on /p.pcPXF fitting. Therefore, it is still useful even if it lacks variable [ 𝛼 / Fe ] . Revised grids interpolated by ULYSS (Koleva et al. 2009) in the same way as Kacharov et al. (2018) for PEGASE-HR are also available. The new grids contain 15 steps in metallicity between -2.3 and 0.7 dex and 50 steps in age between 10 Myr and 14 Gyr.', '2.2 Configurations of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc': "To meet different research purposes, we incorporated a wide range of adjustable parameters in G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc. Specifically, the adjustable parameters are divided into three groups: GALPARAMS , SSPPARAMS and INSPARAMS , as listed with descriptions in Table 1. The user can set up their preferred parameters to obtain the expected data cubes. GALPARAMS is for setting the distance, inclination, and position of the mock MW using coordinates transformation; SSPPARAMS governs the spectral properties, such as the selection of SSP model, single or variable [ 𝛼 / Fe ] , spectral resolution, age, [ M / H ] and wavelength ranges and interpolation method to be used to assign the spectrum according to particles parameters (see details in Section 2.3). Many interpolation options could be used such as 'nearest', 'linear', etc. The INSPARAMS controls the instrument spatial sampling, atmospheric effects (PSF), and the number of spatial bins in each coordinate. We also provide an option to derive the data cube in the format of a specified instrument. Alternatively, the user can also design a hypothetical instrument that does not exist by manually setting up these parameters, which will be useful for future instrument designs.", '2.3 Procedures of making a mock data-cube': "The procedures of generating mock MW data cubes are described in detail as the following steps, along with a flowchart in Fig. 1:  \n- (i) Use E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc to generate a mock stellar catalog of the MW, with particles containing position ( x ) , velocity ( v ) , age ( 𝜏 ) , extinction 𝐸 ( 𝐵 -𝑉 ) , metallicity ( [ Fe / H ] ) and [ 𝛼 / Fe ] , transfer [ Fe / H ] to [ M / H ] .\n- (ii) Loadtheconfigurations of Table 1 provided by the user, as well as a list of data cubes to be generated with their center coordinates in ( 𝑙, 𝑏 ) or ( 𝑟𝑎, 𝑑𝑒𝑐 ) .\n- (iii) Apply coordinates transformation based on GALPARAMS parameters to project this Galaxy on the Celestial sphere with a certain distance and inclination.\n- (iv) Apply the spatial binning based on IFS instrument properties given by INSPARAMS , and find the particles included in each bin, then note the bin index for each particle for later use.\n- (v) Load the SSP templates with defined age and [ M / H ] ranges by SSPPARAMS , then oversample the spectra by a factor of SSPPARAMS:factor using spline interpolation with the order of  \nthree. Next, construct the interpolator to be used to generate integrated spectra in the next step.  \n- (vi) Select particles in a spatial bin, assign each particle an SSP spectrum based on its age, [ M / H ] , and [ 𝛼 / Fe ] by interpolating the templates with the method defined by SSPPARAMS:spec\\_interpolator . Multiply the spectrum by the particle's stellar mass because SSP spectral templates are normalized to 1 M ⊙ . Shift the spectrum due to its line-of-sight velocity ( 𝑉 LOS ) using the Doppler equation. Then apply a flux calibration due to the particle's distance and extinction.\n- (vii) Loop the procedure (vi) over all the spatial bins to stack integrated spectra. In the meantime, generate the mass/light-fraction distribution of each spatial bin. The light weights are obtained within the wavelength range given by SSPPARAMS:wave\\_range .\n- (viii) After generating integrated spectra for all the spatial bins, apply the atmosphere effects by convolving each wavelength slice with a point spread function (PSF), this can be either a Gaussian or Moffat kernel function with the given INSPARAMS:FWHM\\_spatial .\n- (ix) Degrade the stacked spectrum to the instrument resolution given by SSPPARAMS:FWHM\\_gal using convolution with a Gaussian line-spread-function.\n- (x) Re-bin the oversampled flux array into the original wave grid or the user-defined wavelength interval and wavelength range, according to SSPPARAMS:dlam and SSPPARAMS:wave\\_range .\n- (xi) Write the data cube flux array in the format of FITS file with G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc configuration parameters in the header.  \nOther than the above procedures, there are a few points that need to be clarified to the users as follows:  \n- · The original mock stellar catalog in step (i) should have two coordinate systems, Cartesian coordinates ( 𝑥, 𝑦, 𝑧 ) and Galactic coordinates ( 𝑙, 𝑏 ) , where ( 𝑙, 𝑏 ) are overlapped with ( 𝑦, 𝑧 ) , respectively. Therefore, by adjusting ( 𝜃 𝑧𝑥 , 𝜃 𝑦𝑧 ) , users can rotate the mock MW into a defined inclination. This is convenient when compared with real observations.\n- · The oversampling in step (v) has two purposes: one is to ensure the validity of degrading - when degrading the SSP templates from the original spectral resolution to instrumental resolution (e.g., from MILESto MUSE), the 𝜎 of the Gaussian kernel to be convolved with the templates will be smaller than wavelength interval Δ 𝜆 without oversampling. In this case, the kernel function array becomes invalid with only one element having a value of 1, and the rest having values of 0. Then the degraded spectrum will be still in its original resolution. Another reason is to minimize the particle sampling error when stacking the spectrum with different 𝑉 LOS . Fig. 2 shows an example mock integrated spectrum generated by using the non-oversampled (in red and light-red) and oversampled (in blue and light-blue) MILES SSP templates. The light-color lines are the spectra convolved with a Gaussian kernel by the given mean velocity and dispersion ( 𝜇, 𝜎 ) , which is the manner of /p.pcPXF during the kinematics measurements. The dark-color lines are stacked from 2000 spectra shifted by Gaussian distributed 𝑉 LOS using the same ( 𝜇, 𝜎 ) , which is the manner of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc. By calculating the residuals of these two lines (shown in the bottom panel), we find that oversampling can reduce the deviation between the convolved and stacked spectrum significantly by a factor of ∼ 25, which will be helpful to reduce potential issues when /p.pcPXF performs Gauss-Hermite convolution to fit on the particle-stacked spectrum.\n- · This package can select a portion of particles within the field of view (FoV) based on the user-defined instrument to generate the data cube, rather than taking a whole catalog into account. It helps  \nFigure 1. A flowchart demonstrating the workflow of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code. The numbers on some processes correspond to the steps listed in Section 2.3.  \n<!-- image -->  \nFigure 2. Comparison of the mock integrated spectrum generated by using the non-oversampled (in red and light-red) and the oversampled (in blue and light-blue) MILES templates. The light-color lines are spectra convolved by a Gaussian kernel with given ( 𝜇, 𝜎 ) . The dark-color lines are spectra stacked by 2000 shifted spectrum with different line-of-sight velocities which follow a Gaussian distribution with the same ( 𝜇, 𝜎 ) . The bottom panel shows residuals of convolved and stacked spectra. We find that oversampling can reduce the deviation between the convolved and stacked spectrum significantly by a factor of ∼ 25.  \n<!-- image -->  \nto mimic a more realistic IFS observation and reduce computational expenses.  \n- · G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc can be executed in the batch mode where users can provide a list of cubes with the central coordinates in ( 𝑙, 𝑏 ) or ( 𝑟𝑎, 𝑑𝑒𝑐 ) , and all the cubes can be automatically generated in one execution;\n- · Other than the data cube FITS file, this package also generates some by-products including mass/light-fraction and number of particle distribution arrays for each spatial bin, mass/light-weighted [ M / H ] , [ 𝛼 / Fe ] , age and kinematic properties. These by-products are obtained from properties of particles in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc and can be used to calculate the expected true answers that /p.pcPXF is expected to recover from full-spectrum fitting. This will be described in detail in Section 3.", '2.4 Estimate the sampling error': "Compared to the real catalogs of the MW (e.g. Gaia), one shortcoming of the E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc mock stellar catalog is the limited number of particles. Although there are 10 -100 times more particles ( 10 8 ) in the mock stellar catalog compared to MWAs in Nbody/hydrodynamical simulations, some spatial bins inevitably contain too few particles for robust measurements of their properties. Even for spaxels or Voronoi bins containing more than ∼ 10 4 particles, the spectral noise due to finite star particles is still nonnegligible. We call this type of spectral noise 'sampling error' ( 𝑒 𝑓 ,𝑆 ) . One way to reduce sampling error is to generate more particles from E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc. However, this will increase the computational expenses and RAM/memory usage of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc significantly, and exceed most of the HPCs' limitation ( ∼ 22 GB for a catalog containing 10 8 particles, hence ∼ 220 GBfor 10 9 particles). Therefore, when mocking IFS observations, one has to ensure that for each integrated  \nFigure 3. Thesampling signal-to-noise ratio (SNR) of an integrated spectrum as a function of the number of particles in a spatial bin. The shaded region is the standard deviation of sampling SNRs for spaxels with the same number of particles. This is calculated by using original and bootstrapped mock data cubes described in Section 3.1. This figure provides an estimation of spectral sampling noise due to the limited number of particles in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc.  \n<!-- image -->  \nspectrum, the sampling error ( 𝑒 𝑓 ,𝑆 ) is smaller than the observational flux error ( 𝑒 𝑓 ,𝑂 ) at given instrument conditions and exposure time. In this way, it is safe to apply the data reduction pipeline on this mock data cube and allow /p.pcPXF to derive mathematically reasonable results. This is particularly important in kinematics because the LOSVD effect in G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc integrated spectra is implemented by stacking individual spectra of particles with each shifted by its own 𝑉 LOS . But /p.pcPXF employs the Gauss-Hermite function to convolve with SSP templates and determines the kinematics moments. In Section 3.2, we will provide a detailed example of how to deal with the sampling error for each Voronoi bin.  \nTo estimate the sampling error, the G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code provides an option to apply bootstrapping. First, it randomly re-selects particles from the original E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog and generates a bootstrapped catalog with the same particle number ( 𝑁 𝑝 ). Then, G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc uses this catalog to generate bootstrapped data cubes by repeating the above procedures a certain number of times, and the sampling error ( 𝑒 𝑓 ,𝑆 ) for each spaxel can be represented by the standard deviation of these bootstrapped fluxes. This sampling flux error will be used as the lower limit of the acceptable Gaussian noise when mimicking the real observations and obtaining the final integrated spectra of Voronoi bins. Fig. 3 illustrates how the sampling SNR (flux divided by sampling noise) varies as a function of the number of particles in a spatial bin. The shaded region is the standard deviation of sampling SNRs for spaxels having the same number of particles. We obtain this figure by using the original and bootstrapped mock data cubes later in Section 3.1. It can be seen that a spatial bin with ∼ 10 3 particles can generate a spectrum with sampling SNR ∼ 25 pix -1 , and a spatial bin with ∼ 10 4 particles can generate a spectrum with sampling SNR ∼ 80 pix -1 .", '2.5 Computational expenses and multiprocessing strategy': "The execution time of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc to generate mock data cubes depends mostly on the number of particles included in the instrument FoV and the spectral interpolation method. In general, the execution time is proportional to the number of particles used, and the 'nearest' interpolation is three times faster than the 'linear' interpolation. In this code, we apply python-multiprocessing techniques to speed  \nup the computing. For a typical MUSE FoV containing 6 × 10 6 particles, the execution time spent with a 24-core CPU (2.50GHz) is ∼ 1 . 4 hour. Based on this, the users can roughly estimate the execution time they will spend. If taking into account the bootstrapped cubes, the total execution time will be multiplied by the number of bootstrapping times plus 1. Therefore, it is highly recommended to run it on high-performance computers (HPC) or a Cluster where these 21 jobs can be executed simultaneously.", '3 RECOVERY OF THE GALAXY CHEMODYNAMICAL PROPERTIES': 'In this section, we take /p.pcPXF as a representation software to test the ability of full spectrum fitting methods to recover galaxy properties by applying it to mock cubes generated from G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc. We measure kinematics ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ), stellar population parameters (age, [ M / H ] , [ 𝛼 / Fe ] ) and light/mass fraction distributions of different structural components. The analysis is performed in the same way as extragalactic studies. Then we compare the results with the input true values that are obtained by properties of SSP particles from E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog. This test allows us to access the consistency of parameters measured via broadly applied software in other studies (e.g. Pinna et al. 2019a,b; Scott et al. 2021; Martig et al. 2021), which was not possible previously as the true values of external galaxies are unknown. Furthermore, it also provides standard references for the future to better understand extragalactic results (e.g., gradient, flaring) by distinguishing real distributions from artificial effects due to the spectral fitting methods, projected view, and integrated light. We note that our goal is to explore the general performances of the full-spectrum fitting method using template weighting and regularization, i.e., its underlying mathematical framework. We choose to test /p.pcPXF because it has been widely used by many previous studies. However, any systematic bias found in this study is not specific to /p.pcPXF, but would be equally applicable to other software using the same framework. And we leave testing on other software for future studies.', '3.1 Mock cube generation for MUSE instrument': "We generate a mock MUSE observation by G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc, using the E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog that contains 10 8 particles. We remove particles with stellar age less than 0.25 Gyr because their position and kinematics are erroneous in the current version of E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc, and we confirm that removing these particles does not affect our conclusions. The mock MW catalog is assumed to have a distance of 26 . 5 Mpc and inclination of 86 · to the observer, which is the same as the projection of NGC 5746 observed by MUSE with comprehensive analysis in Martig et al. 2021 (hereafter M21). We use MILES 𝛼 -variable SSP templates (Vazdekis et al. 2015) with the BaSTI stellar isochrones (Pietrinferni et al. 2004) and Kroupa Universal IMF (Kroupa 2001). The templates have 53 bins in age, 12 bins in [ M / H ] and 2 bins in [ 𝛼 / Fe ] and we apply a 'linear' interpolation to assign each particle a spectrum based on its age, [ M / H ] and [ 𝛼 / Fe ] , and then degrade the stacked spectra to MUSE spectral resolution (FWHM ∼ 2 . 65 Å ). Following procedures in Section 2.3, we obtain two mock cubes focusing on the central ( 𝑁 𝑝 = 61575676 ) and the disk ( 𝑁 𝑝 = 7379847 ) regions, as shown in Fig. 4. This observation strategy is also the same as MUSE observations on NGC 5746 in M21. The total execution time spent by G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc on a 24-core CPU for these two cubes is ∼ 14 . 5 hours. We also generate 2 × 20 bootstrapped cubes and use 16th and 84th percentiles to calculate the sampling error of each spaxel. We  \nFigure 4. Two mock data cube FoV (in purple) relative to the number distribution of particles in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc mock stellar catalog. The model was set with a distance of 26.5 Mpc and inclination of 86 · . This observation strategy is the same as the MUSE observation on NGC 5746.  \n<!-- image -->  \ndo not apply extinction in these cubes because here we only focus on full-spectrum fitting validation. Adding extinction would blend all the effects and make it difficult to differentiate their individual impacts. Therefore, we reserve this topic for future studies.  \nThe next procedure is to add Gaussian flux error to the spectra. We first derive the observational flux error ( 𝑒 𝑓 ,𝑂 ) of the mock cubes. The observational flux error depends on many aspects but can be classified into two main categories: the observation conditions (seeing, airmass, exposure time, etc.) and the instrumental properties (telescope aperture, system efficiency, dark current, read-out noise, etc). For simplicity, we ignore the sky conditions, dark current, and read-out noise which only contribute a few percent to total received photons, and assume the spectral SNR and received photons are defined by  \nSNR = 𝑓 𝑒 𝑓 ,𝑂 = √ 𝑁, 𝑁 = 𝑎 𝑓 𝑡, (1)  \nwhere 𝑓 is the flux of the target; 𝑒 𝑓 ,𝑂 is the observational flux error (by MUSE in this case); 𝑁 is the received number of photons; 𝑡 is the exposure time; 𝑎 is an overall reaction of sky transmission, efficiency and telescope aperture. Therefore, 𝑎 should only depend on wavelength for the same instrument. By substituting the above equations, the parameter 𝑎 can be calculated using 𝑓 , 𝑒 𝑓 ,𝑂 and 𝑡 from an observation by  \n𝑎 = 𝑓 𝑒 2 𝑓 ,𝑂 × 𝑡 . (2)  \nIn this work, we take all the bulge and disk observations of NGC 5746 from MUSE in M21 and fit equation 3 as a function of wavelength ( 𝜆 ) using a 4-degree polynomial, which is described by  \n𝑎 ( 𝜆 ) = 4 . 34274826 -18 𝜆 4 -1 . 43263443 -13 𝜆 3 + 1 . 61141240 -9 𝜆 2 -7 . 29164505 -6 𝜆 + 1 . 17213217 -2 , (3)  \nwhere 𝑎 ( 𝜆 ) is in the unit of 1 / 10 -20 erg cm -2 Å -1 . Next, we set the bulge and disk mock cubes to have an exposure time of 1729.39 s and 6221.84 s, respectively, and use the equation 2 and 3 to estimate the observational flux error 𝑒 𝑓 ,𝑂 of each spaxel. Then we use this error to add Gaussian noise to all the spaxels. Here the disc exposure time is chosen to satisfy the upper bound of 𝑒 𝑓 ,𝑂 ≥ 𝑒 𝑓 ,𝑆 . The bulge exposure time is then determined by assuming the same bulge-to-disk exposure time ratio in MUSE observations of NGC 5746 by M21, and we confirm that 𝑒 𝑓 ,𝑂 ≥ 𝑒 𝑓 ,𝑆 also applies. These two values are slightly smaller than those used by M21. Finally, the two mock cubes are stitched together.  \nFigure 5. A spectral fit example of a Voronoi bin by /p.pcPXF in Section 3.2. The 'observed' spectrum is in black while the best-fit spectrum is in red. Residuals are in green. The fitting area is within the grey lines on both sides.  \n<!-- image -->", '3.2 Extracting galaxy properties': "WeapplytheGISTpipeline 1 (Bittner et al. 2019) on the stitched mock cube to measure the kinematics and stellar population parameters. The GIST pipeline combines all the tools needed to process the data and users can obtain final results in a single execution. Here we use a modified version to implement some functionalities that the current public version (v3.1.0) does not have but are needed in this work. A detailed list of added features is given in Appendix A.  \nWerunthe GIST pipeline in the following steps: First, we apply the Voronoi tessellation software (Cappellari & Copin 2003) to spatially re-bin the mock cube and increase the spectral SNR to 80 pix -1 , which results in 1477 Voronoi bins. Most of the Voronoi bins contain 𝑁 𝑝 > 10 4 . For the other Voronoi bins, 𝑁 𝑝 is also very close to this number. To ensure the sampling noise 𝑒 𝑓 ,𝑆 is less than the MUSE flux error 𝑒 𝑓 ,𝑂 , we apply the same Voronoi binning arrangement to all 20 bootstrapped cubes. Then we calculate the sampling noise 𝑒 𝑓 ,𝑆 of each Voronoi bin by taking half of the difference between the 16th and 84th percentiles of its 20 bootstrapped integrated spectra. We confirm that on average 𝑒 𝑓 ,𝑂 > 𝑒 𝑓 ,𝑆 applies for all the bins.  \nNext, we apply the full-spectrum fitting method /p.pcPXF (Cappellari & Emsellem 2004; Cappellari 2017, 2023) to measure the stellar kinematics of each Voronoi bin. Here we use the same MILES 𝛼 -variable templates (Vazdekis et al. 2015) as when we generated the mock cubes in the same MUSE spectral resolution. In the following analysis, we always use the whole grid of MILES templates in age, [ M / H ] , and [ 𝛼 / Fe ] . We are aware of some unsafe regions as listed in Table. 2 of Vazdekis et al. (2015). However, because no particles fromtheE-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog have ages or [ M / H ] that fall within these regions, it does not affect our conclusions. We also confirm that our conclusions remain the same after removing templates in the unsafe regions. Since there is no emission line or atmosphere effect, we use a wide wavelength range of ( 4750 , 6550 ) Å to fit with the Voronoi binned spectra and remove the first and last 75 Å to avoid effects caused by spectral oversampling, rebinning and Doppler shift, etc. During the fitting, the MILES templates are convolved with a line-ofsight velocity distribution (LOSVD) described by a Gauss-Hermite equation to match the Voronoi binned spectrum. We parameterize the LOSVD using four moments, which are mean line-of-sight velocity 𝑉 LOS , line-of-sight velocity dispersion 𝜎 , and the third- and fourthorder moments ℎ 3 , ℎ 4 . No regularization is applied in this process,  \nand we include a fourth-order multiplicative Legendre polynomial during the fitting.  \nAfter measuring kinematics, we employ /p.pcPXF again to obtain the stellar population parameters for each Voronoi bin. We choose the same templates, spectral resolution, and fitting wavelength range as the previous step. We use the LOSVDs ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ) measured in the last step as input, and fix them during the fitting to obtain the weight of each template. The best-fit spectrum is the weighted sum of all the templates. For the initial fitting, we set no regularization to obtain the initial 𝜒 2 . To avoid the ill-conditioned inverse problem due to severe degeneracies in mathematics (as discussed in Section 3.5 of Cappellari 2017), we then follow the approach in McDermid et al. (2015) and iterate the fitting by increasing regul until 𝜒 2 is increased by Δ 𝜒 2 = √ 2 𝑁 , where 𝑁 is the number of wavelength pixels considered for the fit. This iteration process allows us to obtain the smoothest solution that is still compatible with the data within 1 𝜎 level. At this stage, we note regul reaches the maximum regularization regul 𝑚𝑎𝑥 . Next, we choose regul = 5 which is between 0 and regul 𝑚𝑎𝑥 ( 30 -100 ) to keep smooth solutions while still allowing for a variation on short timescales of the star formation, which will disappear if regul is too large (see similar discussions in Pinna et al. 2019a and M21). Here we only apply the first order of regularization for all the fittings. We employ this fitting strategy using templates normalized to both 1 M ⊙ and 1 L ⊙ to obtain mass- and light-weighted results, respectively. Fig. 5 is an example of the fitting on one Voronoi binned spectrum. Finally, using the weights, we can calculate light-weighted and mass-weighted age, [ M / H ] , [ 𝛼 / Fe ] , and light/mass fraction distributions.  \nTo estimate uncertainties of template weights and weighted-mean age, [ M / H ] , and [ 𝛼 / Fe ] , we apply Monte Carlo (MC) realizations similar to the way described in Kacharov et al. (2018), M21, and Cappellari (2023). For each Voronoi-binned spectrum, we first perform a fitting with no regularization and obtain the best-fit spectrum and residuals compared to the initial input spectrum. Next, we add residuals with signs randomly flipped to the best-fit spectrum 100 times to obtain resampled spectra. Each resampled spectrum is then fitted by /p.pcPXF with regul = 0 . 1 to reduce bias towards smooth solutions. Finally, we calculate uncertainties of template weights and weighted mean properties by taking the half of difference between the 16th and 84th percentiles of 100 MC realizations.  \nIn addition, we also apply the LS module in the GIST pipeline to compute line-strength indices of each Voronoi bin spectrum in the LIS system (Vazdekis et al. 2010) choosing the definitions provided at a spectral resolution of 8.4 Å. This routine was presented by Kuntschner et al. (2006) and Martín-Navarro et al. (2018). Next, given the relationship between templates' properties (age, [ M / H ] , [ 𝛼 / Fe ] ) and line strengths, the measured line strengths can be matched to SSP-equivalent parameters by using the MCMC implementation from the package emcee (Foreman-Mackey et al. 2013). In this work, we follow M21 and use H 𝛽 o as an age indicator and Fe5015, Fe5270 and Fe5335, and Mg 𝑏 to trace metallicity and [ 𝛼 / Fe ] , respectively. The Monte-Carlo simulation is run 15 times and each one uses 100 walkers and 1000 iterations to obtain uncertainties.  \nAfter measuring all the parameters, in the next sub-sections, we compare the kinematics, light- and mass-weighted stellar population parameters, and light/mass fraction distributions of the mock cube with the true values to verify the recovery ability of full-spectrum fitting method /p.pcPXF.  \nFigure 6. Stellar kinematics maps ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ) of the mock MUSE cube generated in Section 3.1. The scale of the color bar is given in the second line of the upper left corner for each panel. Top row : True values calculated by fitting a Gauss-Hermite equation with particles' velocity distribution weighted by their total flux in /p.pcPXF fitted wavelength region for each Voronoi bin. Middle row : Results from /p.pcPXF. Bottom row : Residuals of the /p.pcPXF results and the true values. The residual panels indicate results from /p.pcPXF have systematic offsets compared to true values.  \n<!-- image -->", '3.3 Kinematic maps': 'Fig. 6 shows the kinematics maps in four moments ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ) of the mock MUSE cubes. The scale of the color bar is given in the second line of the upper left corner of each panel. We calculate the true values shown in the first row with the following procedures: First, for each Voronoi bin, we calculate the total flux of each particle in the /p.pcPXF fitted wavelength region. Then we plot flux-weighted 𝑉 LOS histogram distribution using all the particles included in this bin. Next, we fit this histogram with a Gauss-Hermite equation and obtain four best-fit moments ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ). This method is consistent with the definition of light-weighted kinematics which /p.pcPXF is expected to recover during the fitting process. The second row is the results from /p.pcPXF and the bottom row shows the residuals between /p.pcPXF results and the true values.  \nIn this figure, the kinematic moments obtained by /p.pcPXF have the same trend compared to the true values. Both show two different structures: one is aligned to 𝑦 ∼ 0 and the thickness increases with 𝑥 , which has larger absolute 𝑉 LOS and ℎ 3 , and smaller 𝜎 ; the other is in a similar projected radius but vertically higher and thicker, and it has smaller absolute 𝑉 LOS and ℎ 3 , but larger 𝜎 . An anti-correlation of ℎ 3 with 𝑉 LOS which are usually associated with disk-like components (e.g. Krajnović et al. 2008; Guérou et al. 2016; van de Sande et al. 2017) is seen and are similar to MUSE edge-on galaxies studies of Pinna et al. (2019a,b) and M21. These two components are mostly likely to be associated with thin and thick disks. We will explore this in detail in Section 3.6.  \nHowever, in the residual panels, all these four moments show systematic offsets. Compared to the true values, 𝑉 LOS from /p.pcPXF is around 17 km s -1 lower above and below the very thin midplane ( 𝑦 ∼ 0) and shows a more significant difference around 𝑥 ∼ [ 10 , 30 ] arcsec. Around the galaxy center, the residual of 𝑉 LOS also shows a continuous decrease from negative to positive 𝑥 ; 𝜎 is generally overestimated everywhere in the galaxy with few light blue residuals. ℎ 3 is overestimated in regions of 𝑥 ∼ [ 10 , 60 ] arcsec and 𝑦 ∼ [ 10 , 25 ] arcsec and underestimated in the outer region of 𝑥 ∼ [ 60 , 110 ] arcsec; ℎ 4 from /p.pcPXF has no significant structures like  \n𝜎 and ℎ 3 maps, which is also seen in real galaxies results (e.g., Pinna et al. 2019a,b and M21), but the true ℎ 4 map clearly shows kinematic differences. The clear structures in these residual panels indicate that it is not because of the fitting uncertainties. We will investigate this in detail in Section 4.1.', '3.4 Stellar population property maps': "Fig. 7 and Fig. 8 show the light- and mass-weighted age, [ M / H ] and [ 𝛼 / Fe ] maps of the mock MUSE cubes, respectively. The first row is the true values by calculating the median age, [ M / H ] , and [ 𝛼 / Fe ] of E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc particles weighted by luminosity or mass, which are equivalent to light- or mass-weighted values. The second and third rows are results from /p.pcPXF with regul = regul 𝑚𝑎𝑥 and regul = 5, respectively. We also write the average uncertainty from MC realizations in the top left corner of the third row. The last row shows the residuals of the /p.pcPXF results with regul = 5 and the true values, with the average of absolute residual written in the top left corner. For both figures, the overall distributions of these three parameters obtained by /p.pcPXF are very close to the true values, and the residuals are within the order of uncertainties on average. This confirms the reliability of spectral fitting methods to measure the weighted age and chemical compositions. Especially, the [ 𝛼 / Fe ] map from /p.pcPXF with regul = 5 indicates the capability of /p.pcPXF to identify distinct [ 𝛼 / Fe ] -rich and [ 𝛼 / Fe ] -poor populations in the thick and thin disk, respectively, even though only two [ 𝛼 / Fe ] bins are available. The residuals of [ 𝛼 / Fe ] from /p.pcPXF with regul = 5 and the true values are flat and no systematic pattern is found.  \nHowever, mass-weighted age from /p.pcPXF is slightly overestimated in the outer regions with more yellow and red Voronoi bins in the residual panel, where the residuals are larger than uncertainties. This overestimation is much more obvious in mass-weighted results. In addition, even though residuals of light-weighted [ M / H ] are mostly close to 0, the mass-weighted [ M / H ] are overestimated in the central regions. This means the age gradient from /p.pcPXF is underestimated but [ M / H ] gradient is overestimated, and such effects are more dom-  \nFigure 7. Light-weighted stellar population property maps (age, [ M / H ] and [ 𝛼 / Fe ] ) of the mock MUSE cube generated in Section 3.1. The scale of the color bar is given in the second line of the upper left corner for each panel. First row : True values calculated by the light-weighted average of particles' age, [ M / H ] and [ 𝛼 / Fe ] for each Voronoi bin. Second row : Results from /p.pcPXF with regul = regul 𝑚𝑎𝑥 , calculated using strategies of McDermid et al. 2015. Third row : Results from /p.pcPXF with regul = 5 and the average uncertainty of all the Voronoi bins from MC realizations is written in the left top corner. Last row : Residuals of /p.pcPXF results with regul = 5 and the true values. The average absolute residual of all the Voronoi bins is written in the left top corner. This figure indicates that /p.pcPXF results with proper regularization can identify different galaxy components by their stellar population parameters, which are consistent with true values and the residuals are within the order of uncertainties on average. When applying regul = regul 𝑚𝑎𝑥 , the distributions are smoothed and the [ 𝛼 / Fe ] panel becomes inconsistent with the true values.  \n<!-- image -->  \ninant in mass-weighted results. In addition, the [ 𝛼 / Fe ] distribution from /p.pcPXF results with regul = regul 𝑚𝑎𝑥 is almost uniformly high and much larger than the true values for all the Voronoi bins in both light- and mass-weighted results. This is because when regul is very large, the /p.pcPXF algorithm forces the result to have very smooth template weights in three parameter dimensions (age, [ M / H ] , [ 𝛼 / Fe ] ). Since there are only two [ 𝛼 / Fe ] grids, regularization will force them to have similar weights to achieve smoothness requirements and does not permit large deviations (e.g., more than 2%). Therefore, it will be challenging to identify [ 𝛼 / Fe ] bimodality. Results from /p.pcPXF with regul = regul 𝑚𝑎𝑥 also show much underestimation for age gradients along the x-axis than results with regul = 5. The age and [ M / H ] gradients are essential properties to help understand the star formation and chemical enrichment processes. Therefore, a wrong choice of regularization will then easily lead to wrong conclusions. We will explore the reasons for these offsets in more detail in Section 3.6 using light and mass fraction distributions and the effect of regularization in Section 4.2 and Section 4.3.", '3.5 SSP-equivalent maps from line-strength indices': 'Fig. 9 shows the SSP-equivalent age, [ M / H ] and [ 𝛼 / Fe ] maps of the mock MUSE cubes measured by line-strength indices. This figure shows that the main structures we derived from /p.pcPXF are also recovered by the line-strength analysis with consistent trends. In the age panel, young populations are closer to the mid-plane, and old populations are further to the mid-plane or above/below the central region. In the [ M / H ] panel, we see the metallicity gradient from the inner center to the outer Galaxy. In the [ 𝛼 / Fe ] panel, we could see the 𝛼 -rich bins in the center and 𝛼 -poor bins in the outer region, even though the differences are not as obvious as /p.pcPXF results in Fig. 7 and Fig. 8. The main difference compared with /p.pcPXF results is that the age panel shows a very low range of [ 1 -5 ] Gyr. This is also seen in M21 (Fig. B1) and because of the Balmer line indices being dominated by young stars. Therefore, the SSP-equivalent ages only reflect the fraction of stars formed within the past Gyr (Serra & Trager 2007; Trager & Somerville 2009). The SSP-equivalent [ M / H ] and [ 𝛼 / Fe ] range are much closer to /p.pcPXF results because young populations do not contribute much to the metal lines, which is also indicated in M21. This figure confirms that both line-strength indices and /p.pcPXF analysis can identify 𝛼 -rich and 𝛼 -poor populations.  \nFigure 8. Same as Fig. 7 but for mass-weighted stellar population properties.  \n<!-- image -->  \nFigure 9. SSP-equivalent age, [ M / H ] and [ 𝛼 / Fe ] maps of the mock MUSE cube generated in Section 3.1. The parameters are measured by line-strength indices. The scale of the color bar is given in the second row of the upper left corner for each panel. Similar to Fig. 7, galaxy components with different ages, [ M / H ] , and [ 𝛼 / Fe ] can be identified. However, the parameter ranges differ from /p.pcPXF results and the true values.  \n<!-- image -->', '3.6 Weight fraction distributions of different galaxy components': "In addition to calculating the light- and mass-weighted properties, we can also study the light/mass fraction distribution of stellar populations along the age and [ M / H ] dimension. This is done by using weights of templates from /p.pcPXF. Because the flux of each template is normalized to either 1 L ⊙ or 1 M ⊙ , the weights array from /p.pcPXF outputs in our tests are equivalent to stellar population light or mass fractions. Therefore, we can study the weight distribution of any component of the mock MW. M21 employed multiple components morphological fitting to a Spitzer 3.6 micron image of NGC 5746 to obtain regions dominated by the boxy/peanut bulge, nuclear disk, and thin and thick disks. In Fig. 10, we artificially select similar regions based on the locations ( 𝑥, 𝑦 ) of different components of M21, and  \nname them 'upper central', 'inner central', 'thin disk' and 'thick disk', as shown in different colors. We call them 'upper central' and 'inner central' because there is no boxy/peanut bulge and nuclear disk in the GCE of S21. Note these component definitions are purely following those in M21 to mock their data analysis. In reality, radial scale lengths of the thin disk ( 𝑅 𝑡 ) and thick disk ( 𝑅 𝑇 ) in NGC 5746 and the MW are very different ( 𝑅 𝑡 MW = 2 . 6 ± 0 . 5 kpc and 𝑅 𝑇 MW = 2 . 0 ± 0 . 2 kpc from Bland-Hawthorn & Gerhard 2016; 𝑅 𝑡 NGC 5746 = 6 . 1 kpc and 𝑅 𝑇 NGC 5746 = 8 . 2 kpc from M21). For each component, the light and mass weights of all the Voronoi bins are combined to represent its light and mass fraction distributions.  \nFig. 11 shows the light and mass fraction distributions of these four components, respectively. The total weights are normalized to one for each panel. The left two columns are light fractions and  \nFigure 10. Demonstration of the projected thin disk (blue), thick disk (red), upper central (purple), and inner central (orange) regions on top of a greyscale image of the mock data cube generated in Section 3.1.  \n<!-- image -->  \nthe right two columns are mass fractions, respectively. The /p.pcPXF results (first and third columns) are obtained with regul = 5. The true values (second and last columns) are calculated using particles' properties in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog weighted by luminosity or mass. Weonly plot the weights above 0.001 Gyr -1 dex -1 . For the thin disk, true light fractions are dominated by the youngest populations with ages less than 2 Gyr, and /p.pcPXF can consistently recover its distribution. True mass fractions indicate a rapid metallicity enrichment history ∼ 10 Gyr ago, and it slowly increases later on. However, such a metallicity enrichment trend is indistinguishable in /p.pcPXF mass fraction distributions. Moreover, /p.pcPXF mass fractions are dominated at relatively young (2 -4 Gyr) and old (12 -14 Gyr) stellar populations. The same features are also seen in mass fraction distributions of other Galaxy components. We will investigate these findings in detail in Section 4.2 to 4.5. For the thick disk, both true light and mass fractions show the metallicity enrichment trend which is again not seen in /p.pcPXF results. In addition, the thick disk contains populations that are young and more metal-poor than the thin disk. One reason is the geometrical definition of the thick disk, which contains young and relatively metal-poor stars that are flared in the outer disk ( 𝑥 ∼ [ 60 , 100 ] arcsec) due to radial gradients of age and metallicity. Given that NGC 5746 in M21 is analysed to be four times more massive than the MW in total stellar mass, it has a larger scale length and the definition of its thick disk might not apply to the MW. Another reason could be the projection effect, where young stars flared in the outer disk can have large Galactocentric (intrinsic) radius but small projected radius, so they could appear at the front and back of the line of sight in the region of 𝑥 ∼ [ 30 , 60 ] arcsec. The upper central mass fractions show a similar metallicity enrichment trend with the thick disk and is more dominant in the old populations, but this domination is smoothed out in the /p.pcPXF results. Moreover, light fractions from /p.pcPXF are dominated at the age around 8 Gyr, which is inconsistent with the true light fractions. For the inner central, true light and mass fractions are showing a clearer chemical enrichment trend and there is no new population born with [ M / H ] < -0 . 2 dex in the young region, while the /p.pcPXF mass fractions show again two overdensities at 2 -4 Gyr and 12 -14 Gyr. Therefore, except for the overestimation in the young and old population regions in mass fraction distributions, and the inconsistency in light fractions of the upper central region, the light and mass fraction distributions of /p.pcPXF are generally in agreement with true values.  \nIn Fig. 12, we integrate the mass fraction distributions in Fig. 11 along the two axes and derive age and [ M / H ] distributions for each component. The top panels are mass distributions as a function of age which is the definition of star formation history (SFH) or star  \nformation rate (SFR), and the bottom panels are defined as metallicity distribution function (MDF). Results from /p.pcPXF are in red lines and the true values are in black lines. For each panel, we calculate the correlation between these two lines to quantify their similarity. For the age distributions, we find the same as in Fig. 11. Compared to the true values, the /p.pcPXF results of all the components demonstrate an overestimation of weights in the ranges of 2 -4 Gyr and 12 -14 Gyr and underestimation in the range of 4 -11 Gyr. The underestimated regions seem to compensate for the overestimated regions. And the thin disk, inner central, and global panels show a peaked feature with age < 1 Gyr. For the [ M / H ] distributions, /p.pcPXF results are consistent with true values for most regions, as indicated by the correlation coefficients. However, mass fractions in the metal-rich region are overestimated. Other than that, the overall trend of results from /p.pcPXF is consistent with true values.  \nIn conclusion, we find that spectral fitting methods can recover the broad trends of 2-D light and mass fraction distributions for different components, but with mass fractions overestimated in 2 -4 Gyr, 12 -14 Gyr and most metal-rich regions. When integrating into 1-D age and metallicity distributions, these inconsistencies are significant. According to correlation coefficients, [ M / H ] distributions are more consistent with the true values than age distributions. We will investigate the reasons for such differences in detail in Section 4.2 to 4.5.", '3.7 Distributions of [ 𝛼 / Fe ] -[ M / H ] along different galaxy locations ( 𝑅, 𝑧 )': "According to S21, different radial migration and kinematic heating efficiency will cause different fractions of stars in the distinct [ 𝛼 / Fe ] -rich and [ 𝛼 / Fe ] -poor sequences. A recent study by Scott et al. (2021) derived this distribution from an external galaxy UGC 10738 using MILES 𝛼 -variable templates (Vazdekis et al. 2015) and concluded it has similar bimodality distributions to the MW. Given we know the true values from the mock stellar catalog particles, here we explore the recovery ability of spectral fitting methods on measuring the changes of this bimodality at different projected 𝑅 and | 𝑧 | . Different from Scott et al. (2021) where they compared integrated values of UGC10738 with individual stellar values of the MW, our integratedto-integrated value comparison on stellar population level is more direct and there is no systematic bias due to different methodologies. This direct comparison also provides an example for future studies on comparing an integrated version of the MW [ 𝛼 / Fe ] bimodality with MW-like edge-on galaxies from MUSE observations (e.g., GECKOS survey van de Sande et al. 2024).  \nThe results are shown in Fig. 13 where we separate different locations in the same way as Hayden et al. (2015). For each panel, the total mass fraction given by true (dotted lines) and /p.pcPXF (solid lines) is normalized to 1, respectively. Blue lines are mass fraction distributions for [ 𝛼 / Fe ] = 0 . 0 dex while red lines represent [ 𝛼 / Fe ] = 0 . 4 dex. The overall trends for [ 𝛼 / Fe ] -rich and [ 𝛼 / Fe ] -poorsequencesfrom/p.pcPXF are consistent with the true values for most of the panels, which indicates both sequences are well recovered by /p.pcPXF, as also shown by correlation coefficients. However, there are some discrepancies such as in the inner regions ( 𝑅 < 5 kpc) and thin disk ( | 𝑧 | <0.5 kpc), where /p.pcPXF show smoother metallicity distributions in both the blue and red lines compared to true values. In addition, the mass fraction is underestimated at -1 < [ M / H ] < 0 . 0 dex and overestimated at the most metal-rich regions, which is the same as metallicity distributions in Fig. 12.  \nWhen comparing different [ 𝛼 / Fe ] sequences, the peaks of distributions are at the same metallicity position in both /p.pcPXF and true  \nWeight fraction (Gyr  \nFigure 11. Light and mass fraction distributions of the thin disk, thick disk, inner central, and upper central of the mock data cube generated in Section 3.1. The total template weights are normalized to 1 for each panel. First column : Light fraction distributions from /p.pcPXF with regul = 5. Second column : True light fraction values calculated using particles' properties in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog. Third column : Mass fraction distributions from /p.pcPXF with regul = 5. Last column : True mass fraction values calculated using particles' properties in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog. The color bar is shown on the top. We only plot the weights with values above 0.001 Gyr -1 dex -1 . This figure indicates the broad trends of results from /p.pcPXF are consistent with the true values for the thin/thick disk and inner central regions. However, both the light and mass fraction distribution from /p.pcPXF are different from the true values for the upper central region. In addition, there is an overdensity in 2 -4 Gyr in /p.pcPXF mass distributions, and populations around 12 Gyr are smoothed towards older, more metal-rich regions, which are not seen in the true values.  \n<!-- image -->  \n1  \ndex  \n1  \n)  \nresults. This differs from our understanding of [ 𝛼 / Fe ] -[ M / H ] relation of the Milky Way. It is more likely due to the limited number of [ M / H ] and [ 𝛼 / Fe ] bins in MILES templates that cause the differences of metallicity distribution in different [ 𝛼 / Fe ] bins to be indistinguishable. Similar findings are also mentioned by Scott et al. (2021) when they analyzed the [ 𝛼 / Fe ] -bimodality of UGC 10738. Therefore, we emphasize that more [ M / H ] and [ 𝛼 / Fe ] bins in the spectral templates are necessary to obtain more detailed distributions of this bimodality. It can help to make it more feasible to identify the effect of different radial migration and kinematic heating efficiency from the relative fractions of these two sequences when compared with other IFS observations.", '4.1 Systematic offsets in the kinematics recovery': "In Fig. 6, we show the deviations between the kinematics map from /p.pcPXF and the true values. To better illustrate the offsets of each moment, we plot their residuals as a function of true velocity dispersion 𝜎 true in Fig. 14, in which the y-axis is calculated by /p.pcPXF fitted value subtracting the true value. Each data point represents one Voronoi bin. The blue dotted lines represent the instrumental dispersion ( 𝜎 inst ) based on spectral resolution. We also plot the zero-line in red to guide the eyes. The first row demonstrates the results in Fig. 6, and it clearly shows the systematic offsets for each moment: Δ 𝑉 LOS  \n1  \nFigure 12. Age distribution (top row) and [ M / H ] distribution (bottom row) of the thin disk, thick disk, inner central, and upper central of the mock data cube generated in Section 3.1. These panels are obtained by integrating the mass fraction distributions in Fig. 11 along [ M / H ] and age axis, respectively. Red lines are results from /p.pcPXF with regul = 5. Black lines are the true values calculated using particles' properties in the mock E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog. The correlation coefficient of these two lines is written in the left top corner for each panel. The overall trends from /p.pcPXF are consistent with the true values for [ M / H ] . However, same as Fig. 11, the mass fraction is overestimated in 2 -4 Gyr and 12 -14 Gyr and underestimated in 5 -12 Gyr for all the components, and the coefficients indicate inconsistency for age distributions. We also find peaked features in regions less than 1 Gyr. The mass weights of the most metal-rich bin are overestimated for all the components.  \n<!-- image -->  \n1.0  \nFigure 13. [ M / H ] distributions for different [ 𝛼 / Fe ] components distributed along projected 𝑅 and | 𝑧 | of the mock cube generated in Section 3.1. For each panel, the total mass fractions of true (dotted lines) and /p.pcPXF (solid lines) are normalized to 1, respectively. Blue lines represent mass fraction distributions for [ 𝛼 / Fe ] = 0 . 0 dex while red lines represent the distributions with [ 𝛼 / Fe ] = 0 . 4 dex. Solid lines represent the results from /p.pcPXF and dotted lines are the true values. This figure is similar to the [ 𝛼 / Fe ] -[ Fe / H ] distributions in Hayden et al. 2015 and Scott et al. 2021. The correlation coefficient of the dotted and solid lines are written in the same color in the left bottom corner. For most of the panels, /p.pcPXF results are consistent with the true values. The difference appears at the inner central and thin disk, where both lines show a mass fraction underestimation at [ M / H ] ∼ 0 . 0 dex and overestimation at the most metal-rich regions.  \n<!-- image -->  \nincreases with 𝜎 true , and 𝜎 is overestimated for most of the Voronoi bins; as for ℎ 3 and ℎ 4 , the residuals have a slight positive gradient.  \nFrom Fig. 6, we find that several Voronoi bins have large absolute ℎ 3 and ℎ 4 true values (e.g., central region and the thick disk), which will not be well constrained when SNR is relatively low (Cappellari & Emsellem 2004; Pinna et al. 2021). Therefore, /p.pcPXF with keyword bias will penalize them towards lower absolute values during the fitting to make the LOSVD towards a Gaussian shape. However, using different penalizations can result in different kinematic measurements (see full analysis using SAMI galaxies by van de Sande et al. 2017). In our case, when the LOSVD is strongly non-Gaussian, the Gauss-Hermite approximation no longer works well and large ℎ 3 and ℎ 4 could be obtained. Therefore, we turn off the penalization of these higher-order moments and re-measure the kinematics, and  \nshow the results in the second row of Fig. 14. A spatially distributed map can also be seen in Fig. B1 of the Appendix B. After turning off penalization, both ℎ 3 and ℎ 4 from /p.pcPXF are visually closer to the true values in Fig. B1. However, the second row of Fig. 14 still shows similar trends to the first row, and the four moments are more scattered, which indicates penalization is not the cause of the systematic offsets.  \nAnotherpossible reason causing the kinematic inconsistency could be the low spectral resolution of the MUSE instrument. As shown in the first row of Fig. 14, 𝜎 true values for most of the Voronoi bins are smaller than the MUSE instrumental dispersion 𝜎 inst ∼ 62 km s -1 . In this case, the recovery of kinematics will be uncertain. This has been pointed out in Fig. 3 of Cappellari (2017), and the only way to avoid it is to increase the instrumental spectral resolution. The  \nFigure 14. Residuals of four Gauss-Hermite moments Δ ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ) as a function of true velocity dispersion. Each data point represents a Voronoi bin. The red line at zero is to guide the eyes. The vertical dashed line indicates the instrumental dispersion 𝜎 inst . First row : Results from /p.pcPXF using MILES templates in MUSE spectral resolution ( 𝜎 inst ∼ 62 km s -1 ), same as the residual panels in Fig. 6. Second row : Similar to the first row but results are fitted by /p.pcPXF with no penalization on ℎ 3 and ℎ 4 . Third row : Results from /p.pcPXF on a mock cube similar to Section 3 but generated and fitted using PEGASE-HR (Le Borgne et al. 2004) templates ( 𝜎 inst ∼ 15 km s -1 ). Fourth row : Similar to the third row but penalization was turned off during /p.pcPXF fitting. Fifth row : Similar to the third row but for LOSVD-Convolved mock cubes. Sixth row : Similar to the first row but for LOSVD-Convolved mock cubes. Seventh row : Similar to the fourth row but all the analysis is in MUSE spectral resolution (FWHM ∼ 2 . 65 Å ). This figure indicates that both the spectral resolution and the relation between LOSVD and age, [ M / H ] create the kinematics systematics.  \n<!-- image -->  \nexplanation is when 𝜎 < 𝜎 inst , during /p.pcPXF fitting, the broadening of one sharp spectral feature is less than the distance to its nearby wavelength pixels. In this case, the nearby wavelength pixels have minor changes and this brings difficulties for /p.pcPXF to measure 𝜎 correctly, and ℎ 3 and ℎ 4 will go towards zero because there are also not enough pixels to identify the skewness and kurtosis.  \nTo remove the effect of low spectral resolution, we tested the kinematics recovery again by using PEGASE-HR templates (Le Borgne et al. 2004) generated by Kroupa IMF and PADOVA 1994 isochrones, which has a higher spectral resolution (FWHM = 0 . 55 Å ) than MUSE (FWHM ∼ 2 . 65 Å ). The instrumental velocity scale of PEGASE-HR is 𝜎 inst ∼ 15 km s -1 , which is smaller than the minimum 𝜎 true of our Voronoi bins. We repeat all the procedures in Section 3.2 to generate new cubes using PEGASE-HR, and then apply them to the GIST pipeline to measure the kinematics. We keep the original PEGASEHR spectral resolution throughout the whole process. The third row  \nof Fig. 14 shows the kinematics recovery using PEGASE-HR (also see the kinematics map in Fig. B2). Compared to the first row, the residuals of each panel are slightly better or clearer when using higher spectral resolution templates. However, all these four moments from /p.pcPXF still have the same systematic offsets to the true values. Most improvements are for ℎ 3 and ℎ 4 , with less bias to zero values because the higher spectral resolution helps identify skewness and kurtosis. We also show the results by turning off penalization in the fourth row (also see the kinematics map in Fig. B3), and find the ℎ 3 and ℎ 4 residuals are identical to the third row. Therefore, turning off penalization helps on low-spectral-resolution kinematics measurements, but not for higher-resolution spectra. The higher spectral resolution improves the kinematics recovery but does not fix the systematics offsets.  \nFundamentally, during the spectral fitting process, SSP templates are loaded and convolved with a Gauss-Hermite function by FFT  \nFigure 15. LOSVD of particles with different ages and [ M / H ] in one Voronoi bin. Particles are divided into three [ M / H ] and age groups, then their 𝑉 LOS histograms are plotted on the left and right panels in different colors, respectively. The grey histograms are the total 𝑉 LOS distributions of the three groups. The top rows are results using the original mock cubes. while the bottom rows use LOSVD-Convolved mock cubes. The red lines are /p.pcPXF results for mock cubes generated and fitted using PEGASE-HR templates in PEGASE-HR spectral resolution (FWHM = 0 . 55 Å ). The black lines are the true values calculated by fitting the Gauss-Hermite function with particles' light-weighted LOSVD in the black histogram. The kinematic parameters for these two lines are written in the top right corner. This figure indicates LOSVD changes with age and [ M / H ] . After using LOSVD-Convolved cubes, this relationship can be artificially eliminated, and /p.pcPXF can measure kinematics consistently.  \n<!-- image -->  \nto match the observed spectrum (Cappellari 2017), and the code returns one series of LOSVD moments which are best fitted to the observed spectrum. This means all the templates are assumed to have the same kinematics ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ) in a line-of-sight, i.e., stellar populations with different [ M / H ] and ages are assumed to have the same LOSVD. But in real galaxies, due to the joint process of chemical enrichment and dynamical movements of stars, populations with different [ M / H ] and ages are different in LOSVDs. This effect is non-negligible for edge-on projected galaxies because metallicity and age gradients along the disk are the strongest. Therefore, 𝑉 LOS should change with [ M / H ] and age, which disagrees with the analysis used in most of the previous studies when employing /p.pcPXF, even though /p.pcPXF allows such variations. Because of the relation between [ 𝛼 / Fe ] and [ M / H ] , stars with different [ 𝛼 / Fe ] should also have different 𝑉 LOS .  \nTo investigate it in our mock stellar catalog, we select one Voronoi bin and split the included particles into three groups using [ M / H ] and age, respectively. We reserve this investigation on [ 𝛼 / Fe ] for the future due to the limited number of [ 𝛼 / Fe ] templates. Then, we plot the LOSVDs of these groups of particles, weighted by each particle's total flux, and show them in the first row of Fig. 15. On the top left panel, particles with -0 . 8 < [ M / H ] < -0 . 4 dex have the sharpest distribution, while particles with -0 . 2 < [ M / H ] < 0 . 4 dex show a nearly normal distribution. In the top-right panel, particles with 12 < age < 14 Gyr have the broadest distribution, while particles with 0 < age < 8 Gyr have the narrowest and most peaked distribution. The red line and black line represent /p.pcPXF fitted and the true LOSVD for this Voronoi bin, with kinematic parameters written in the top right corner, respectively. Compared to the true  \nLOSVD,thewidth of /p.pcPXF fitted curve in red is close to the youngest populations (0 < age < 8 Gyr), which is reasonable because these populations dominate the light. However, the top right panel shows that /p.pcPXF tends to also care for the oldest populations by having a more skewed tail of their LOSVD at 𝑉 LOS around 150 -300 km s -1 . Even though the youngest populations dominate the spectral light, the old populations have the most features along the whole fitted wavelength region. Therefore, the differences of the skewness in the region of 𝑉 LOS at 150 -300 km s -1 indicate the effect caused by different spectral features having different LOSVDs. This figure strongly suggests there are limitations on current techniques to obtain unbiased kinematics due to the dependence of 𝑉 LOS on age and [ M / H ] .  \nWe made a new experiment to verify this effect on spectral fitting results. Firstly, we re-generate mock MUSE cubes following the same procedures as above, but not using each particle's 𝑉 LOS to shift the spectrum. Then we obtain non-kinematic mock data cubes. Next, for each Voronoi bin, we select all the particles included and obtain its 𝑉 LOS histogram, which is the true LOSVD as shown in grey in Fig. 15. Then, we directly use this histogram as the LOSVD kernel function and convolve it with all the spectra in this Voronoi bin to add kinematic effect. We do it for all the Voronoi bins and obtain new mock cubes, where spectral templates with different [ M / H ] and age have the same LOSVDs. Therefore, in this way, we eliminated the relation of 𝑉 LOS with [ M / H ] and age as mentioned above. Finally, we employ /p.pcPXF to the new cubes (hereafter LOSVD-Convolved cubes) and measure the kinematics. The fifth row of Fig. 14 shows the kinematics recovery results. It is clear that all the Gauss-Hermite moments are well recovered with all the data points aligned to the  \nzero line. The red and black lines in the bottom row of Fig. 15 also indicate the correction of /p.pcPXF measurements. We also plot the kinematics maps in the same way as Fig. 6 in Appendix (Fig. B4) to better demonstrate the consistency.  \nWe also tested the kinematics recovery on LOSVD-Convolved cubes generated using MILES templates, as shown in the sixth row of Fig. 14 (see also Fig. B5). In this low spectral resolution, even though we make different populations have the same LOSVD for the spectra, there are still systematic offsets appearing in all the panels. We also demonstrate in the last row the results on LOSVDConvolved cubes generated using PEGASE-HR templates, but in MUSE spectral resolution to remove the effects of different SSP templates (see kinematics maps in Fig. B6). The appearance of offsets indicates the effect purely due to insufficient spectral resolution.  \nIn conclusion, all these results strongly suggest that both the low spectral resolution and variance of LOSVDs in different age and [ M / H ] contribute to the systematic offsets we see in residual panels of Fig. 6. In the future, to remove the systematic offsets of stellar kinematics, one needs to use an instrument with higher spectral resolution than MUSE. This means more exposure time is needed to achieve the same SNR as MUSE observations. As we know 𝑉 LOS changes with [ M / H ] and age in the line-of-sight of real galaxy observations, it might be necessary to allow different templates to have different LOSVDs during the spectrum fitting analysis to obtain more accurate results. However, kinematic and stellar population measurements from spectral fitting methods have already been a highly degenerate problem. Allowing variable LOSVDs will increase the degeneracy by a factor of the number of templates. In this case, it is easy for the algorithm to obtain incorrect results. One possible way would be simply assuming a quantitative relation between LOSVD and metallicity and age at different locations of the galaxy and that this relation can be expressed by an analytical equation. This is equivalent to adding a prior to the fitting process to improve the accuracy of kinematics measurement without adding degeneracy. Even though this way has many limitations, it will be still useful for the analysis of Milky-Waylike galaxies, and G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc can help provide this prior for galaxies in different projections.  \nApart from the above analysis on inconsistencies caused by spectral fitting methods, Fig. 15 also shows the deviations between 𝑉 LOS histogram (grey line) and the parametrically best-fitted true LOSVD (black line). In the future, we will explore using non-parametric techniques (e.g., BAYES-LOSVD Falcón-Barroso & Martig 2021) to recover the true kinematic parameters rather than Gauss-Hermite equation, which might result in better fits when the LOSVD is strongly non-Gaussian.", '4.2 Inconsistency of mass fraction distributions recovery': "In general spectral fitting, measuring mass/light fraction distributions of stellar populations is a highly degenerate problem, especially when one wants to obtain the star formation history at different components. In this study, we follow the data analysis procedures in M21 to divide the template's age and metallicity bin size due to the definition of star formation history. However, most of the previous studies did not consider bin size of template grids. Nevertheless, the mass fractions increasing from young to old and metal-poor to metal-rich populations found in Fig. 12 is still seen in some components of the galaxies (e.g., Pinna et al. 2019a,b). In the mock and real data tests of Kacharov et al. (2018), the mass fraction distributions are also smoothed towards old and metal-rich regions (see their Fig. 5). The over-dense feature of young populations at around 2 -4 and 12 -14 Gyr in our Fig. 12 is also seen in some of the  \nabove-mentioned studies. In this section, we re-investigate all these tests more comprehensively and discuss the differences in the mass fraction distribution from /p.pcPXF, including whether or not to divide the bin size, effects caused by SNR and different SSP models, and try to address potential reasons causing the inconsistency to the true values. From now on, we refer 'mass fraction' as /p.pcPXF obtained mass weights of SSP templates normalized by the total weights (for each panel in the figures). And we refer 'mass fraction (Gyr -1 )' and 'mass fraction (dex -1 )' as mass fraction divided by the bin size of age and metallicity grids, respectively, where mass fraction (Gyr -1 ) is also equivalent to star formation history (SFH).  \nWefirst convert the age distribution x-axis of Fig. 12 from linear to log-scale, because MILES templates are closer to being evenly scaled in log age than linear age, then the width of each step is close to being equal. Next, we re-plot the age distribution recovery in Fig. 16 and calculate the correlation coefficients using the original age grids. The first row shows age distribution in mass fraction (Gyr -1 ), which is the SFH of each component. The thin disk, inner central, and global panels clearly show two overdensity regions, which are pointed out by two green boxes. The thick disk and upper central panels also have peaked features in these boxes, but not significant. The second row shows just the mass fraction without dividing the bin size, and all the panels show relatively smoother distributions than the first row. Therefore, given the differences between the first and second rows are whether or not dividing bin size, we speculate the two overdensity regions are artifacts due to the bin size arrangements of the MILES template grids.  \nTo verify our speculation, we plot in Fig. 17 a zoom-in of the two overdensity regions pointed out on the top right panel of Fig. 16. The blue and orange lines are mass fraction (Gyr -1 ) and mass fraction distributions, respectively, which are equivalent to the first and second row of Fig. 16, respectively. The mass fraction for each panel is multiplied by a factor for better visualization. We also plot the age bin size of MILES templates as the grey line and the scale in the right y-axis. From these two panels, we find the peaked features of mass fraction (Gyr -1 ) (blue line) appear when the age bin size experiences a decrease from old to young ages, and the mass fraction (orange line) is still relatively smooth. Therefore, we can confirm the peaked features of mass fraction (Gyr -1 ) in our results are due to the non-regularized spacing of MILES age grid. This also explains the better recovery of [ M / H ] distributions (bottom row of Fig. 12) which is due to the almost linearly spaced [ M / H ] grids.  \nPrevious studies applied regularization (McDermid et al. 2015) in /p.pcPXF to obtain smooth weight fraction distributions. The regularization works as an extra term in the equation to be minimized during the optimization to damp high-frequency variations in the weight distributions along spectral templates grid (see details in Section 3.5 of Cappellari 2017) and leads to smooth weights distributions. However, this smoothness has an effect on weights without dividing the bin size. This explains why mass fraction distributions are smooth in our results, and whether mass fraction (Gyr -1 ) distributions are smooth depends on how the template grids are spaced. For the rest of this section and Section 4.3, we only use the mass fraction without dividing the bin size to justify the recovery consistency (Fig. 18 and Fig. 20), because this is expected to be achieved by spectral fitting methods with regularization. And we leave investigations on mass fraction (Gyr -1 ) to Section 4.4 and 4.5.  \nAlthough the current /p.pcPXF regularization algorithm does not consider the bin size of the template age grid and different regularization strategies (e.g., regul values, order of regularization) can affect the mass fraction recovery (will be discussed in Section 4.3), we still see the inconsistency between spectral fitting results and true mass  \nFigure 16. Age distribution recovery of different components of the galaxy from /p.pcPXF on the mock MUSE cube generated in Section 3.1. This figure is similar to Fig. 12 but the x-axis is in logarithmic scale. Black lines are the true values and red lines are the results from /p.pcPXF. The first row shows mass fraction (Gyr -1 ), which is the representation of star formation history, and the two green boxes in the 'Global' panel are the two overdensity regions. The second row shows just the mass fraction, which is the direct output from /p.pcPXF. The correlation coefficients of the red and black lines are written on the top right corner for each panel.  \n<!-- image -->  \n<!-- image -->  \nFigure 17. Zoom-in of the two overdensity regions pointed out in Fig. 16. The blue and orange lines are mass fraction (Gyr -1 ) and mass fraction distributions, equivalent to the first and second row of Fig. 16, respectively. The mass fraction for each panel is multiplied by a factor for better visualization. We also plot MILES templates' age bin size as the grey line and the scale in the right y-axis.  \n<!-- image -->  \nfraction distributions of the oldest populations in the bottom row of Fig. 16. The main reason is that spectral templates in the oldest and metal-rich regions are very similar and indistinguishable at the current SNR level. Therefore, it is difficult for /p.pcPXF to correctly recover the mass distributions in these age-[ M / H ] regions when performing fitting using templates normalized to 1 M ⊙ , even though the algorithm has already reached the global minimum. The mass fraction uncertainties in these regions can be much larger (relative uncertainties from 1 𝜎 values dividing the mean of MC realizations are ∼ 20% for regions less than 8 Gyr and 50 -65% in 10 -14 Gyr). We also confirm in Appendix Fig. B7 that using LOSVD-Convolved cubes does not improve the results significantly, so the effect of systematic offsets in kinematics measurement in Fig. 6 is not the cause of mass fraction distribution inconsistency.  \nOne way to increase the quality of the mass fraction recovery in old and metal-rich population regions could be to increase the SNR of the spectrum. This can be achieved by requiring a higher target SNR during the process of Voronoi binning but with the cost of reducing the number of Voronoi bins, i.e., details in the spatial distributions. We show in Fig. 18 the age and [ M / H ] distributions at different SNRs (shown in different colors), compared to the true  \nvalues shown in black lines. This figure is obtained by running the GIST pipeline on the mock MUSE cube generated in Section 3.1 with different target SNRs during Voronoi binning, and then the mass fraction distributions of all the Voronoi bins are added together to represent the global mass distributions of the mock galaxy, and finally integrated through age and [ M / H ] axes, respectively. To make it a fair comparison, we make sure that the spatial pixels used to generate Voronoi bins are the same. According to the correlation coefficients, this figure indicates that with the increase of Voronoi bin SNR, the global mass fraction distributions in both panels are more consistent with the true values. Specifically, the weights of metalrich and old stellar populations are going lower with higher SNR. Therefore, Fig. 18 indicates that increasing SNR can improve SFH recovery. However, we note that in reality, an integrated spectrum with SNR larger than 1000 is impractical. In addition, even with SNR at 200 pix -1 is pointless because in this case numerical issues and spectral library uncertainty can become the dominant factor.  \nTo study the effect of different choice of SSP templates, we re-plot age and [ M / H ] distributions of mock cubes generated by PEGASEHR templates in Fig 19. PEGASE-HR templates are evenly spaced in log age. During the /p.pcPXF fitting, we use the same templates and  \nFigure 18. Global age (top panel) and [ M / H ] (bottom panel) distributions at different SNRs, which are obtained by running /p.pcPXF on the mock MUSE cube generated in Section 3.1 with different target SNRs of Voronoi bins, and then the mass fraction of all the bins are added together to represent the global mass distributions of the galaxy. Black lines are the true values. 'Corr' is the correlation coefficient value of the /p.pcPXF and true values at each SNR. This figure indicates that increasing SNR can help obtain better mass distributions compared to the true values, especially the issue of metal-rich, old stellar population overestimation is alleviated.  \n<!-- image -->  \nfitting strategies (fitting wavelength region, regularization, etc.) as those in Section 3.2 and also degrade to MUSE spectral resolution to remove the resolution differences (even though we find increasing spectral resolution does not improve the recovery significantly after comparing to Fig. B8 of a PEGASE-HR spectral resolution version). The Voronoi bin allocation is the same as mock cubes generated by MILES templates. According to the correlation coefficients, mass fraction using PEGASE-HR in Fig. 19 are more consistent than those using MILES in Fig. 12 and Fig. 16 for both age and [ M / H ] . Moreover, the mass fraction - log age panels (second row) have almost no overdensities than the results using MILES (second row of Fig. 16), and the oldest populations have a more modest increase. Therefore, in our test, PEGASE-HR performs better during the mass fraction recovery and could provide more consistent mass fraction distributions than MILES. We think this is likely due to the templates' perfectly even spacing in log age, which leads to smoother input and is preferred by /p.pcPXF during regularization. However, the PEGASE-HR templates used in our tests are generated using PADOVA 1994 (Bertelli et al. 1994) isochrones but MILES are generated using BASTI isochrones (Pietrinferni et al. 2004, 2006). Whether isochrones differences can affect mass fraction distribution recovery ability requires further investigations. Moreover, the current limitation of PEGASE-HR templates is they only have one bin  \nin [ 𝛼 / Fe ] , we emphasize here again the need for multiple-[ 𝛼 / Fe ] templates.  \nIn the future, non-parametrically, it is worthwhile to test the mass fraction recovery ability using templates with linearly spaced grids (see details in Section 4.4) or modifying regularization to apply smoothness effect in terms of mass weights (Gyr -1 dex -1 ). Another approach is to take the mean of Monte Carlo realization (e.g., Pessa et al. 2023), which is more consistent with the true values than regularization based on simple tests (see Fig. 1 of Cappellari 2023). We can verify whether this is true for a mock MW cube in future studies. Parametrically, one can try to recover both chemical enrichment history and star formation history by taking into account chemical evolution theories, which will help remove counter-intuitive features in the mass fraction distributions. One approach is to use Bayesian spectral fitting codes, such as Prospector (Leja et al. 2017) or BAGPIPES (Carnall et al. 2018). Another example is the semi-analytic spectral fitting method from Zhou et al. (2022). This method applies spectrum fitting with the predicted best-fit spectrum from chemical evolution models, which is similar to adding a prior on /p.pcPXF during the fitting. Fig. 6 of Zhou et al. (2022) showed that this method has better consistency on mass fraction recovery than /p.pcPXF, and the over-estimations in metal-rich, 2 -4 and 12 -14 Gyr regions are fixed successfully. Therefore, it provides a way to measure accurate chemical enrichment and SFH for future studies. One important caveat is that their chemical evolution model is a close-box model, which does not apply to galaxies dominated by frequent passive merger events (e.g. NGC 7793 in Kacharov et al. 2018). In addition, this method does not take into account the relation of chemical abundances with velocity, which are important indicators for radial migration and kinematic heating.", '4.3 Effect of regularization': "In this section, we focus on the effects of different regularization strategies on mass fraction distributions. There are two parameters controlling the regularization process, which are regul and reg\\_ord , where regul applies linear regularization to the weights during the /p.pcPXF fit and reg\\_ord controls the order of regularization. Kacharov et al. (2018) investigated mass fraction differences of mock data with different order of regularization (see their Fig. 5 and Fig. C2) and concluded that the choice of order of regularization does not affect the results over a significant level. Boecker et al. (2020) tested the mass fraction distributions recovery of the second- and third-order of regularization and compared them to the true values using stellar particles from an EAGLE simulated galaxy. They found results from third-order regularization are more consistent with the true values. Cappellari (2023) compared the mass fraction recovery consistency of non-regularized /p.pcPXF fit, an average of 100 /p.pcPXF fits to MC realizations, and single regularized /p.pcPXF fit (see their Fig. 1), and concluded that regularized results are comparable to that of the average of multiple realizations.  \nThe selection of regul value was briefly discussed in Section 3.2. Here we re-investigate this question in detail by applying different regul and reg\\_ord strategies on our mock cubes generated in Section 3.1. The results are shown in Fig. 20, where we plot the age distributions on the top panels and [ M / H ] distributions on the bottom panels for different components of the galaxy. We also plot the global mass distributions in the last column. We apply five different strategies during the fitting: the first three cases use firstorder regularization with regul = 5, regul = 10, regul = 20, and regul = regul 𝑚𝑎𝑥 ( 30 -100 ) , respectively; the last uses regul = 5 and second-order regularization. Black lines are the true values.  \n)  \n1  \nr  \ny  \nG  \n(  \nn  \no  \ni  \nt  \nc  \na  \nr  \nf  \ns  \ns  \na  \nm  \nFigure 19. Age and [ M / H ] distributions of each component for the mock cubes generated using PEGASE-HR templates and degraded to MUSE resolution (FWHM = 2 . 65 Å ). The first row demonstrates mass fraction (Gyr -1 ) changes with age (SFH). The second row shows the mass fraction changes with age (MDF). The third row shows mass fraction changes with log age. The fourth row shows cumulative mass fraction distributions and the fifth row shows metallicity distributions.  \n<!-- image -->  \nWhen comparing results with the same order of regularization but different regul , we find results with larger regul are smoother, and their weights at old populations are relatively smaller than results with smaller regul . However, in the [ M / H ] distribution panels, weights are going more toward metal-rich populations and the distributions become less consistent with the true values for all the components. For results with the same regul but different orders of regularization, we see results with reg\\_ord = 1 show better consistency with the true values than those with reg\\_ord = 2. Specifically, results with reg\\_ord = 2 have more mass fractions in the oldest populations.  \nIn conclusion, Fig. 20 indicates that for this mock cube, results using the first order of regularization are better than results using the second order of regularization; increasing regul can help obtain smoother results but the weights will go towards metal-rich regions. Given we do not have the codes to perform third-order regularization, it is hard to say which order is the best. Intuitively, it seems higherorder regularization tends to have more mass fraction in the old populations. Therefore, we think the third-order regularization is preferred for EAGLE simulated galaxy test (Boecker et al. 2020) because the particles' mass weights from the EAGLE simulations  \nare naturally dominated by old and metal-rich populations, then the third-order regularization tends to fit better in this region. Our EG/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc mocked catalog has no particle in the metal-rich and oldest regions, so the first-order regularization is better. However, all the above discussions need further verification and we will test the thirdorder regularization results in the future. For the real galaxy fitting results of Pinna et al. (2019a,b) and M21, they all used secondorder regularization and showed an overdensity of mass weights in the oldest, metal-rich populations for different components of the galaxies, which makes it difficult to discuss differences in the formation epoch between these components, and it is not possible to tell if the overdensity is true or instead from full-spectrum fitting artifacts.", '4.4 Mass fraction distribution recovery using interpolated MILES templates': "According to discussions in Section 4.2 and 4.3, the overdensity features of age distributions in mass fraction (Gyr -1 ) shown in Fig. 16 are likely due to not considering templates bin size during regulariza-  \nFigure 20. Age (top panel) and [ M / H ] (bottom panel) distributions for different components by running /p.pcPXF on the mock MUSE cube generated in Section 3.1 with different regularization strategies (shown in different colors). Black lines are the true values. The last column shows global age and [ M / H ] distributions. This figure indicates that results using the first order of regularization are better than results using the second order of regularization; increasing regul can help obtain smoother results but the mass fractions will go towards metal-rich regions for all the components.  \n<!-- image -->  \n1  \nFigure 21. Same as Fig. 12 but for the mock cube generated using MILES templates interpolated linearly in age and [ M / H ] in MUSE spectral resolution.  \n<!-- image -->  \ntion, and using different regularization strategies can barely alleviate the inconsistencies with the true values. Therefore, we investigate the improvements of age and [ M / H ] distributions using interpolated MILES templates in this section. Fig. 21 shows the age and [ M / H ] distributions of mock cubes generated by G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc and fitted by /p.pcPXF using MILES templates that are interpolated linearly in age and [ M / H ] grid. The fitting strategies are the same as those in Section 3.2. We can clearly see that the overdensities in age distributions (in Gyr -1 ) at 2 -4 Gyr are alleviated to a large extent in the 'thick disk', 'inner central', and 'global' panels compared to the results using original MILES templates in Fig. 12. In addition, [ M / H ] distributions (in dex -1 ) from /p.pcPXF are also more consistent to the true values in terms of correlation coefficients. The increase of old populations still exists, which will be discussed in detail in Section 4.5. We also investigate the results using MILES templates interpolated linearly in log age and [ M / H ] (shown in Fig. B9) and find the overdensities still exist with the shape being more similar to the results using PEGASE-HR (Fig. 19). Given PEGASE-HR template grid is also linearly spaced in log age, this similarity is expected. In conclu-  \nsion, we confirm that performing spectral fitting using linearly interpolated templates in age and [ M / H ] can alleviate the overdensities in 2 -4 Gyr. In addition, the new version of /p.pcPXF considering bin size during regularization is under development (private communication with Michele Cappellari). Given this is equivalent to interpolating templates as we did in Fig. 21, one can choose this approach in the future to avoid additional template uncertainties due to interpolation.", '4.5 Potential effects of effective age resolution on star formation history recovery': "In Section 4.2 to 4.4, we investigated the effects of template grid spacing and regularization on the recovery of mass fraction distribution. The comparison between Fig. 12 and Fig. 21 confirms that the over-estimation of age distribution (SFH) in the 2 -4 Gyr range is significantly alleviated by using templates with a linearly spaced age grid. However, in Fig. 21, which uses the linearly spaced age grid, narrow peaks at around 1 and 2 Gyr can still be seen prominently in the 'thin disk', and moderately in the 'thick disk', and 'inner  \ncentral'. These peaks differ from the broad 2 -4 Gyr overdensities seen in Fig. 12, but are more similar to peaks seen for the case of the log age grid (Fig. 19 and Fig. B9), although they suffer from the additional effect of inconsistent regularization. When using the template grid with linear age spacing, there is an over-estimation in 12 -14 Gyr range that is present across all components (Fig. 21). These results suggest that even after fixing the issue of inconsistent regularization with the linearly spaced age grid and trying out different SSP models, the recovery of SFH is far from perfect and a possible reason for this could be the fact that weights of some SSP templates are degenerate with respect to the observed spectra.  \nTo examine the effect of template weight degeneracies on SFH recovery, Zibetti et al. (2024) estimated and explored the age resolution of SSP templates that can be achieved through non-parametric fullspectrum fitting, which we refer to as effective age resolution. The effective age resolution is defined as the smallest width of independent age bins in a reconstructed SFH, and it is limited by the amount of information that is possible to extract from a spectrum. We refer to their Section 2 for a detailed analytical definition. Zibetti et al. (2024) found that the effective age resolution is generally no better than 40% for most age grids older than 100 Myr (see their Fig. 4). This fraction is roughly consistent with the mass weights uncertainties in our results. More importantly, Zibetti et al. (2024) found that the age resolution does not vary smoothly but shows some well defined maxima, e.g. at 10 Myr and 1 Gyr. The origin of this could be related to the fact that this is the stage where some unique spectral features develop due to significant changes in stellar evolution. They speculated (their Sec 5.5) that these maxima in age resolution could be responsible for the peaks seen in SFH of galaxies extracted using non-parametric full-spectrum fitting on IFS data (e.g., CALIFA Cid Fernandes et al. 2013; González Delgado et al. 2017 and MaNGA Sánchez et al. 2019). These peaks in SFH appear across a wide variety of galaxies and span of cosmic time, which lends further support to the idea. Our work where we try to recover the SFH of simulated galaxies further supports this idea- even though the SFH of the simulated galaxy is smooth, the recovered SFH has peaks/overdensities (e.g., see the age distributions in mass fraction (Gyr -1 ) shown in 16, 19, 21, B8, B9). The over-estimation of SFH at 12 -14 Gyr in our results could be due to a combination of a number of factors. Our input SFH as specified in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc drops sharply to zero at 12.8 Gyr which makes it challenging to recover the SFH. Unfortunately, the spectra of templates at old ages are quite similar to each other (low effective age resolution), this would lead the fitting machinery to assign non-zero weights for old templates. The regularisation which favours smoothness will also favour non-zero weights for old populations, even with the input truncating at 12.8 Gyr. Additionally, due to the low luminosity of these old SSPs, to account for a given flux in spectra, a significant mass fraction can be assigned to these old ages as compared to young ages. The linear age grid has more templates at old ages as compared to a log age grid. This suggests that bias in SFH (which is mass fraction per unit age) due to the above effects is expected to be higher for a grid with linear age as compared to log age. This is evident in our results and can be seen by comparing the SFH at old ages in Fig. 21 and Fig. B9.", '4.6 Caveats when using G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc': 'In this section, we give some caveats when using the G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code. The E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc mock stellar catalog is generated based on the GCE model of S21. The success of this model is that it can reproduce [ 𝛼 / Fe ] -[ Fe / H ] bimodality which is quantitatively consistent with the APOGEE observations of the MW (Hayden et al. 2015) by  \nemploying radial migration and kinematic heating. However, more complex kinematic and dynamical processes are likely to exist in real galaxies, such as non-asymmetric perturbations due to spiral arms, bars, interlopers, etc. The current model does not have these processes and features, which should be added in the future. In addition, the interpretation of [ 𝛼 / Fe ] bimodality in S21 due to radial migration and kinematic heating is opposite to findings from simulations which suggest the need of mergers (e.g., Clarke et al. 2019; Buck 2020), and this model lacks stars in the Galactic halo. Therefore, this code is not applicable for comparing the merger features or the halo distributions; This model also lacks descriptions of parameter distributions in the bulge, especially the nuclear disk in the MW center (e.g., Schultheis et al. 2021) and chemodynamics of the B/P bulge (e.g., Wylie et al. 2021), which could be improved in the future.', '5 FUTURE STUDIES': "In this study, we did not apply extinction on the mock data cube because this work is purely on testing the kinematics and stellar population recovery via full-spectrum fitting. In real edge-on galaxies, extinction can have an essential effect because it can obscure nearly 50% of the total light from stars in the thin disk, making SNR in this region lower. With the G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code, we can test the kinematics and stellar population recovery of mock cubes with extinction. The extinction model in E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc is assumed as a double exponential distribution mainly in the thin disk, which can be integrated through the line-of-sight, and each particle will have a 𝐸 ( 𝐵 -𝑉 ) value. Then, we can add the extinction effect to the spectrum by using a specific reddening curve (e.g., Cardelli et al. 1989; Calzetti et al. 2000), which is also the strategy that /p.pcPXF applies to estimate extinction. Ge et al. (2018) has already tested the effect of extinction on fullspectrum fitting results using one stellar population template, we can re-investigate this question with comprehensive mock IFS MilkyWay data from the Galactic chemical evolution model S21 and also study how different extinction models can affect the projected galaxy properties. In addition, several non-asymmetric features should be added to the chemical evolution model like the spiral arms, bars, mergers, and halo particles to mimic a more realistic MW. It is also necessary to add gas emissions spectral features.  \nThere are several improvements can be made to the choice of SSP models. The current version of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc has embedded MILES and PEGASE-HR templates. In the future, we will add other empirical models such as BC03 (Bruzual & Charlot 2003), X-Shooter (Verro et al. 2022), and theoretical models like STARBURST99 (Leitherer et al. 1999) and BPASS (Eldridge et al. 2017), which can satisfy different science goals since different models have different advantages. Some interesting tests can also be executed, e.g., using one SSP model to generate mock cubes and fit them with another model, which can help study the effect of SSP template uncertainties on full-spectrum fitting results.  \nIn addition, even though we showed in Section 3.4 that the [ 𝛼 / Fe ] -rich and [ 𝛼 / Fe ] -poor populations can be well identified using two [ 𝛼 / Fe ] bins, it is still essential to add more bins to obtain detailed [ 𝛼 / Fe ] distributions. The MW results in Hayden et al. (2015) have shown that stars in the thin and thick disk have [ 𝛼 / Fe ] values about 0 . 0 and 0 . 2 dex, respectively, and there are some stars with negative [ 𝛼 / Fe ] , which is also shown in the E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc catalog. However, it is challenging to use current MILES 𝛼 -variable templates to identify these two [ 𝛼 / Fe ] sequences because most of the weights are in the bin of [ 𝛼 / Fe ] = 0, and it is also difficult to assign a spectrum to particles with negative [ 𝛼 / Fe ] . Therefore, templates with more [ 𝛼 / Fe ] bins  \nsuch as sMILES (Knowles et al. 2023) and FSPS (Conroy et al. 2009; Conroy & Gunn 2010) will help on the study of [ 𝛼 / Fe ] -[ M / H ] distributions for external galaxies at different components.  \nAccording to Section 4.1, given that spectral resolution is important in stellar kinematics measurements, it is also necessary to increase both the instruments' and stellar population models' spectral resolution. Current IFS observations such as MUSE perform well on kinematics measurements of galaxies in the central region. However, to go deeper and obtain distributions of the outer regions, especially the outer thin disk, which normally has lower dispersion, one might need to use instruments and SSP templates with a higher spectral resolution to obtain more accurate kinematics maps.  \nAs discussed in Section 4.1, other than increasing spectral resolution, a prior can be added to the spectral fitting process which assumes the relation between LOSVD and metallicity and age at different locations of the galaxy. Then it can allow templates to have different LOSVDs and the degeneracy is not increased. The other way is to add the evolution of stellar kinematics to the semi-analytic model like the one in Zhou et al. (2022) and make it derive the chemical evolution histories along with kinematic processes directly by fitting with the integrated spectra, where the fitting process is constrained by several model parameters to be more physical. We will test the feasibility of these methods in the future.  \nAs for weight fraction distributions, Section 4.2 demonstrated how the bin size can artificially create peak features in the star formation history, and in Section 4.4 we verified that fitting the spectra using templates that are interpolated with grids linearly spaced in age and [ M / H ] can alleviate this effect. Given interpolating semi-empirical spectral libraries is not recommended when fitting with real observational spectra, our tests indicate the need for considering bin size during regularization. In the future, we will also test modifying the regularization algorithm by considering the bin size in age and metallicity. Parametrical methods like Zhou et al. (2022) can also be used to test recovery abilities. Moreover, in this work, we input a smooth SFH for the tests. Whether the mass fraction distribution inconsistencies we found also apply to merger-dominated galaxies with multiple starburst phases can be investigated using G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc generated mock cubes of N-body/hydrodynamical simulations.  \nThe aim of G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code is to apply the knowledge we learned from the MW to other MW-like galaxies by mocking the IFS observations as a bridge. Especially, we want to interpret physical processes such as radial migration, kinematic heating, and [ 𝛼 / Fe ] -[ M / H ] distributions of other galaxies to answer the question of whether all the MWAs have similar processes to the MW, and if not, how much difference can be. However, we have not yet performed any comparison of the mock MW IFS data cube with existing MWA observations due to the lack of extinction and some Galactic components in the current version of E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc, and the artifacts created by full-spectrum fitting methods we found in this work. Therefore, when all these aspects get improved in the future, we will compare our mock MW data cube from G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc with real IFS observations from the GECKOS survey (van de Sande et al. 2024), which will observe 35 edge-on MW-like galaxies, and see whether the MW is similar to the socalled MWAs according to their kinematics and stellar population properties.", '6 SUMMARY AND CONCLUSIONS': "In this work, we present the G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code which uses simplestellar population models and the mock MW catalog from EG/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc to generate a mock MW IFS data cube in the same format  \nas extragalactic observations. We aim to eliminate the differences in analysis techniques between the Galactic and extragalactic studies such as the use of individual stars vs. integrated stellar populations, the number density distributions of stellar parameters vs. integrated mass- or light-weighted distributions, the results with or without projection effect, etc. The mock data cube can be put into the GIST pipeline to perform data analysis in the same way as extragalactic observations and the results can be compared directly to external galaxies to study the similarities and differences between the MW and its analogs. Therefore, this code is a bridge to link the Galactic and extragalactic studies to understand the formation and evolution of disk galaxies.  \nThe G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code is flexible, allowing users to choose their preferred SSP templates, galaxy distance and inclination, spatial/spectral resolution, and field-of-view of the instrument. G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc is also designed to mock current existing instrumental observations such as MUSE, SAMI, or MaNGA, as well as test the performance of future instrument designs. Moreover, it can also be applied to Nbody/hydrodynamical simulations that contain MW-like galaxies and generate mock observations for comparisons.  \nFor the rest of the paper, we applied /p.pcPXF (included in the GIST pipeline) as a representation of full-spectrum fitting methods on G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc generated mock cubes to test the ability of /p.pcPXF to recover kinematics and stellar population properties. The mock MUSE data cubes have a distance of 26 . 5 Mpc and inclination 86 · , which is the same as MUSE observations on NGC 5746 from M21.  \nAfter comparing the true values calculated from E-G/a.pc/l.pc/a.pc/x.pc/i.pc/a.pc particles' properties with /p.pcPXF results, we found that there are systematic offsets between kinematic moments ( 𝑉 LOS , 𝜎 , ℎ 3 , ℎ 4 ) and the true values. We confirm two reasons are causing it:  \n- · The velocity dispersion 𝜎 of most Voronoi bins are smaller than the instrument velocity dispersion 𝜎 inst in MUSE spectral resolution (FWHM = 2 . 65 Å ).\n- · 𝑉 LOS changes with age and [ M / H ] for particles through the line-of-sight, but most previous studies assumed all stellar population templates have the same LOSVD during the spectral fitting.  \nBy using the higher spectral resolution templates PEGASE-HR (FWHM = 0 . 55 Å ) and applying the LOSVD-Convolved cubes which eliminate degeneracies of 𝑉 LOS with age and [ M / H ] , we can obtain consistent kinematic results with the true values. Therefore, we indicate the need to allow different templates to have different LOSVDs or assume a quantitative relation between LOSVD and metallicity and age at different locations of the galaxy, where the equation coefficients can be measured during spectral fitting. The latter method is equivalent to adding a prior to the fitting process to measure kinematics without adding degeneracy. We will perform tests in the future to verify if kinematics recovery can be improved. In addition, our tests also indicate the need to use non-parametric methods such as BAYES-LOSVD (Falcón-Barroso & Martig 2021) rather than Gauss-Hermite equation for more accurate measurements of 𝑉 LOS distributions.  \nIn terms of stellar population properties, we verified that spectral fitting methods can recover light- and mass-weighted age, [ M / H ] , and [ 𝛼 / Fe ] with good consistency, except the mass-weighed age and [ M / H ] are overestimated in the outer and inner regions of the galaxy, respectively. Both the [ 𝛼 / Fe ] -rich and [ 𝛼 / Fe ] -poor disk structures can be identified with reasonable regularization values during the fitting. Line-strength indices method can also identify these structures. Wefoundthe weight fraction distributions of stellar populations from spectral fitting using MILES templates on MILES-generated mock cubes are mostly consistent with the true values. However, some devi-  \ne shown in mass fraction distributions, where mass fractions normalized by the bin size are overestimated in regions of 2 -4 Gyr, 12 -14 Gyr, the most metal-rich regions, and underestimated in regions of 5 -10 Gyr. This is be due to many reasons including having low SNR that the flux error is larger than templates' similarities in the region of old and metal-rich populations, the uneven spacing of age grids in MILES templates, not considering bin size in the current regularization algorithm. We verified that having a regularly spaced age grid (linear or log) resolved the sharp discontinuities originally observed in the SFH with the default MILES library. The log age grid is still incompatible with the way the regularization is implemented in /p.pcPXF, so one should either alter the regularization scheme or use linear age spacing. Even when using the linear age grid, we see narrow peaks at around 1 -2 Gyr in the SFH. This is most likely due to the fact that the precision of estimating the age of an SSP varies significantly with age as pointed out by Zibetti et al. (2024). In addition, we found that the first-order regularization can obtain better mass fraction distributions than the second order regularization. Our tests and conclusions are helpful in identifying limitations of extragalactic data analysis using current methods and provide a reference for potential improvements in the future.  \nEven though G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc provides a bridge to link the Galactic and extragalactic studies by transferring MW to mock IFS observations, there remains a need for future improvements to facilitate more accurate measurements of external galaxy properties such as [ 𝛼 / Fe ] -bimodality and to enable detailed comparisons of the MW with MW-like galaxies. These improvements include more accurate SSP models with higher spectral resolution and more [ 𝛼 / Fe ] grids, more advanced spectral fitting algorithms, and instruments with deeper observations. Future tests using G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc are required including modifying spectral fitting codes to improve the recovery ability of kinematics and stellar population properties, using templates with more [ 𝛼 / Fe ] grids, and adding extinction, to achieve measurements of known parameters such as radial migration and kinematic heating efficiencies of the MW in external galaxies.", 'ACKNOWLEDGEMENTS': "Wethank James Binney, Alina Boecker, Martin Bureau, Andy Casey, Michele Cappellari, Scott Croom, Eric Emsellem, Jesus FalconBarroso, Jianhui Lian, Richard McDermid, Martin Roth, Anil Seth, Juntai Shen, Yuan-Sen Ting, Glenn Van de Ven, Gail Zasowski, Ling Zhu and Zhuyun Zhuang for useful discussions. We thank all the attendees who joined the discussions of this work during the conference Linking the Galactic and Extragalactic (Wollongong, Australia) and GECKOS-Oxford conference (Oxford, UK). We also thank the referee for providing useful comments.  \nZWacknowledgestheHPCserviceatTheUniversityof Sydney for providing HPC resources that have contributed to the research results reported in this paper. Works in this paper were also done using Yoga 2 , a privately built Linux server for astronomical computing.  \nZW is supported by Australian Research Council Centre of Excellence for All Sky Astrophysics in Three Dimensions (ASTRO3D) through project number CE170100013. MRH acknowledges support from ARC DP grant DP160103747 and ASTRO-3D. SS is funded by ASTRO-3D Research Fellowship and JBH's Laureate Fellowship from the Australian Research Council. JBH is supported by an ARC Australian Laureate Fellowship (FL140100278)  \nand ASTRO-3D. JvdS acknowledges the support of an Australian Research Council Discovery Early Career Research Award (project number DE200100461) funded by the Australian Government. FP acknowledges support from the Spanish Ministry of Science and Innovation (MICINN) through the project PID2021-128131NB-I00 / 10.13039/501100011033.  \nThis research has also made use of Astropy 3 , a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007) and SpectRes (Carnall 2017).", 'DATA AVAILABILITY': 'G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc source code is publicly available via https://github. com/purmortal/galcraft . The test data and figures can be shared with reasonable requests.', 'REFERENCES': 'Spitoni E., Silva Aguirre V., Matteucci F., Calura F., Grisoni V., 2019, A&A, 623, A60', 'APPENDIX A: MODIFICATIONS OF THE GIST PIPELINE': 'We provide a list of modifications on the GIST pipeline (Bittner et al. 2019) to improve the flexibility of this software: We add more templates such as the original and interpolated PEGASE-HR (Le Borgne et al. 2004) to the software and allow the option to oversample the spectra when degrading to lower spectral resolution, which is for the same reason as we do for the G/a.pc/l.pcC/r.pc/a.pc/f.pc/t.pc code in Section 2.3; We also added the option to select SSP templates with a certain age and [ M / H ] range in some special cases. When measuring stellar kinematics and stellar population properties, we add the options to choose velscale\\_ratio , reg\\_ord , and use real spectra noise during the fitting and normalizing the integrated spectrum and noise by the median of the spectrum, which will be important to perform the iteration for estimating regul 𝑚𝑎𝑥 ; For regularization, we add an option to estimate regul 𝑚𝑎𝑥 following the procedures in McDermid et al. (2015) and save the stellar population results with regul = regul 𝑚𝑎𝑥 . We also allow further degrading to lower spectral resolution during the fitting when measuring SFH, which is mainly for the reduction procedures in M21 where they wanted to compare the mass-weighted stellar population parameters from /p.pcPXF with the SPP-equivalent parameters from line-strength indices in the same spectral resolution; We add the option to change the penalization value bias during /p.pcPXF fitting and measure templates weights uncertainties using Monte Carlo realizations similar to procedures in Kacharov et al. (2018), M21 and Cappellari (2023).  \nFigure B1. Same as Fig. 6 but turning off the kinematics penalization during /p.pcPXF fitting.  \n<!-- image -->  \nFigure B2. Same as Fig. 6 but for the mock cube generated and fitted by PEGASE-HR templates in PEGASE-HR spectral resolution (FWHM = 0 . 55 Å ).  \n<!-- image -->', 'APPENDIX B: COMPLEMENTARY FIGURES': 'This paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure B3. Same as Fig. B2 but turning off the kinematics penalization during /p.pcPXF fitting.  \n<!-- image -->  \nFigure B4. Same as Fig. 6 but for LOSVD-Convolved mock cubes generated and fitted by PEGASE-HR templates in PEGASE-HR spectral resolution (FWHM = 0 . 55 Å ).  \n<!-- image -->  \nFigure B5. Same as Fig. 6 but for LOSVD-Convolved mock cubes generated and fitted using MILES templates in MUSE spectral resolution (FWHM = 2 . 65 Å ).  \n<!-- image -->  \nFigureB6. SameasFig.6butforLOSVD-ConvolvedmockcubesgeneratedandfittedbyPEGASE-HRtemplatesinMUSEspectralresolution(FWHM = 2 . 65 Å ).  \n<!-- image -->  \nFigure B7. SameasFig.12andFig.16butforLOSVD-ConvolvedmockcubesgeneratedusingMILEStemplatesinMUSEspectralresolution(FWHM = 2 . 65 Å ).  \n<!-- image -->  \n1  \nFigure B8. Same as Fig. 19 but for mock cubes generated using PEGASE-HR templates in PEGASE-HR spectral resolution (FWHM = 0 . 55 Å ).  \n<!-- image -->  \n1  \nFigure B9. Same as Fig. 12 but for mock cubes generated using MILES templates interpolated linearly in log age and [ M / H ] in MUSE spectral resolution.  \n<!-- image -->'}