Ages and metallicities of globular clusters in M81 using GTC/OSIRIS spectra

We here present the results of an analysis of the optical spectroscopy of 42 globular cluster (GC) candidates in the nearby spiral galaxy M81 (3.61~Mpc). The spectra were obtained using the long-slit and MOS modes of the OSIRIS instrument at the 10.4~m Gran Telescopio Canarias (GTC) at a spectral resolution of $\sim$1000. We used the classical H$\beta$ vs [MgFe]$'$ index diagram to separate genuine old GCs from clusters younger than 3 Gyr. Of the 30 spectra with continuum signal-to-noise ratio $>10$, we confirm 17 objects to be classical GCs (age $>10$~Gyr, $-1.4<$[Fe/H]$<-$0.4), with the remaining 13 being intermediate-age clusters (1-7.5~Gyr). We combined age and metallicity data of other nearby spiral galaxies ($\lesssim18$~Mpc) obtained using similar methodology like the one we have used here to understand the origin of GCs in spiral galaxies in the cosmological context. We find that the metal-poor ([Fe/H]<$-$1) GCs continued to form up to 6~Gyr after the first GCs were formed, with all younger systems (age $<8$~Gyr) being metal-rich.


INTRODUCTION
Globular Clusters (GCs) are among the oldest objects in the universe, which make them a key component to understand the formation and assembly history of galaxies (Brodie & Strader 2006).The availability of wide-field digital detectors on ground-based telescopes and the increased spatialresolution offered by the Hubble Space Telescope (HST) in the nineties opened up the study of GCs systems in external galaxies (Harris 1991).Unlike the GCs in the Milky Way, which are selected based on their morphological appearance, the samples of GCs in external galaxies are defined in terms of the photometric properties of compact sources in images (see e.g.Chandar et al. 2004;Muñoz et al. 2014;González-Lópezlira et al. 2017;Lomelí-Núñez et al. 2022).
⋆ E-mail: luislomeli@ov.ufrj.br(LLN) At present, there exists good quality photometric studies of large samples of GCs and spectroscopic studies for a limited number of GCs in elliptical and spiral galaxies.
(see e.g.Alves-Brito et al. 2011;Ko et al. 2018;Wang et al. 2021;Kim et al. 2021;Escudero et al. 2022).A common result found in these studies is the bimodal distribution of colours (Zepf & Ashman 1993, Gebhardt & Kissler-Patig 1999, Larsen et al. 2001).A correlation between colour and metallicity was also observed in different galaxies (e.g.Peng et al. 2006;Alves-Brito et al. 2011).Given this correlation, the bimodal colour distribution has been attributed to a bimodality in the abundance of metals (see, Brodie & Strader 2006).The bimodality would indicate the presence of two different populations of GCs, the metal-poor (commonly referred to as blue), and the metal-rich (or red).The two kinds of GCs also seem to be spatially segregated, with the distribution of metal-poor blue GCs having a larger scale length  1.The GC candidates from S10 are identified by red circles.The slitlets in the MOS fields that do not include GCs are SSCs and suspected Fuzzy Clusters, which are not part of the discussions in this study.
Two interpretations are in common use to explain the origin of the bimodality of metallicity in GC systems: a) a major-merger scenario proposed by Ashman & Zepf (1992), where the metal-poor GCs are accreted by progenitor merging spiral galaxies and the metal-rich GCs are formed insitu in metal-enriched gas following the merger event, b) an accretion scenario (Côté et al. 1998) invoking a massmetallicity relation, where the metal-rich GCs were formed in situ within a massive seed galaxy, while the metal-poor GCs were formed in less-massive satellite galaxies which are subsequently acquired by the host galaxy during the accretion process.In the major-merger scenario, the metal-poor and metal-rich GCs were formed before and after the merger event, respectively, and hence the metal-poor GCs are expected to be systematically older than the metal-rich GCs.This model predicts an age spread among the GCs.More importantly, clusters of intermediate-age (∼5-10 Gyr) are expected to be present among the objects classified as GCs from photometric data.On the other hand, in the accretion scenario both the metal-rich and metal-poor GCs were present in the host galaxies before the merger event.Under such circumstances, the relative ages of the metal-rich and metal-poor GCs depend on the formation epochs of lowmass GCs and massive halos (see, Peebles 1984;Rosenblatt et al. 1988).
Most of the analysis of spectroscopic data of GCs hitherto have been focused on obtaining their metallicity.Modern techniques of analysis of spectral data of unresolved stel-lar populations allow the determination of ages and metallicities with enough precision to be able to distinguish the predictions of these two cosmological scenarios of galaxy formation.We use new spectroscopic data of the GC systems in M81, a nearby giant spiral (Sab) galaxy located at 3.61 Mpc (Tully et al. 2013), along with available data on other GC systems, to look for the presence of intermediate-age clusters and also explore the age-metallicity relation in GC systems.
The population of GCs in M81 has been studied in the past by several authors.Perelmuter & Racine (1995) found ∼70 objects classified as cluster candidates in the inner 11 kpc.In a similar way Nantais et al. (2010a) and Santiago-Cortés et al. (2010) studied the star clusters population in M81 using the available data from the Hubble Space Telescope/Advance Camera for Surveys (HST/ACS), which consisted of 29 adjacent fields covering a field of view of ∼340 arcmin 2 with a sampling of 0.05 arcsec pixel −1 (0.88 pc pixel −1 ), reporting 233 and 172 GC candidates, respectively.More recently, Chies-Santos et al. ( 2022) studied a sample of GC candidates around the M81 group (M81, M82 and NGC 3077), using wide-field ground-based images, reporting 642 new GC candidates in a region of 3.5 deg 2 around the triplet.Some of the GC candidates have been targets of follow-up spectroscopic observations.Perelmuter et al. (1995) analyzed spectra of 82 GC candidates using the relative strengths of Hδ, CaI λ4227, and FeI λ4045, confirming 25 of them as bonafide GCs.Nantais et al. 2010b, obtained spectra of 74 GCs, and determined metallicities for GCs using an empirical calibration based on Milky Way (MW) GCs.They found a mean metallicity of [Fe/H]=−1.06± 0.07 dex.
The main goal of this work is the determination of age and metallicity for the GC candidates in M81 using spectroscopic data with the OSIRIS instrument at the 10.4-m Gran Telescopio Canarias (GTC).In §2, we describe the sample of GC candidates, the spectroscopic data, reduction and extraction.We explain the age and metallicity determination methods in §3.In §4, we carry out a detailed analysis of the determined age and metallicity of the GC candidates.In §5, we combine the M81 data with similar data from the literature for other spiral galaxies to address the age and metallicities of the GCs in spiral galaxies in the cosmological context.We carried out a spectroscopic campaign to observe  Gómez-González et al. 2016).Some slits passed through nebular regions which helped to obtain nebular abundance and its gradient in M81 (see Arellano-Córdova et al. 2016).Spectra of the brightest GC, M81-GC1 was analyzed by Mayya et al. (2013).

Spectroscopic observational details
Our sample has spectra of 42 of the 172 GC candidates reported in S10.All observations were carried out with the OSIRIS spectrograph at the Nasmyth-B focus using the R1000B grism with a slit-width 1.0 to 1.23 arcsec, covering a spectral range from 3700 to 7500 Å with a spectral resolution of ∼7 Å measured at λC = 5455 Å.We selected the observation mode with a binning of 2×2, obtaining a spatial scale of 0.254 arcsec pixel −1 (in the horizontal axis), and a spectral sampling of ∼2 Å pixel −1 (in the vertical axis).
In Table 2, we tabulate the photometric and spectroscopic parameters, of our sample of GC candidates.Object IDs (column 1) are taken from Lomelí-Núñez et al. (2022) (hereafter L22).Coordinates (RA, DEC) in J2000 epoch (column 2-3), B-magnitude (column 4) and (B − I)0 colour (column 5) taken from S10, (u − g)0 colour (column 6) are taken from L22, observed velocity (column 7), rootmean-square (RMS) and signal-to-noise ratio (SNR) of the continuum are measured in our spectra at 4190 Å (column 8-9).Quality (Q) of spectra of each GC candidate, which is obtained by visual inspection and is mainly based on a combination of the SNR, sky subtraction accuracy and position of the object in the slitlet (10).A more detailed description about Q parameter is given in Section 2.3.Observation run for each GC candidate (11), and its name from S10 (column 12) are also given.

Data reduction and spectral extraction
Both the long-slit and MOS spectral images were reduced using GTCMOS2 , which is a dedicated pipeline devoted to the reduction of OSIRIS/GTC observational blocks (OBs).The main steps in the reduction process are the following: (i) joining the two CCD images into a single mosaic image; ii) bias subtraction; iii) wavelength calibration and iv) flux calibration using standard stars (column (10) of Table 1).Wavelength calibration was carried out for each slitlet using the arc lamp images using the same slit mask as that used for source observations.The final product is a 2D-spectra corrected for geometric distortions, wavelength and flux calibrated.For a more detailed description of the software and reduction method see Gómez-González et al. (2016).
Given the relative faintness of the extragalactic GCs compared to the combined brightness of sky and disk background, it is crucial for the purpose of this study to maximise the SNR ratio while extracting 1-D spectra from the wavelength and flux calibrated 2D-images.Hence each of the 42 spectra was extracted interactively using the IRAF task apall.In all cases, the continuum was strong enough for a reliable tracing of the spectrum.Extraction of the spectra was carried out using a fixed window size of ∼4 pixels (1 arcsec) centered on the trace axis.For each object, a sky+background spectrum is extracted in windows of 3  (1) GC candidate name.(2,3) Right ascension and declination coordinates (J2000).( 4) B-magnitude from HST-F 435W band corrected for galactic extinction.
to 6 pixel widths on either side of the object, which is then scaled to the size of object extraction window and subtracted to get pure GC spectra.For MOS observations, the background window was chosen in slitlets specially reserved for that purpose, when the slitlet containing the object did not have enough background pixels.In either case, we ensured by examining visually the HST/ACS multiband images that the surface brightness of the disk in the zone selected for the extraction of the sky+background spectrum is as close to that in the vicinity of the target GC.
In Figure 2 we show an example of an extracted spectrum of a GC candidate, where we indicate the position of some of the prominent features observed in GC spectra.For the sake of clarity, the spectrum is shown in two sections: the bluer range (3900-6000 Å) which has most of the GC features (CaIIK,CaIIH,CNR,MgH,Mg2,Mgb,Fe52,Fe53,and F54), and a red window (6450-6800 Å) for the Hα line.The undisplayed middle part is featureless in all cases.The Hα and Hβ are used to obtain an average radial velocity of each GC candidate.All spectra shown in this work are in the rest frame after correcting for the Doppler shift using the average radial velocity.

Spectral quality
In the whole sample we have 9 objects with SNR at 4190 Å<10, 11 objects with 10≤SNR≤20, and 22 with SNR>20.In Figure 3, we show examples for three values (low, intermediate and high) of continuum SNR.The spectra were normalized to their flux at 5500 Å and vertically shifted adding an arbitrary constant.For the medium and low SNR spectra, we also show smoothed versions of the spectra.It can be noted that the absorption features can be identified in the smoothed spectrum even though they are buried in the noise in the case of low SNR spectrum.On the other hand, absorption features are identifiable without the need for smoothing for SNR>10.We assigned a quality index of Good (G) or Bad (B) for each extracted spectrum, which is based on the SNR, accuracy of sky subtraction and  the positioning of the object in the slitlet.For a spectrum to be assigned a G quality, it has to have an SNR of at least 10.All spectra with SNR<10, and some spectra with SNR>10, are assigned a B quality if the sky subtracted spectrum showed unphysical continuum shape (eg., wavy nature, negative fluxes etc.), or the object was at the edge of a slitlet, in which case the spectrum does not cover the entire targetted wavelength range.In column (10) of Table 2 we list the spectral quality for each GC spectrum.In total, spectra for 12 GCs with Quality B are not useful for our analysis.Thus the spectroscopic sample consists of 30 objects with Quality parameter G.

ANALYSIS
In this section, we describe in detail the methods we have used to determine the metallicity and ages of the GC candidates.

Hβ † vs [MgFe] ′ grid method
To begin with, we follow the classical grid method suggested by Thomas et al. (2004) (see also Brodie & Huchra 1990) in the Hβ † vs [MgFe] ′ for a first approximation of age and metallicity.This involves measurement of hydrogen (Hβ † ), magnesium (Mgb † ) and iron (Fe5270 † and Fe5335 † ) line strengths which were measured using the Lick/IDS index definitions from Trager et al. (1998) for equivalent width: where FI is the mean flux measured in the center of the index (column (3) in Table 3), FBC and FRC are the mean flux in the continuum bandpass in both sides of the index, blue continuum (BC, column (2) in Table 3) and red continuum   (RC, column (4) in Table 3).The spectral ranges of the indices measured with this definition are listed in the last four rows of Table 3. Thomas et al. (2003) realized that the [MgFe] ′ is independent of [α/Fe] and is a good tracer of metallicity.This index is determined as:

GC Hβ
(2) We used those indices for constructing the index-index diagram for a first-approximation age and metallicity values.
The calculated indices for all our GC candidates are tabulated in Table 4. Error on each index was calculated by performing 1000 Monte Carlo simulations using the rms errors on the fluxes as the sigma of the Gaussian error distribution.In Figure 4, we show our sample of GC candidates in the Hβ † vs [MgFe] ′ diagram, where we also plot the theoretical age-metallicity grids using the SSP models of Bruzual & Charlot (2003).Most of the GC candidates in the figure have Hβ † < 2.5 Å, and are located in the grid where the SSP ages are older than 8 Gyr and metallicity sub-solar, similar to the values for the Galactic GCs.Hence, for these objects, we obtained the metallicities using the empirical relation between the Fe index and metallicity suitable for the Galactic GCs, as explained in the next subsection.However, there are 8 GC candidates whose location in Figure 4 suggests that they are relatively young as compared to the classical Galactic GCs to be able to apply the Galactic Fe index-metallicity relation.Six of these have Hβ † ≥ 2.5 Å, suggesting ages less than ∼3 Gyr.The remaining two objects (GC29, GC136) show Hα in emission, which is most likely originating from the diffuse ionized gas in the disk, suggesting that the Hβ † line is contaminated by the nebular line and hence for these the determined Hβ † index is a lower limit, and hence likely to be younger than classical Galactic GCs.For these eight objects, we assign the metallicity value of the closest SSP models.
Three of these (GC44, GC60, GC114) have [MgFe] ′ ≥ 2.5 Å, which corresponds to super-solar metallicity SSP models ([Fe/H]=0.288).The remaining GCs are closest to SSP tracks corresponding to [Fe/H]=−1.249(GC11, GC87), −1.647 (GC47) and 0.066 (GC29, GC 136).The SSP metallicity obtained from the grid method is expected to be less accurate than that obtained using the Fe indices described below.Hence, the eight objects to which we assign their metallicity using the grid method are distinguished from the rest by a separate symbol in the corresponding Figures, and are also marked by a dagger in Tables 5 and 7.

Metallicity from Fe indices
We calculated the well-known Lick indices (Burstein et al. 1984) for all the GCs that have quality parameter G in Table 2.For the calculation of spectral indices of our sample of GC candidates observed with GTC/OSIRIS (R1000B), we used the definition from Brodie & Huchra (1990): where FI is the mean flux measured in the center of the index (column (3) in Table 3), FBC and FRC are the mean flux in the continuum bandpass in both sides of the index, blue continuum (BC, column (2) in Table 3) and red continuum (RC, column (4) in Table 3).The first eight indices listed in Table 3 are measured using Eq. 3, whereas the remaining four are measured according to the definition of Trager et al. (1998).Error on each index was calculated performing 1000 Monte Carlo simulations using Eq. 3, by taking into account the uncertainties of the fluxes FI , FBC and FRC .The calculated spectral indices, along with their errors, are listed in Table 5.We used the empirical relations between three iron indices (Fe52, Fe53 and Fe54) and the metallicity used in Mayya et al. (2013) to infer the chemical abundance of the GC candidates.The behaviour of each spectral index as a function of the metallicity is well described by a second order equation: where k, p and h are empirically derived coefficients.The fits are displayed in the windows: Fe52, Fe53 and Fe54 in Figure 6.To calibrate the above relationship, the MW-GCs (Schiavon et al. 2005) and SSP models were used.Because the SSP spectrum (see Section 3.3 for SSP description) have a higher resolution than the observed, the indices and metallicities determination were obtained after degrading the model to the same wavelength step of the observed spectrum.Once our observed and SSP spectra were matched in resolution and sampling, we used the Eq. 4 (coefficients k, p and h are given in Table 6) to calculate the metallicity.
The metallicity reported in Table 5 is the weighted mean from the three iron indices.The error in [Fe/H] reported in Table 5 was calculated propagating the error of the Fe indices, along with the uncertainty of the coefficients: k, p, and h.
In Figure 5, we compare the metallicity obtained here with those determined by Nantais et al. (2010b) (hereafter N10) for 12 GC candidates that are common between the two samples.The two measurements agree between each other (solid line of unit slope), within the quoted errors.However, we note here that N10 reported a metallicity using the indices even for some of the candidates (eg.GC87) for which we have measured Hβ index>2.5.As noted in the previous section, these GC candidates are younger than 5 Gyr and the index-metallicity relations are not reliable.The [Fe/H] values measured by us as the mean of the Fe indices are plotted against all the eight indices (CNR, G, MgH, Mg2, Mgb, Fe52, F53, and F54) in Figure 6 (black dots).We also show in this plot the corresponding values reported by N10 as diffuse red dots.Both the datasets show the expected correlation, with our measurements having, in general smaller spread, as compared to that of N10.Also in last panel of Figure 6, we show the [Fe/H] vs I-band magnitude, from which it can be inferred that we have [Fe/H] measurements for candidates that are 1 magnitude fainter, and the dispersion is smaller over the entire range of magnitudes as compared to that of N10.

Spectral fit: age and extinction determination
We determined the age for each GC candidate by selecting an SSP spectrum that best fits our observed spectrum.The fitting is carried out using a χ 2 procedure, which is explained For a given GC candidate, we used the SSP model that has metallicity closest to the [Fe/H] value measured in the previous section.These model spectra are smoothed to match the resolution of the observed spectra of ∼7 Å.Both the observed and model spectra are then divided by average flux in the wavelength range from 5490 to 5510 Å, to get the normalized spectra.Subsequently, the spectra were interpolated linearly in order to have the same sampling in the observed λ range.The χ 2 -fitting is carried out in intervals of ∆λ ∼ 300 Å over the spectral range I λ = [3650, 5600] Å.The spectra redward of 5600 Å do not contain age-sensitive absorption features, which is the reason for excluding this part in our χ 2 analysis.The piece-wise, rather than the entire spectral, fitting procedure ensures that the χ 2 values are not affected by differences in the shapes of the observed and model spectra over scales larger than ∆λ.Such differences are common due to the wavelength-dependent flux calibration errors, and/or due to residual errors in the subtraction of the sky in case of fainter objects.The χ 2 values obtained in the adopted window size is only weakly dependent on the value of the extinction, which allows us to obtain AV in an iterative manner as explained later in this section.
Then, for a given SSP model, we compute the reduced χ 2 k value within each k th -window using the formula where n k is the number of spectral points with wellmeasured fluxes, F obs,k is the observed spectrum, F mod,k is the model spectrum, σ obs,k is the standard deviation of the flux and ν = n k − 1 is the number of degrees of freedom.We defined masks to ignore spectral regions that are affected by the subtraction of bright sky lines.The value of n k is the number of unmasked pixels within the window.The σ obs,k is given by where F obs is the mean of observed flux in the k th window.
The value of σ obs,k is obtained iteratively, clipping the observed points that lie above and below 3σ and re-calculating a new F obs,k , F obs and σ obs,k in each iteration.The σ obs,k values obtained after ∼5 iterations do not change much, and hence we use the value after 5 iterations as the typical standard deviation of the spectrum in the window under analysis.4).‡Contains two populations.The two ages correspond to that of the metal-rich (young) and metal-poor (old) populations that best fit the spectra, with the [Fe/H] in column 2 corresponding to that of the older of the two populations.
The final value of χ 2 over the entire spectral range is calculated as the average of all χ 2 k .This process is repeated for all ages of the SSPs between 1 and 13.75 Gyr and the model age with the minimum χ 2 is taken as the best age of the GC candidate.
For the three cases (GC11, GC47, GC87) that have EW(Hβ)>2.5Å and sub-solar metallicities, the residual clearly shows Hβ in absorption suggesting a necessity for a second younger population.We hence found a new χ 2 solution where the input SSP is a linear combination of an old, metal-poor and young, solar metallicity populations.The χ 2 values for the combined SSP model were found to be lower than with only a single population for these three GCs.
We clarify that the extinction differences within the 300 Å window used for obtaining χ 2 k are ignored in Eq. 5, which is justifiable given the relatively small wavelength interval.However, we take into account the extinction correction in an iterative way.The observed spectrum over the entire observed wavelength is systematically redder as compared to the best-fit model spectrum.We used the slope between λB=4400, the normalized wavelength, λV , and λR=6400 Å, to obtain the AV using the following formula: A where the F B obs , F R obs , F B mod and F R mod are the average values over ∼200 Å windows centered at λB and λR measured in the observed and best-fit model spectrum, respectively, and EMW is the MW extinction curve from Cardelli et al. (1989).The A B V and A R V values are averaged to get the mean AV and error in the corresponding AV .All model spectra are reddened using the determined AV , normalized to its flux at λV and the entire fitting-procedure to obtain the minimum χ 2 is repeated.The iterative procedure is stopped when the age and AV values in two consecutive iterations converge.We find that the solutions converge after 2 and 5 iterations for high and low SNR spectra, respectively.
The ages of the SSP model and extinction for the bestfit case are tabulated in Column 3 and 5 in Table 7 respectively.According to the χ 2 -statistics all models that have ν(χ 2 − χ 2 min ) = 1 are acceptable solutions within 1−σ with 1 parameter involved (Wall & Jenkins 2012).We used this criterion to estimate the error on the best-fit ages.The possible range of ages, which are obtained by applying the lower and upper errors, are given in Column 4.
In Figure 7 we show the fitted spectrum for one of the representative GCs, GC16.The measured [Fe/H]= −0.54 ± 0.08 for which we used the [Fe/H] model = −1.2486.The bestfit age=12.75+1.0 −0.25 Gyr and AV =1.32±0.02.The residual is shown in the bottom panel.In order to illustrate the goodness of the fit, we zoom in on the Ca H+K feature, and the Hα.It can be appreciated that the best-fit SSP for the determined AV reproduces the observed strength of the Hα line as well as the continuum level very well, in spite of the fact that the fitting was carried out blueward of 5600 Å.

DISCUSSION
Our spectroscopic sample consists of 42 GC candidates, among which 30 have good quality spectra for a reliable determination of their metallicities, ages and extinctions.Spec-  Bruzual & Charlot (2003) with Kroupa IMF between 0.1 and 100 M ⊙ are over-plotted for four fixed ages between 1 and 14 Gyr and six metallicities between −1.64 ≤[Fe/H]≤ +0.28.The solid coloured lines connecting grid points at fixed age are second order polynomials fitted to the SSP models.The GCs whose metallicity was estimated using the grid method from Figure 4 are encircled in magenta.
tra of all the 42 however are good enough for the measurement of recessional velocity (see column 7 in Table 2).Being a sub-sample of the S10 sample, all also have well-measured colours and magnitudes.Given that GC populations of very few spiral galaxies have all these data available, we use these data to both characterise the properties of GC samples as well as to address the origin of GCs in a Cosmological context.In the first two sections, we discuss colour-metallicity relation and age distribution of our spectroscopic sample.In section 4.3, we combine the velocity of our sample objects with that from N10 to explore whether the majority show kinematics expected for halo or disk objects.In section 4.4, we combine our age and metallicity data with that from other samples that are analysed using methods similar to ours to discuss GCs from a Cosmological perspective.

Metallicity-colour relation
One of the most frequently discussed topics in the study of GCs is the bimodal colour distribution, which is thought to be originated from a metallicity-colour relation (see e.g.Brodie & Strader 2006).It may be recalled that our spectroscopic sample is drawn from S10, which is a photometric sample of 172 GC candidates covering the entire FoV offered by the HST/ACS footprints displayed in Figure 1.Our spectroscopic sample is not a representative subset of the original sample in colour and is relatively of small size, which prevents us an exploration of the issue of bimodality in metallicity.Besides, the original S10 (black solid histogram in top panel of Figure 8) sample itself does not demonstrate an obvious colour bimodality.We confirmed an absence of evidence of colour-bimodality in the S10 sample using the Gaussian Mixture Modeling (GMM) code (Muratov & Gnedin 2010) which carries out a robust statistical test for evaluating bimodality, and uses the likelihood-ratios to compare the goodness of the fit for double-Gaussian versus a single-Gaussian.This method is independent of the binning of the sample.The results from the GMM test in the colour distribution are shown in Table 8 (see the footnote for the explanation of the column headers).For a distribution to be considered bimodal, the Kurtosis must be negative, D> 2 and p-values must be small.We found a Kurtosis of 2.48 and a separation (D) of 1.42, values that strongly discard bimodality in the distribution.
We also used "dip test" (Hartigan & Hartigan 1985), which calculates the maximum distance between the cumulative input distribution and the best-fitting unimodal distribution.This test estimates the significance level with which a unimodal distribution can be rejected.The probability given by the "dip test" to be a unimodal distribution is of 89%.After carrying out these two independent statistical tests, we conclude that the M81 GCs of the S10 sample do not show evidence for bimodality in colour.This may not necessarily mean an absence of bimodality among all M81 GCs, as the S10 sample is more biased towards the disk GCs and lacks halo GCs, which are known to be systematically metal-poor and blue, as compared to the disk GCs.
Several studies have addressed the issue of a relation between metallicity and colour of GCs from both observational (e.g.Larsen et al. 2001;Alves-Brito et al. 2011;Fahrion et al. 2020;Kim et al. 2021) and theoretical (e.g.Yoon et al. 2006;Cantiello & Blakeslee 2007) perspectives.The availability of AV intrinsic to the objects in our spectroscopic sample allows investigating this relation with the reddening-corrected colours.The latter correction was carried out using the AV in column 5 of Table 7 and the Cardelli et al. (1989) reddening curve.In the bottom panel of Figure 8, we show the metallicity-colour relation for the spectroscopic sample, before (red circles) and after (blue circles) reddening correction.The colour distribution of the spectroscopic sample before and after reddening correction is shown in the top panel.We also show the theoretical age-metallicity grids using the SSP models of 1, 3, 8 and 14 Gyr (black, green, cyan and magenta triangles, respectively) from Bruzual & Charlot (2003), with Kroupa IMF (from 0.1 to 100 M⊙),  and six metallicities: Z=0.0004, 0.001, 0.004, 0.008, 0.019 and 0.030.At a given age, the metallicities and SSP model colours are well-fitted by second order polynomial curves, which are shown by solid lines.The SSP model curves for the two extreme ages (1 and 14 Gyr) enclose the distribution of observed points after dereddening the colours.Note that the observed colours before reddening correction for all objects are redder than the values expected for the oldest SSPs.We find a mean extinction of AV = 0.80 ± 0.13 mag for our spectroscopic sample of GCs in M81.This relatively high values illustrate the importance of reddening corrections before any robust statistical analysis of colours, such as the exploration of bimodality, colour-metallicity relation etc.
It is interesting to note that the uncorrected (i.e.observed) colours of the majority of clusters of our spectroscopic sample are redder than the SSP values expected for their spectroscopic ages.The reddening correction, however, makes the colour bluer for their derived ages.Errors in both the axis (horizontal error bars take into account the errors on our AV measurements) can account for some of this inconsistency.It is also likely that a hidden younger second population is making the colours bluer for the old GCs (solid blue dots).A detailed multi-colour study would be required to simultaneously explain the metallicity, age and colours for each object, which is beyond the scope of the current discussion.Notwithstanding, the figure suffices to illustrate the existence of an empirical metallicity-colour relation.

Age distribution of candidate GCs
Ages of GC candidates were determined using the χ 2 fitting between the observed and SSP spectra as discussed in section 3.3.In Figure 9, we show the age distribution of our spectroscopic GC sample.The distribution peaks at the oldest age bin, with 17 GCs being older than 8 Gyr.We consider these 17 cases as genuine "classical GCs" similar to the ones in the MW.The remaining 13 cases, including seven cases for which we estimated metallicity using the grid-method, have ages between 1 and 8 Gyr, which are clearly younger than the typical Galactic GCs.We refer to these 13 cases, 9 of which are younger than 3 Gyr, as "intermediate age star clusters".
Ages of our spectroscopic sample suggests that as much as 43% (13/30) of GC candidates are not classical GCs like those in the MW.The main reason for this relatively large contamination fraction is the age-extinction degeneracy which results in the reddened relatively younger clusters to have similar F 435W − F 814W colours as that of unreddened classical GCs.L22 used u − g vs F 435W − F 814W colour-colour diagram to break this degeneracy and determine the contamination of GC candidate samples from reddened young clusters (age <3 Gyr) for a sample of five spiral galaxies, including M81.They found a contamination fraction of 32% in M81.Our spectroscopic sample covered the very inner part of M81, where the extinction and contamination from younger populations are expected to be higher, which is the most likely reason for slightly higher contamination factor for the spectroscopic sample.

Velocities and spatial distribution
In early-type galaxies the distribution of metal-poor GCs extends to the outer halos while metal-rich ones are more concentrated to the internal parts (see e.g.Hargis & Rhode 2014;Kartha et al. 2014).However, this is not entirely proven for late-type galaxies, since studies of GCs in these kind of galaxies are scarce (see e.g.Hargis & Rhode 2014).Metal-poor outer halo GCs and metal-rich disk GCs also distinguish in their kinematical behaviour with the former being pressure supported, whereas the rotation dominates the kinematics of the latter.Testing these scenarios require GC Black ellipse is 0.5 times R 25 .We show our sample as big circles (6 ′′ ) and the N10 sample as small squares (3 ′′ ).The GCs whose metallicity was estimated using the grid method from Figure 4 are encircled in yellow.North up and east left.
samples with velocity and metallicity measurements that cover both halo and disk components.Unfortunately GCs in our spectroscopic sample reside in the inner parts of the galaxies, and lack halo GCs (for simplicity in our sample we called GCs to both intermediate age star clusters and classic GCs, unless a distinction is made between them).In order to overcome this bias, we complement the velocities and metallicities for the 30 GC candidates of our spectroscopic sample with these quantities for 108 GCs from N10 sample, which contains the outer halo GCs.The combined sample allows us to examine whether the metal-poor and metal-rich GCs differ in their kinematical properties.
In Figure 10, we show the spatial distribution of GCs for the combined sample superposed on a HST/ACS F 814W image.The metal-poor and metal-rich GCs are shown in separate colours.We also show an ellipse with semi-major axes a=6.725 ′ and b=3.525 ′ which encloses the region defined by 0.5R25 (Corwin et al. 1994).Unfortunately, even the combined sample has very few GCs outside the ellipse, and hence the most complete sample that can be constructed today for M81 is skewed towards disk GCs.
In left panel of Figure 11, we show the 1D phase-space diagram from this work (big circles) and N10 (small squares) samples of GC candidates, where the galaxy systemic velocity, −38.97±2.69km s −1 (Speights & Westpfahl 2012) is shown as a dashed line.Metal-poor and metal-rich GCs are shown in blue and red colours, respectively.In such diagram, the object location is related to the kinematics of the component it belongs to (i.e.disk or bulge/halo) and to the time since first in-fall into the galaxy halo (Rocha et al. 2012).In fact, the outer envelope in the R vs v obs space is dominated by recently accreted objects, while objects with smaller difference to the galaxy systemic velocity and lower distance to the galaxy center, have been accreted further back in time (Rocha et al. 2012).
From Figure 11 (left), it can be noted that the GCs lying farther than 10 kpc are mostly metal-poor, whereas at smaller distances both kinds of GCs are found.Majority of the GCs have velocities less than ±200 km s −1 .In the right panel of Figure 11, we show the velocity versus distance along the major axis.In this plot it is possible to observe that most of the GCs follow a pattern of disk rotation, i.e. they present a gradient in velocity typical of rotating systems.

Cosmological implications from metallicity vs age diagram for spiral and lenticular galaxies
Most of the age determinations for extragalactic GCs are based on the analysis of colours and Spectral Energy Distributions (SEDs), which have the advantage of studying the properties for the whole population.However, photometrically-based quantities suffer from age, reddening and metallicity degeneracy, which makes the obtained results unreliable.On the other hand, the availability in recent years for Multi-Object Spectroscopic (MOS) and/or Integral Field Unit (IFU) facilities at large telescopes covering reasonably large FoVs have enabled the determination of the ages without the effects of degeneracy.For example, the main part of our analysis is based on three MOS pointing observations each covering a FoV of 7 × 2 arcmin 2 (7.5 × 6.0 arcmin 2 ) at the 10.4-m GTC.We combine our measurements with spectroscopic measurements for other galaxies from the literature and analyse them in a metallicity-age diagram in Figure 12.The top panel shows the ages on a logarithmic scale, whereas the bottom panel is zoomed in on a linear scale for ages >8 Gyr to show metallicity-age trends for classical GCs.The sources of literature data are: MW (empty black triangles) and M31 (empty blue stars) data from Cezario et al. (2013), where they used the ULySS package (Koleva et al. 2009) to derive the metallicity and age determination, and M85 GCs (empty red squares) from Ko et al. (2018), where they used different methods for estimating age and metallicity, but for consistency, we only show the results that they obtained with the ULySS package.Also we show the M85 GCs (empty cyan squares) from Escudero et al. (2022), who also used ULySS package to derive the metallicity and age determination.We refer to the latter sample as M85*.A horizontal dashed line divides metal-poor ([Fe/H]≤ −1.0) GCs from the metal-rich ([Fe/H]> −1.0) GCs.
It can be noticed from the figure that the extragalactic GCs with spectroscopic ages are dominated by the metalrich population, with most of the measurements for the metal-poor GCs coming from the Galactic (from a sample of 41 GCs, 19 have an Age> 8 Gyr and [Fe/H]≤ −1.0) and M85* (from a sample of 46 GCs, 14 have an Age> 8 Gyr and [Fe/H]≤ −1.0) GCs, whose spectroscopic sample covered halo GCs faraway to the galaxy center.This is due to an observational bias originated due to a tendency to maximize the objects per pointing, which invariably leads to select zones of high object density for spectroscopic observa- In both plots we show our spectroscopic sample as big circles and the N10 sample as small squares.Also in both plots GCs whose metallicity was estimated using the grid method from Figure 4 are encircled in magenta.
tions.Such strategies end up sampling the disk GCs which are metal-rich.Differently, from other works, in Escudero et al. (2022) the outer part of the galaxies were observed, combining part of the disk and possibly halo of the galaxy.Interestingly, some of the GCs selected in that work present halo like characteristics, metal-poor and older than 12 Gyr.We can draw some general remarks from the combined sample.
Metal-poor GCs are concentrated in the age range from 8 to 12 Gyr, with a mean age that is ∼3 Gyr younger than that for metal-rich GCs.According to the Hierarchical scenario of galaxy formation (e.g.Côté et al. 1998;Strader et al. 2005;Forbes & Remus 2018), metal-poor GCs are formed in the halos of low-mass galaxies, which are then accreted onto the giant galaxies during subsequent merger processes.Under this scenario, metal-poor GCs are expected to be the oldest GCs (e.g.Peng et al. 2004;Beasley et al. 2008) which is opposite to what is inferred from the figure.The observed contradiction is opposite to what is expected if the spectroscopic determinations for some unknown reason suffer from age-metallicity degeneracy.On the other hand, the presence of metal-rich old (12-13 Gyr) GCs could be understood under the Cosmological context as it takes less than a 1 Gyr for metal enrichment in GCs.To better illustrate these statements, in the Figure 13 we show the GCs age distributions.In the top panel of Figure 13 we show the spectroscopic age distribution of 142 GCs, for M81, MW, M31 and M85 galaxies.An age bimodality can be seen at young and old ages.In the middle panel, the mean age distribution for metal-poor GCs ([Fe/H] <= −1), whereas the corresponding value for metal-rich ([Fe/H] > −1) GCs is 12.8 Gyr.The metal-rich GCs also show a peak at relatively a young age of 2 Gyr.
How do we reconcile the absence of metal-poor GCs among the oldest GCs?We believe the answer lies in the observational bias discussed above.The accreted old metalpoor GCs are expected to be in the outer halos of galaxies.In general, these outer-halo GCs are not included in the spec-troscopic samples, which explains the absence of old metalpoor GCs in Figure 12.The work of Escudero et al. (2022) mitigates this situation by finding several old, metal-poor GCs in the outer part of M85, suggesting that interesting insights on galaxy formation and its connection with its GC system will arise by an extensive study of GCs in galaxy outskirts, as well a systemic compilation of archival data.The presence of metal-poor, but marginally younger GCs among the spectroscopic samples suggests that not all the metal-poor GCs are formed at the earliest epochs.This can be explained under the hierarchical scenario of galaxy formation, if GCs formation in some of the accreting low-mass galaxies took place after a delay of a few gigayears.Another possibility is that these metal-poor GCs are formed in situ when gas-rich galaxies merge with a giant galaxy.The fact that these metal-poor GCs are found in spectroscopic samples that cover the inner parts of galaxies favours the later hypothesis.The acceptance of this hypothesis requires a change in the current paradigm that all in situ formed GCs are metal-rich.Surprisingly there are very few metalrich GCs of age between 12 and 8 Gyr.As well as, the formation of metal-poor GCs seems to have stopped suddenly at ∼8 Gyr.
M81 GCs follow another trend seen in other galaxies, namely the presence of intermediate-age clusters that have photometric and morphological characteristics similar to the classical GCs.Metal-rich objects dominate the sample at these ages.The metal-poor GCs are scarce at ages younger than 8 Gyr.This suggest that if metal-poor GCs in the inner disk are formed in situ during accretion of gas-rich dwarf galaxies, such process stopped suddenly at ∼6 Gyr after the Big Bang.This implies that most of the accreting galaxies have sufficiently enriched their metal content at this epoch.The accretion process and formation of GC-like clusters continued after that, but the newly formed clusters were mostly metal-rich.In both panels the horizontal black dashed line indicates the limit between metal-rich ([Fe/H]> −1.0) and metal-poor ([Fe/H]≤ −1.0) GCs.In top panel the vertical black dashed line indicates the limit between intermediate and old age clusters.Top: X-axis in logarithmic scale for all the range of ages.The GCs whose metallicity was estimated using the grid method from Figure 4 are circled in magenta.Bottom: X-axis in linear scale with ages from 7.6 to 18 Gyr.

CONCLUSIONS
In this work, we presented the results of the analysis of new spectroscopic data of 42 GC candidates in the nearby spiral galaxy M81, obtained with GTC-OSIRIS in the long-slit and MOS modes.Spectra for 30 of these objects have a S/N > 10 providing an accurate determination of metallicity, age and extinction.We used the classical Hβ vs MgFe diagram to separate clusters that are younger than ∼3 Gyr.For objects older than 3 Gyr, we used the iron indices (Fe52, Fe53 and Fe54) to determine the metallicity.We used the χ 2 -technique to determine the best-fit SSP model age of each GC candidate.We find that 17 GC candidates are genuine classical GCs with ages >8 Gyr.The remaining 13 candidates (43%) are intermediate-age clusters.M81 GCs experience a mean AV = 0.80 ± 0.13 mag which illustrates the importance of reddening corrections while inferring any physical quantity from optical colours and/or SEDs.We combine the spectroscopically determined age and metallicity of M81 GCs with those obtained using similar techniques for other samples of nearby galaxies to discuss the trends seen in the metallicityage diagram.
Most of the metal-rich ([Fe/H] > −1) GCs in spectroscopic samples of spiral galaxies are either older than 12 Gyr, or younger than 8 Gyr, with the age range of 8 to 12 Gyr being populated by metal-poor GCs.We postulate that these metal-poor GCs were formed in situ in the disks of galaxies during the accretion of gas-rich dwarf galaxies.Clusters younger than 8 Gyr are systematically metal-rich, suggesting the lack of accreting metal-poor dwarf galaxies at these relatively recent epochs.The absence in spectroscopic samples of old (>12 Gyr) metal-poor GCs is most likely due an observational bias of not including outer halo GCs, which are expected to be metal-poor.We recommend spectroscopic studies to include more outer halo GC candidates to address this question.

Figure 1 .
Figure 1.M81 footprint in the HST/ACS F 814W -band, showing the five long-slits and three MOS fields used for the data analysed in this work.Slits and MOS fields are identified by observing run code listed in Table1.The GC candidates from S10 are identified by red circles.The slitlets in the MOS fields that do not include GCs are SSCs and suspected Fuzzy Clusters, which are not part of the discussions in this study.
-Cortés et al. (2010) (hereafter S10) analyzed HST images of M81 using multi pointing observations with the Advanced Camera for Surveys (ACS) in F 435W (B), F 606W (V ) and F 814W (I) bands.The population of stellar clusters was classified in two subpopulations: young (blue) 3 Myr to 2 Gyr and old (red) ≳5 Gyr.They obtained 263 blue and 172 red star clusters.The selection criteria of stellar clusters was based on structural parameters obtained with SExtractor 1 (Bertin & Arnouts 1996): fwhm, area and ellipticity; and an evolutionary parameter, the photometric colour.

Figure 2 .
Figure 2. Example of the extracted spectrum of a GC candidate.In this case we show the GC named GC16.The observed flux (normalized at 5500 Å) from: i) 3900 to 6000 Å (blue line) and ii) 6450 to 6800 Å (red line) is plotted.

Figure 3 .
Figure 3.Comparison of spectral quality: spectrum with low SNR (red line), medium SNR (blue line) and high SNR (black line).In the cases of low and medium SNR we superposed their smoothed spectra in black color.

Figure 5 .
Figure 5. Metallicity comparison for M81 GCs.In the X-axis we plot our M81-GCs metallicities estimation, whereas in Y-axis are plotted the GCs metallicities estimation from N10.

Figure 6 .
Figure 6.[Fe/H] vs spectral indices.Black dots: metallicity and spectral indices for each GC candidate.Red diffuse dots: metallicity and spectral indices reported in N10.In the iron (Fe52, Fe53 and Fe54 vs [Fe/H]) windows the black solid lines are the curves tracing Equation4.In last panel we show the [Fe/H] vs I-band magnitude.In all the panels the GCs whose metallicity was estimated using the grid method from Figure4are encircled in magenta.

Figure 7 .
Figure 7. Observed spectrum (black line) and the spectral fit (red line).Gray bands indicate the wavelength coverage of the spectral indices.The fitted parameters are: [Fe/H]= −0.54 ± 0.10, Age=12.75Gyr, A V =1.32 and χ 2 = 2.53.The bottom plots show the residuals between the observation and the model.Interior windows are zooming the regions of the CaH+K (left) and Hα (right) lines.

Figure 8 .
Figure 8. Top: (B−I) 0 colour distribution.Black histogram: whole sample from S10 corrected for the Galactic reddening.Blue histogram: sample of GC candidates corrected for the Galactic and intrinsic reddening.Red histogram: sample of GC candidates corrected only for the Galactic reddening.Bottom: metallicity-colour relation for the sample of GC candidates with spectroscopic observations after (blue circles) and before (red circles) intrinsic reddening correction.Also the blue circles are coded by their spectroscopic ages as indicated by the legends at the bottom right corner.SSP model grids fromBruzual & Charlot (2003) with Kroupa IMF between 0.1 and 100 M ⊙ are over-plotted for four fixed ages between 1 and 14 Gyr and six metallicities between −1.64 ≤[Fe/H]≤ +0.28.The solid coloured lines connecting grid points at fixed age are second order polynomials fitted to the SSP models.The GCs whose metallicity was estimated using the grid method from Figure4are encircled in magenta.
Galaxy.(2-3) Mean and standard deviation of the first peak in the double-Gaussian model.(4-5) Mean and sigma of the second peak in the double-Gaussian model.(6) Total number of GCs.(7) Fration of N GC associated with the second peak.(8) Separation of the means relative to their widths.(9) GMM p-values based on the likelihood-ratio test p(χ 2 ), peak separation p(DD), and kurtosis p(kurt) (lower p-values are more significant).(10) Kurtosis of the colours distribution.(11) Dip value is the significance level with which a unimodal distribution can be rejected.(12) Bimodality final evaluation: Y (Yes), N(No).We note that it was only possible to fit 153 of the 172 GC candidates reported in S10.

Figure 9 .
Figure 9. Age distribution of GC candidates.The bulk of the sample has an old age typical of GC.The magenta histogram corresponds to GCs whose metallicity was estimated using the grid method from Figure 4.

Figure 11 .
Figure 11.Phase-space diagram of the GC candidates, with metal-poor [Fe/H] ≤ −1.0 (blue color) and metal-rich [Fe/H] > −1.0 (red color) distinguished.The observed velocity is plotted against galactocentric distance (left) and distance along the major axis (right).The horizontal dashed line in both panels is the galaxy systemic velocity.In right panel the vertical line shows the centre of major axis.In both plots we show our spectroscopic sample as big circles and the N10 sample as small squares.Also in both plots GCs whose metallicity was estimated using the grid method from Figure4are encircled in magenta.

Figure 12 .
Figure 12.Metallicity vs age for different galaxies: M81 (green solid points), MW (black empty triangles), M31 (blue empty stars), M85 (red empty squares) and M85 * (cyan empty squares).In both panels the horizontal black dashed line indicates the limit between metal-rich ([Fe/H]> −1.0) and metal-poor ([Fe/H]≤ −1.0) GCs.In top panel the vertical black dashed line indicates the limit between intermediate and old age clusters.Top: X-axis in logarithmic scale for all the range of ages.The GCs whose metallicity was estimated using the grid method from Figure4are circled in magenta.Bottom: X-axis in linear scale with ages from 7.6 to 18 Gyr.

Table 1 .
Log of the long-slit and MOS spectroscopic observations with GTC/OSIRIS in M81.
a subsample of star clusters catalogued by S10 using the OSIRIS instrument at the 10.4-m Gran Telescopio Canarias (GTC).The observational campaign consisted of three runs, one each in 2010, 2012 and 2014, totalling eight independent pointings.The footprints of these pointings are shown in Figure1, with the observational details for each pointing given in Table1.The 2010 and 2012 observations were carried out using long-slits whereas the 2014 observations were carried out in the MOS mode.The number of GCs NGC in each pointing is given in column 11 in Table1.The targets for observations included blue objects, which were super star cluster (SSCs), GCs and suspected Fuzzy clusters.The pointing 2014A-MOS3 was specially designed for the observation of GC candidates, where 25 GCs were observed.The pointing 2014A-MOS1 and 2014A-MOS3 included blue objects, which were SSCs and suspected Fuzzy clusters.In these different observational runs, 47 GC candidates were targeted.However, five of them were too close to the slit borders for a reliable spectral extraction.Hence, the final sample of GC candidates with extracted spectra contained 42 objects.Results obtained from spectra targeted towards SSCs helped to detect individal Wolf Rayet stars in M81 (see

Table 2 .
Sample of GC candidates.

Table 6 .
Parameters of the second order equation which was used to estimate the metallicity.

Table 8 .
Results from GMM and dip test.