Chemical abundances of Young Massive Clusters in NGC 1313

We analyze spectroscopic observations of five young massive clusters (YMCs) in the barred spiral galaxy NGC 1313 to obtain detailed abundances from their integrated light. Our sample of YMCs was observed with the X-Shooter spectrograph on the Very Large Telescope (VLT). We make use of theoretical isochrones to generate synthetic integrated-light spectra, iterating on the individual elemental abundances until converging on the best fit to the observations. We measure abundance ratios for [Ca/Fe], [Ti/Fe], [Mg/Fe], [Cr/Fe], and [Ni/Fe]. We estimate an Fe abundance gradient of $-$0.124 $\pm$ 0.034 dex kpc$^{-1}$, and a slightly shallower $\alpha$ gradient of $-$0.093 $\pm$ 0.009 dex kpc$^{-1}$. This is in contrast to previous metallicity studies that focused on the gas-phase abundances, which have found NGC 1313 to be the highest-mass barred galaxy known not to have a radial abundance gradient. We propose that the gradient discrepancy between the different studies originates from the metallicity calibrations used to study the gas-phase abundances. We also observe an age-metallicity trend which supports a scenario of constant star formation throughout the galaxy, with a possible burst in star formation in the south-west region where YMC NGC 1313-379 is located.


INTRODUCTION
Understanding how galaxies evolve chemically continues to be a fundamentally important topic in modern astrophysics. Stellar abundances are emerging as powerful tools for tracing the chemical composition of Galactic and extragalactic stellar populations, as the signatures of the gas reservoirs that formed such populations are preserved in their atmospheres. Through the analysis of abundance patterns of individual elements, particularly those which are sensitive to the different timescales of galactic evolution, one can reliably trace the enrichment history of individual galaxies. For example, the ratio of α elements (i.e. Ca, Ti, Mg, Si) to Fe-peak elements (i.e. Fe, Mn, Cr, Sc) can be used as indicator of the relative contribution from core-collapsed supernovae (SNe) and type Ia SNe (Tinsley 1979;McWilliam 1997).
In the last decades, chemical abundance studies of individual stars in the Milky Way (MW) have provided us with an extremely detailed picture of various Galactic components. The analysis of high-dispersion spectroscopy of individual stars in our Galaxy have shown that halo and bulge stars are primarily old with enhanced abundance ratios of α over Fe (Sneden et al. 1979;Fulbright 2000). Stars in the disk, on the other hand, are observed to be young with abundance trends similar to those observed in the sun (Edvardsson et al. 1993;Bensby et al. 2003;Reddy et al. 2003;Bensby et al. 2014). A natural next step is to perform similarly detailed studies on individual stars outside of the MW, and beyond the Local Group.
In star-forming galaxies, most of the abundance work is done through the analysis of emission lines from H II regions (Searle 1971). The gas-phase abundances obtained from these objects, however, do not usually pro-vide constraints on ratios such as [α/Fe]. An alternative approach to characterize the chemical patterns in extragalactic environments is the spectroscopic technique developed in the last couple decades, relying on both red (RSG, Davies et al. 2010) and blue supergiants (BSG, Bresolin et al. 2006Evans et al. 2007). Such studies have shown that accurate metallicities can be obtained through the analysis of spectroscopic observations from RSGs/BSGs in agreement (∼0.1-0.3 dex) with those inferred from H II regions (Bresolin et al. 2009;Davies et al. 2017).
Another promising probe of detailed abundances in extragalactic environments are star clusters. Since the spectra of the integrated light (IL) of stellar clusters are broadened by a few km s −1 , the individual stellar absorption lines are easily resolved. High-resolution (R ∼ 30,000) spectra of globular clusters (GCs) have been reliably used to study the chemical histories of the MW and nearby galaxies (Bernstein & McWilliam 2002;McWilliam & Bernstein 2008;Colucci et al. 2009Colucci et al. , 2011Colucci et al. , 2012Colucci et al. , 2017Larsen et al. 2012Larsen et al. , 2017Larsen et al. , 2018a. More recently the analysis of the IL of star clusters was expanded not only to younger stellar populations than GCs (Larsen et al. 2006(Larsen et al. , 2008 but also to lower resolution observations (R 8,000, Gazak et al. 2014;Lardo et al. 2015;Hernandez et al. 2017Hernandez et al. , 2018a. The work presented here aims at exploiting these recently developed IL analysis tools to dissect the chemical history of the nearby barred spiral galaxy, NGC 1313. At a distance of 4.4 Mpc (Jacobs et al. 2009), NGC 1313 has caught the attention of many studies due to its peculiar morphology (Figure 1), possibly indicating an interaction with a satellite companion (Sandage & Brucato 1979;Blackman 1981). Larsen et al. (2007) and Silva-Villa & Larsen (2012) have concluded that there is a particular increase in star formation activity in the south-west region of this galaxy supporting tidal interaction scenarios previously proposed by Peters et al. (1994) and Blackman (1981). Interestingly, one of the most massive star clusters, NGC1313-379, is located in this precise region and was previously studied by Hernandez et al. (2017). The detailed abundances of this YMC showed close-to-solar [α/Fe] = +0.06 ± 0.11 dex, also similar to the nebular metallicities measured for two nearby H II regions (Walsh & Roy 1997). We list in Table 1 some of the properties of this nearby galaxy.
Our analysis aims at expanding the exploratory work of Hernandez et al. (2017) by adding five more YMCs distributed throughout NGC 1313. Using the IL analysis tools first developed by Larsen et al. (2012) along with observations taken with the X-Shooter spectrograph on the Very Large Telescope (VLT), we perform  (Ryder et al. 1995) a detailed abundance study of the young stellar populations in this nearby barred spiral galaxy. The contents of this paper are structured as follows: Section 2 describes the spectroscopic observations and the data reduction, Section 3 details the analysis done, in Section 4 we present our results, and Sections 5 and 6 present our findings and our conclusions.

OBSERVATIONS AND DATA REDUCTION
The analysis presented here relies on X-Shooter spectroscopic observations taken as part of the VLT guaranteed time observation (GTO) program ID: 084.B-0468A. The data were collected on November 2009, executing in standard nodding mode with an ABBA sequence. The X-Shooter spectrograph provides observations spanning a wavelength coverage ranging between 3000 and 24800 A at resolving power of R =3000-17000, depending on the configuration used. This instrument makes use of a three-arm system to simultaneously observe in three bands, UVB (3000-5600Å), VIS (5500-10200Å), and NIR (10200-24800Å). Our observations use slit widths of 1.0 , 0.9 , and 0.9 obtaining resolutions of R ∼ 5100, 8800, and 5100 for the UVB, VIS, and NIR arms, respectively.
The YMCs in our sample were primarily chosen from a cluster compilation by Larsen (2004). The selection sample was focused on isolated targets. We highlight the location of the YMCs studied here in Figure 1, also including the cluster studied in Hernandez et al. (2017. We list in Table 2 the physical properties of the clusters in our sample. In Table 3 we present information on the X-Shooter observations including exposure times used for the different arms, exposure dates, and signal-to-noise (S/N) ratios for the individual arms. We point out that given the low S/N ratios observed in the NIR data we based the analysis presented here solely on the UVB and VIS observations. The observations were calibrated using the X-Shooter pipeline (v.2.5.2) with the European Southern Observatory (ESO) Recipe Execution Tool (ESOREX) v.2.11.1. The software allowed us to perform the basic data reduction steps including dark and bias corrections, flatfielding, wavelength calibration and sky subtraction. Similar to the approach in Hernandez et al. (2017), we make use of the IDL routines by Chen et al. (2014) to extract the science spectrum from the calibrated 2D images and combine the different orders using a varianceweighted average in the overlapping regions. We flux calibrate the extracted 1D spectra using the spectrophotometric standard Feige 110. Furthermore, we create response curves for each of the Feige 110 exposures using the same bias and master flats as those used in the reduction of the science observations, correcting for exposure time and atmospheric extinction. As highlighted in Hernandez et al. (2018a), this last step of applying the same master flats in both the science and the response curves is critical to remove any contemporaneous flatfield features.
Telluric contamination from the Earth's atmosphere strongly affect the X-Shooter observations particularly the data from the VIS and NIR arm. We used the routines created by Chen et al. (2014) to remove any contribution from the telluric absorption features. The routines by Chen et al. (2014) are based on Principal Component Analysis (PCA) which rely on a telluric library of 152 spectra to remove and reconstruct the strongest telluric absorption. We highlight that although we apply the PCA technique by Chen et al. (2014) to correct for this contaminating features, throughout our analysis we aim at avoiding these strongly affected regions: 5800-6100Å, 6800-7400Å, 7550-7750Å, 7800-8500Å, 8850-10000Å.
We show in Figure 2 the individual spectra for each of our targets. We note that in all of our clusters we observe strong dichroic features around 5700Å. According to Chen et al. (2014) these features typically appear in the 1D spectra in slightly different locations, which makes removing them a complicated task. Similar to those wavelength regions affected by telluric contamination, we also exclude from our abundance analysis those wavelength ranges affected by this instrumental feature. Fully calibrated X-Shooter spectra for a sample of YMCs in NGC 1313. Our analysis focuses on the UVB and VIS X-Shooter arms, spanning a wavelength coverage of 3000 to 10,000Å. Note the strong dichroic feature around 5700Å observed in the VIS arm. We show in grey the regions contaminated by telluric absorption and dichroic features. These regions have been excluded from the fitting procedure.  Hernandez et al. (2017) b Taken from the LEGUS cluster catalog (Adamo et al. 2017) The work presented here makes use of the IL analysis tool ISPy3 (Larsen 2020). ISPy3 relies heavily on spectral synthesis and full spectral fitting. This technique was primarily developed to analyze high-resolution spectroscopic observations (R ∼ 40,000). However, it was recently shown that the same software can be used to obtain detailed abundances from observations at lower resolutions (R < 10, 000; Hernandez et al. 2017Hernandez et al. , 2018a using predefined wavelength windows which minimize line blending. Briefly described, our analysis involve an iterative process where the abundances are modified each iteration until the best synthetic model is obtained for a given spectroscopic observation. This is done executing four main steps: 1) Derivation of stellar parameters for each of the stars in the clusters, 2) Creation of stellar atmospheric models, 3) Creation of IL synthetic spectrum, 4) Comparison between the modeled spectrum and the spectroscopic observations. In the following sections we expand on each of these steps.

Stellar Parameters
The ISPy3 suite of Python routines rely on the information provided by the Hertzsprung-Russell diagram (HRD) of the star clusters under analysis. Three main approaches have been adopted for deriving the physical parameters from the HRD: [a] Color Magnitude Diagrams (CMD; Larsen et al. 2017), [b] CMD+Isochrones (Larsen et al. , 2018bHernandez et al. 2017), [c] Isochrones (Hernandez et al. 2018a(Hernandez et al. ,b, 2019. Method [a] is optimally applied for nearby or local stellar clusters where resolved CMDs are available covering all evolutionary stages present. Method [b] is commonly used for extragalactic stellar clusters where the resolved CMDs typically include only the brightest stars in the cluster. The HRD is then complemented with theoretical isochrones covering fainter magnitudes. In the analy-sis presented here, we generate HRDs using solar-scaled stellar isochrones from PADOVA v.1.2.S (Bressan et al. 2012), adopting the Isochrone only approach, [c]. The main driver for adopting this approach is the absence of CMDs for the YMCs in our sample. We note, however, that Hernandez et al. (2017) found that the differences between methods [b] and [c] have only a minor effect on the inferred abundances of YMC NGC 1313-379, of the order of 0.1 dex.
For each of the YMCs, we select the initial isochrones assuming a Large Magellanic Cloud (LMC)-like metallicity (Walsh & Roy 1997), [Z] = −0.3 dex, and the ages listed in Table 2 published in the Legacy Extragalactic UV Survey (LEGUS) cluster catalog (Adamo et al. 2017). From the corresponding isochrones we extract the stellar parameters, e.g. effective temperatures (T eff ), stellar masses (M), and surface gravities (log g), following a Salpeter (1955) initial mass function (IMF) and a lower mass limit of 0.4 M . We note that Hernandez et al. (2019) found that changing the choice of IMF (e.g. from Salpeter to Kroupa IMF) modifies the inferred metallicities on average by 0.1 dex. Lastly, similar to previous IL studies of YMCs we adopt microturbulent velocity values (v t ) depending on the effective temperature of the stars: v t = 2 km s −1 for stars with T eff < 6000 K, v t = 4 km s −1 for stars with 6000 < T eff < 22,000 K, and v t = 8 km s −1 for stars with T eff > 22,000 K.

Stellar Atmospheres
Using the parameters described above, we generate atmospheric models for each stellar type. Similar to the work by Hernandez et al. (2018a,b), we make use of the local thermodynamic equilibrium (LTE) onedimensional (1D) plane-parallel atmosphere modeling code ATLAS9 (Kurucz 1970) for any stars with T eff > 5000 K. For stars with T eff < 5000 K we then use MARCS atmospheric models instead, as these have been primarily developed with a special focus on late-type stars (Gustafsson et al. 2008). These MARCS models are LTE 1D plane-parallel or spherical. One important distinction between the two modeling codes is that we generate the ATLAS9 atmospheres on-the-fly, whereas the MARCS models have been precomputed and downloaded from their official website 1 . At this point in the analysis we set the initial abundances of the stellar population. We note that in the work presented here we adopt the Solar composition from Grevesse & Sauval (1998).  Hernandez et al. (2017) We also highlight that the entirety of our analysis is based on LTE modelling. We currently do not apply any non-LTE (NLTE) corrections as these are dependent on the individual stellar types, and our work focuses on the integrated light of the star clusters. Previous studies have looked into the effects of NLTE treatment in the analysis of IL spectroscopic observations of GCs and found that the overall metallicities can change by ∼0.05 dex (Young & Short 2019). A different study by Conroy et al. (2018), also focused on older stellar populations (> 1 Gyr), concluded that the assumption of LTE is adequate when studying Na and Ca. Lastly, Eitner et al. (2019) find that in many cases these NLTE corrections, specifically for Ba, Mg, and Mn, are significant (> 0.1 dex) for metal-poor GCs. Most of these IL NLTE studies have so far focused on old populations of stars, where our focus here is on YMCs.

Synthetic Stellar Spectra
As it has been extensively described in our previous studies, our analysis requires the creation of individual synthetic spectra for each stellar type in the population under study. We continue to use the line lists by Castelli & Hubrig (2004). For those atmospheres generated with ALTAS9 we use the SYNTHE software (Kurucz & Furenlid 1979;Kurucz & Avrett 1981), and TURBOSPECTRUM (Plez 2012) for the MARCS atmospheric models. We then coadd the individual spectra to create a single IL model spectrum.
The IL model produced by the spectral synthesis codes is created at resolutions of R ∼ 500,000. To match the resolution of the X-Shooter data, R ∼ 5000-9000, we fit for the best Gaussian dispersion, σ sm , and convolve it with our IL synthetic spectrum. This smoothing parameter, σ sm , accounts for the instrumental resolution, σ inst , and the line-of-sight velocity dispersion of the star cluster, σ 1D .

Synthetic Spectra vs X-Shooter Spectra
Once the resolution of the synthetic spectrum and that of the observations are similar, the ISPy3 tool matches the continuum of the X-shooter spectrum with that of the model using a cubic spline or polynomial, depending on the size of the wavelength window being analyzed. Our software allows us to define specific weights for a given pixel in the X-Shooter spectrum, anything between 0.0 to 1.0. To avoid contamination from emission lines or other problematic features in the IL observations we have set the weights of these wavelength regions to 0.0, to exclude them from our analysis.

Modelling uncertainties and systematics
In previous publications we have assessed the various uncertainties and systematics introduced by the different choices made in the modelling. As mentioned in Section 3.1, Hernandez et al. (2017) showed that the differences in the inferred abundances applying method [b] and [c] are very minor, of the order of 0.1 dex. Additionally, studying the integrated light of a large sample of stellar clusters Hernandez et al. (2018a) found overall differences <0.1 dex when comparing abundances inferred using solar-scaled (PARSEC) and α-enhanced isochrones (Dartmouth, Dotter et al. 2007). Regarding the uncertainties introduced by the different IMF selection, Hernandez et al. (2019) calculate abundances adopting a Kroupa IMF (Kroupa 2001) and a Salpeter IMF (Salpeter 1955) and found that the uncertainties introduced by this selection is on average of the order of <<0.1 dex. Overall, it has been found that the systematic uncertainties should be of the order of those quoted in Table 5 or smaller.

RESULTS
Before initiating our abundance analysis described in Sections 3.1-3.4, we first obtain an estimate of the radial velocities (v rv ) of the individual YMCs. We list in the second column of Table 4 our inferred v rv and their uncertainties. We note that the inferred radial velocities for the different YMCs agree with the velocity maps for NGC 1313 by Ryder et al. (1995).
As part of the abundance analysis, the ISPy3 software repeats the steps described in Section 3, modifying the chemical abundances in each iteration. We fit for the individual abundances, converging on the best value by minimizing χ 2 and applying a golden-section search technique.
As a first step in our analysis we fit for the overall metallicity, [Z], and the smoothing parameter, σ sm , simultaneously. This is done by scanning the X-Shooter UVB and VIS wavelength ranges, 200Å at a time. For each YMC we select the initial isochrone adopting the ages in Table 2 and a metallicity of [Z] = −0.3 dex (Walsh & Roy 1997). We find that in general using the stellar parameters from the initial isochrone we successfully converge to an overall metallicity of [Z] ∼ −0.3 (±0.1) dex for each of the clusters. These overall metallicities are solely used as a scaling factor in the following stages where we measure the individual element abundances.
To measure the abundances of the elements of interest we adopt the predefined wavelength windows used in Hernandez et al. (2017), as these have shown to help minimize the degree of blending. We begin by fitting those elements with the largest number of lines in the X-Shooter UVB and VIS wavelength ranges, one at a time and setting the smoothing parameter fixed. We first measure Fe, followed by Ti, Ca, Mg, Ni, and Cr. We attempted to measure the abundances of Sc and Mn, however, the ISPy3 software was not able to converge on a final value. We list the individual elements, wavelength bins, best-fit abundances and their corresponding 1σ uncertainties calculated from the χ 2 fit, in Tables 6 through 10 in the Appendix section. In our previous studies we have found that the standard deviation, σ STD , is more representative of the actual measurement uncertainties, given that the scatter in the individual measurements is larger than the formal errors based on the χ 2 analysis. Similar to the work of Hernandez et al. (2017), we convert the measured σ STD into errors on the mean abundances adopting their equation (5). In Table 5 we list the weighted average abundances and their corresponding uncertainties (σ STD ). Lastly, we show example synthesis fits for each of the targets studied here in Appendix Figures 9-18.

Balmer emission
We note that three of the YMCs analyzed here displayed strong emission lines, particularly in Hα. The three different clusters, NGC 1313-463, NGC 1313-503, and NGC 1313-505, exhibit slightly different emission Cluster vrv σ1D (km s −1 ) (km s −1 ) NGC1313-201 403 ± 2 6.8 ± 3.5 NGC1313-439 483 ± 4 5.6 ± 3.8 NGC1313-463 484 ± 2 5.5 ± 5.3 NGC1313-503 510 ± 7 5.0 ± 3.4 NGC1313-505 469 ± 2 6.4 ± 1.5 profiles ( Figure 3). The Hα profile observed in YMC NGC 1313-503 is comparable to those identified in Hernandez et al. (2017) for YMCs NGC 1313-379 and NGC 1705-1. The nature of the broad profiles, similar to those observed in these three clusters, is briefly discussed by Melnick et al. (1985), where they proposed that one possibility for the origin of this emission in Hα, other than ionized gas, is the presence of Be stars. These objects are identified as B-type non-supergiant stars with rotational velocities of several hundred km s −1 which display strong and broad Hα emission features (Marlborough 1982;Townsend et al. 2004). Additionally, studies have found an enhancement of Be stars in stellar clusters with ages <100 Myr (e.g., Mathew et al. 2008). The broad Hα emission observed in NGC 1313-503, reaching velocities of ∼300 km s −1 , along with its estimated age of 20 Myr, hint at the presence of Be stars.
The Hα profile observed in NGC 1313-463, a cluster with an estimated age of 40 Myr, appears to be slightly more complex. From visual inspection we believe the emission could be originating from both a population of Be stars (creating a slightly broad component, reaching velocities of 200 km s −1 ), and from ionized gas surrounding the cluster (much narrower component). Lastly, the emission features observed in NGC 1313-505 appear to be originating primarily from gas, although a population of Be stars cannot be fully discarded. We highlight that any wavelength regions contaminated by emission features, irrespective of their origin, are masked in our abundance analysis.

DISCUSSION
In this section we discuss our abundance measurements from both α-and Fe-peak elements, and their implications on the enrichment history of NGC 1313. We also discuss in detail the abundance trends observed in this nearby galaxy and compare them to those studied in different environments, such as the MW and the LMC.

α-elements
Known as α elements, Ca, Ti and Mg, are primarily produced by high-mass stars and ejected to the interstellar medium through core-collapse supernovae (Woosley & Weaver 1995). Due to their nature, α-elements are typically used to gain insight into the period of time when core-collapse supernovae shaped the chemical evolution of the host galaxy. Since Fe-peak elements are also produced in massive stars, the result is a constant [α/Fe] ratio during those epochs dominated by corecollapse supernovae. Typically, the [α/Fe] ratios begin to decline when type Ia supernovae (SNIa) start to drive the chemical enrichment, as these release large amounts of Fe-peak elements (Timmes et al. 2003). Overall, these are the general trends we observed in the MW. Additionally, Galactic abundances of α elements appear to correlate with each other. Such a trend can be clearly observed in Figure 4, where we show in blue stars the abundances of field stars in the MW as measured by Reddy et al. (2003Reddy et al. ( , 2006. We also show in Figure 4 as orange stars the observed abundances in the LMC, a high-mass dwarf galaxy, published by Van der Swaelmen et al. (2013). We note that the LMC and NGC 1313 share a few properties. Both galaxies are classified as spiral galaxies, the LMC as SBm and NGC 1313 as SB(s) (Sandage & Tammann 1981), with comparable metallicities (Walsh & Roy 1997;Russell & Dopita 1992). The α-element abundances presented here, are the firstever measured in stellar populations in NGC 1313. The abundances for all the clusters in our sample, including those from Hernandez et al. (2017), are plotted as orange squares in Figure 4 as a function of [Fe/H]. The [Ca/Fe] ratios for all of the clusters in NGC 1313 tend to decrease with increasing [Fe/H], following a trend comparable to that observed in the LMC, and with slightly more depleted Ca abundances than those from the MW field stars. The [Ti/Fe] values observed in the YMC sample show a much higher scatter than that of Ca, however, a similarly high scatter is seen in the field stars of the LMC. This larger scatter in the LMC [Ti/Fe] ratios is primarily introduced by the stars in the inner disk, compared to the abundances in the LMC bar ( Van der Swaelmen et al. 2013). It is also clear from Figure 4 that overall, the Ti abundances are lower in NGC 1313, compared to those in the LMC for a given metallicity. We observe a comparable scatter in the [Mg/Fe] ratios in NGC 1313, to that in the [Ti/Fe] abundances. However, we also highlight that the uncertainties in the Mg abundances are much higher than those from Ti, particularly for NGC 1313-505. For this cluster we infer [Mg/Fe] = +0.70 ± 0.7 dex, with abundances compatible to those in the field in the MW and LMC, within the uncertainties.

Fe-peak elements
In Figure 5 we show the NGC 1313 YMC abundances for the Fe-peak elements Cr and Ni, and compared them with those observed in stars in the MW (Reddy et al. 2003(Reddy et al. , 2006 and the LMC (Van der Swaelmen et al. 2013). Both of these Fe-peak elements are primarily produced by Type Ia SNe. However, in spite of sharing a common origin, they exhibit slightly different patterns.
Cr abundances, both in the LMC and the MW, display a much higher scatter than that observed in Ni. The stellar abundances in these two environments are typically flat around [Cr/Fe]= 0.0 dex. The [Cr/Fe] ratios measured for our sample of YMCs in NGC 1313 appear to follow a comparable trend to those observed in the MW and LMC. The only cluster that appears to slightly deviate from this trend is the most metalpoor YMCs, NGC 1313-379, studied in Hernandez et al. (2017), displaying a subtle Cr enhancement.
Ni also exhibits a flat distribution in both the MW and the LMC. The trend observed in the LMC, however, is consistently sub-solar. Overall, the Ni abundances measured in NGC 1313 follow the flat trends of the MW and LMC, again with the exception of the most-metal poor YMCs, as well as the most-metal rich target.

Stellar Abundance gradients
Over the last couple of decades, studies investigating the chemical variations as a function of galactocentric distances in spiral galaxies have provided essential constraints on the evolution of these systems. The majority of these studies, greatly relying on oxygen abundances of H II regions, hint at an inside-out growth of galactic discs (Pilyugin et al. 2014;Bresolin & Kennicutt 2015;Ho et al. 2015;Lian et al. 2018). In this section we present our results on the stellar abundance gradients for Fe, and the α elements measured in this work.   Reddy et al. (2003Reddy et al. ( , 2006. Orange stars show the LMC abundances by Van der Swaelmen et al. (2013).
In Figure 6 we show in different panels the abundances as a function galactocentric distance for the individual elements, Fe, Ca, Ti, and Mg. We compute linear regressions for each element and show the inferred slopes with a solid line. Given that the theorems underlying least squares regressions assume asymptotic normality (meaning large-N samples), a condition that does not apply to our limited sample, we estimate our abundance gradients and their uncertainties applying a linear regression with non-parametric bootstrapping, where we apply the bootstrapping on the residuals, not the fitted parameters themselves. Briefly summarized, we find the optimal linear regression on our original dataset, extract the residuals from the best fit, generate new equal-size samples using the residuals and fitting linear regressions on the new samples. This resampling method allows for a more accurate estimate of the uncertainties in our inferred parameters. We measure an Fe gradient of −0.124 ± 0.034 dex kpc −1 . We find a slightly shallower gradient for Ca, −0.091 ± 0.026 dex kpc −1 . And finally, the gradients for Ti and Mg show a much larger scatter, with values of −0.117 ± 0.044 dex kpc −1 and −0.069 ± 0.089 dex kpc −1 , respectively. We estimate a weighted average αabundance gradient of −0.096 ± 0.009 dex kpc −1 , which is only slightly flatter than that observed for Fe. It is important to highlight this subtle difference as it is expected that these elements, Fe-peak and α, with different nucleosynthetic origins, imprint a different signature in the observed gradients. For example, in the past such a behaviour has been predicted and observed in the MW (Chiappini et al. 2001;Palla et al. 2020). Overall, negative abundance gradients are obtained when the star formation has been more efficient in the inner regions of the galaxy, compared to those in the outer regions of the disk, which based on our study, appears to be the case in NGC 1313.
Lastly, we perform a final test as the abundance gradients in Figure 6 appear to be primarily driven by the individual point at R/R 25 ∼0.9. We remove this individual point and estimate new gradients with the remaining points at R/R 25 <0.8. Overall we find that the abundance gradients for most elements remain within the uncertainties quoted in Figure 6, with the exception of Mg. The gradient for this particular element becomes slightly positive with much larger uncertainties, strongly highlighting the effects of the large scatter on the inferred Mg gradient.

Comparison to ionized-gas abundances
Past studies have investigated the distribution of metals in NGC 1313. Pagel et al. (1980) analyzed the observations of six H II regions and found a uniform gas-phase abundance distribution, with no clear gradient. Similarly, through the spectroscopic analysis of a much larger sample of 33 H II regions, Walsh & Roy (1997) reported a flat gas-phase abundance distribution across the disk. These results made NGC 1313 the most massive barred galaxy to lack a radial abundance gradient in the ionized-gas component. We highlight that these abundance studies relied on strong-line calibrations, which are based on the flux ratios of some of the strongest forbidden lines. It is well known that different strong-line calibrations return different abundance gradients with systematic offsets in the measured abundances (Kewley & Ellison 2008;Maiolino & Mannucci 2019), particularly in high-metallicity environments (Bresolin et al. 2016). Gas-phase metallicities estimated using the T e -based method are considered to be more reliable than those inferred from strong-line calibrations (Maiolino & Mannucci 2019). In Figure 7 we compare our stellar weighted average α-abundance gradient to the much shallower gradient inferred from the H II regions studied in Walsh & Roy (1997, herein after WR97), −0.016 ± 0.008 dex kpc −1 . The nebular abundances were estimated using the strong-line [O III]/[N II] calibration by Pettini & Pagel (2004) along with the measured fluxes by WR97. From this figure, it is clear that the overall trends imprinted in the stellar and gas-phase abundances are remarkably different. A puzzling effect as both of these components are expected to describe the present-day metallictiy distribution of NGC 1313. Studies such as that by Bresolin et al. (2016) and Hernandez et al. (2021) have previously shown that strong-line indicators produce nebular abundance gradients that are significantly shallower than those inferred from young stellar populations.
We expand our nebular and stellar abundance gradient comparison to include the direct abundances for six H II regions studied in Hadfield & Crowther (2007). In order to perform a meaningful comparison between the direct oxygen abundances and those measured from our YMC sample, we must correct the oxygen gas-phase metallicities for depletion onto interstellar dust grains. Overall, studies estimate depletion factors for oxygen that are between −0.08 and −0. show with square symbols the depletion-corrected [O/H] abundances by Hadfield & Crowther (2007). We note that due to the limited size of their H II-region sample, Hadfield & Crowther (2007) Figure 7. Abundances as a function of galactocentric distances. We show with a blue solid line the linear regression for the gas-phase oxygen abundances inferred with the strong-line O3N2 calibration by Pettini & Pagel (2004), using the fluxes by Walsh & Roy (1997). We show with a dasheddotted orange line our stellar weighted average α-abundance gradient. The blue squares display the depletion-corrected direct oxygen abundances by Hadfield & Crowther (2007).
Age-metallicity relations are important for understanding the chemical evolution of individual galaxies. These types of relations are able to shed some light on correlations between the age of the stellar clusters and possible interactions or close encounters with other systems. Livanou et al. (2013) studied the age-metallicity relation in the LMC by analyzing a sample of small open clusters. They report a metallicity gradient with higher abundances for the younger clusters. Such a trend is of course expected in the normal evolution of stellar populations in galaxies. Livanou et al. (2013) also identify a considerable increase in metallicity around the ages of ∼600 Myr, which they propose could be the result of the most recent encounter between the LMC and the SMC.
Our sample of YMCs in NGC 1313 cover ages between 20 and 300 Myr. In Figure 8 we display the agemetallicity trend for NGC 1313 with orange triangles, and compare it to the relation observed in the LMC (shown as blue squares) for clusters with comparable ages. Overall, there is a clear trend where clusters in NGC 1313 with ages ≥100 Myr began to show a steady decline in their metallicities. One YMC appears to deviate from this subtle trend, NGC 1313-379, with an age of 55 Myr and a metallicity of contrast to the rest of the galaxy, the bar and arms, which exhibits a constant star formation history. Their findings support a scenario where the south-west region of NGC 1313 is experiencing an interaction with a relatively small companion satellite. Previous to the work by Silva-Villa & Larsen (2012), several other studies had found evidence of satellite interactions in this same location (Blackman 1981;Peters et al. 1994;Larsen et al. 2007). The age-metallicity relation observed as part of this work ( Figure 8) supports a scenario of constant star formation history throughout the galaxy, with a possible recent burst in star formation close to the region where NGC 1313-379 is located.

CONCLUSIONS
We carry out the most extensive detailed abundance analysis of a sample of YMCs in the barred spiral galaxy NGC 1313 to date. Our analysis relied on X-Shooter/VLT spectroscopic observations and the method and code developed by Larsen et al. (2012). The main results of our work are summarized as follows: • We obtain the first stellar abundance gradients for NGC 1313. We infer an Fe gradient of −0.124 ± 0.034 dex kpc −1 and a weighted average α gradient of −0.096 ± 0.009 dex kpc −1 . We highlight the importance of this subtle difference between the two values as Fe-peak and α elements have different nucleosynthetic origins.
• We compare our stellar abundance gradient to those in the literature for the gas-phase com-ponent. Past studies of H II regions have reported an absence of a metallicity gradient in NGC 1313. Comparing our stellar abundances with those inferred for the ionized gas using the T ebased method, we find a much better agreement than those studies using strong-line calibrations. We highlight, however, that to confirm the presence of a steep gradient in the gas-phase metallicities, the sample of H II regions with available features needs to be expanded to cover the inner-and outer-most regions of the galaxy.
• We investigate the age-metallicity trends in NGC 1313. We observed a general trend where the metallicity declines with age. YMC NGC 1313-379 appears to slightly deviate from this trend, displaying a lower metallicity than that expected for a 55 Myr cluster. Our work supports a scenario of constant star formation throughout the galaxy, with a possible burst in the south-western region, where NGC 1313-379 is located.
Integrated-light abundance analysis continues to be a critical tool to study the chemical enrichment histories of galaxies outside of the Local Group. We highlight that abundance studies comparing the metallicities from different components, e.g., stellar and multi-phase gas (neutral and ionized; Hernandez et al. 2021), are essential to understand the multiple mechanisms driving the evolution of galaxies.

APPENDIX
In this section we present a complete set of modelled spectra which includes our final measured abundances. In Tables 6 to 10 we list the derived abundances for the different wavelength bins. Observed Wavelength (Å)

Facilities
Flux (10 16 erg s 1 cm 2 Å 1 ) Figure 13. Example synthesis fits for NGC 1313-463 (UVB X-Shooter coverage). In black we show the best-fitting model spectra. These models use the final abundances obtained as part of this work.  Flux (10 16 erg s 1 cm 2 Å 1 ) Figure 18. Example synthesis fits for NGC 1313-505 (VIS X-Shooter coverage). In black we show the best-fitting model spectra. These models use the final abundances obtained as part of this work.