The survey of planetary nebulae in Andromeda (M 31) V. Chemical enrichment of the thin and thicker discs of Andromeda. Oxygen to argon abundance ratios for planetary nebulae and H II regions

We use oxygen and argon abundances for planetary nebulae (PNe) with low internal extinction (progenitor ages of (>4.5 Gyr) and high extinction (progenitor ages<2.5 Gyr), as well as those of the H II regions, to constrain the chemical enrichment and star formation efficiency in the thin and thicker discs of M31. The argon element is produced in larger fraction by Type Ia supernovae (SNe) than oxygen. We find that the mean log(O/Ar) values of PNe as a function of their argon abundances, 12 + log(Ar/H), trace the inter-stellar matter (ISM) conditions at the time of birth of the M 31 disc PN progenitors. Thus the chemical enrichment and star formation efficiency information encoded in the [alpha/Fe] vs. [Fe/H] distribution of stars is also imprinted in the oxygen-to-argon abundance ratio log(O/Ar) vs. argon abundance for the nebular emissions of the different stellar evolution phases. We propose to use the log(O/Ar) vs. (12 + log(Ar/H)) distribution of PNe with different ages to constrain the star-formation histories of the parent stellar populations in the thin and thicker M31 discs. For the inner M31 disc (R_{GC}<14 kpc), the chemical evolution model that reproduces the mean log(O/Ar) values as function of argon abundance for the high- and low-extinction PNe requires a second infall of metal poorer gas during a gas-rich (wet) satellite merger. In M31, the thin disc is younger and less radially extended, formed stars at a higher star formation efficiency, and had a faster chemical enrichment timescale than the more extended, thicker disc. Both the thin and thicker disc in M31 reach similar high argon abundances ( 12 + log(Ar/H) ) ~ 6.7. The chemical and structural properties of the thin/thicker discs in M31 are thus remarkably different from those determined for the Milky Way thin and thick discs.


Introduction
Late-type galaxies may contain kinematically distinct components such as the "cold" thin disc and the "hot" thick disc found in the Milky Way (MW; e.g. Gilmore & Reid 1983) and in nearby galaxies (Yoachim & Dalcanton 2006;Comerón et al. 2019). The MW thick disc is structurally and chemically dis-Article number, page 1 of 14 arXiv:2208.02328v1 [astro-ph.GA] 3 Aug 2022 A&A proofs: manuscript no. main_03082022 tinct from the MW thin disc as well as in age. Differences in the structural parameters (stellar masses, exponential scale lengths, and velocity dispersion) of the MW thin and thick discs are summarized in Bland-Hawthorn & Gerhard (2016). The different chemical properties are most prominent in the measured stellar [α/Fe] ratios as a function of [Fe/H], with the old MW thick disc being more metal-poor and α-enriched, compared to the relatively younger MW thin disc (Hayden et al. 2015;Matteucci 2021). The α-enriched thick disk population is found to be confined within R 9 kpc, and older than ∼ 8 Gyr (Haywood et al. 2013;Belokurov et al. 2020), whereas the outer MW disc is composed of low [α/Fe] stars (Hayden et al. 2015). These properties are believed to have been set by the MW's most recent impactful merger ∼10 Gyr ago (Belokurov et al. 2018;Helmi et al. 2018), after which the MW disc is thought to have evolved mainly by secular evolution (see, e.g., Sellwood 2014).
Differently from the MW, M 31 had a more turbulent history, as vividly illustrated by the many substructures identified in its inner halo by PAndAS (McConnachie et al. 2009, including the Giant Stellar Stream (GSS, Ibata et al. 2001). Its most recent important merger is believed to have happened ∼ 2.5-4.5 Gyr ago (Bhattacharya et al. 2019b, hereafter Paper II). The M31 disc has a significantly steeper age-velocity dispersion (AVD) relation than that of the MW disc (Paper II; Dorman et al. 2015), with the velocity dispersion of the 2.5 Gyr and 4.5 Gyr old stellar populations being almost twice resp. three times those of the MW disc stellar populations of corresponding ages. Paper II used planetary nebulae as kinematic tracers to identify a younger (progenitor ages < 2.5 Gyr), dynamically colder disk, and a distinct, older (progenitor ages > 4.5 Gyr), dynamically hotter, hence thicker disc, with the latter having a velocity dispersion σ M31,thick 3 × σ MW,thick in the radial range 14-20 kpc (equivalent to the solar neighbourhood). At these radial distances (R GC =14-20 kpc), the 4.5 Gyr and older population of stars had velocity dispersion values ≥ 90 kms −1 , which are significantly larger even than the average velocity dispersion measured for strongly turbulent discs at redshift ∼ 1 − 2, 30 and 60 kms −1 , respectively, see Wisnioski et al. (2015). Mergers with satellites can dynamically heat thin discs, i.e., increase their velocity dispersion (Quinn & Goodman 1986) and decrease their rotational velocity, resulting in a thickened disc (Hopkins et al. 2009). Using their results, the AVD relation of the M 31 disc in a radial range R GC =14-20 kpc was found to be consistent with the energy injected in the M 31 disc by a major merger with mass ratio ∼1:5 ∼2.5-4.5 Gyr ago (Bhattacharya et al. 2019b), as predicted in the merger simulations of Hammer et al. (2018).
The different merger histories of these two spiral galaxies are seen in their large separation in the halo metallicity vs. total stellar mass diagram, where the MW and M31 are placed at opposite edges of the distribution of galaxies measured by the GHOST survey (Monachesi et al. 2019). In simulations, the spread in halo masses and [Fe/H] values is found to be indicative of different accretion histories (D'Souza & Bell 2018) and accreted satellite stellar mass. Independent confirmation of a recent major merger event in M31 can be sought through chemical abundances, altered by merger-related processes such as gas accretion and star formation bursts (e.g. Kobayashi & Nakasato 2011).
Planetary nebulae (PNe) are useful tracers to constrain the kinematics (Aniyan et al. 2018(Aniyan et al. , 2021 and chemical abundances (Magrini et al. 2016;Stanghellini & Haywood 2018) over a large radial range in nearby galaxies of different morphological types (e.g. Cortesi et al. 2013;Pulsoni et al. 2018;Hartke et al. 2022). PN elemental abundances shed light on the ISM conditions at the time of formation of their parent stellar population. When the PN ages are also constrained, it becomes possible to map abundance variations across different epochs of star formation in galaxies. Abundance distributions and gradients in galaxies were measured using PNe (Maciel & Koppen 1994;Magrini et al. 2016;Kwitter & Henry 2021). In the MW, negative radial oxygen abundance gradient for both thin and thick disc were constrained using PNe (Stanghellini & Haywood 2018), indicating inside-out disk formation. With this aim a large survey of M31 planetary nebulae was undertaken (Bhattacharya et al. 2019a). While it is not possible to determine the [α/Fe] abundance ratio in PNe for constraining the chemical evolution, different chemical enrichment timescales which result from distinct star-formation histories also leave imprints in the abundance distribution of other elements (see Nomoto et al. 2013, and references therein). Bhattacharya et al. (2022, hereafter Paper IV) measured distinct oxygen and argon abundance distributions for the thin and thicker discs in M31. They measured a flat or slightly positive oxygen and argon gradients for the older, thicker disc, and a negative metallicity gradient for the younger thin disc in M31. These results are consistent with a major merger with 1:5 mass ratio for M31. In the simulations of Hammer et al. (2018), a gasrich satellite is accreted on to M 31 with an orbit along the GSS. The merger then heats the pre-existing M31 disc, generating the observed thicker disc. The cold gas accreted through the "wet" merger would lead to a burst of star formation and the formation of a late, more centrally concentrated thin disc. Because of the merger-driven disc evolution, stars in the thin and thicker discs of M 31 would have formed at different epochs under different chemical conditions. Those in the younger thin disc would have formed out of the pre-enriched interstellar gas in M31, mixed with metal-poorer gas brought in by the satellite.
In this paper we investigate whether the different chemical enrichment of the thicker and thin discs in M31 provide "smoking gun" evidence that the secondary gas infall predicted by such a gas-rich merger took place. We will use the direct measurements of oxygen and argon abundances for the M 31 disc PNe over the 2 − 30 kpc radial range to constrain the chemical enrichment and the star formation efficiency in the thin and thicker disc. We further combine the PNe abundance measurements with those for the HII regions in the M31 disc already available in the literature. The oxygen and argon measurements from Paper IV are briefly presented in Section 2. The use of oxygen and argon abundances as chemical tracers in M 31 is discussed in Section 3, as well as the chemical enrichment timescales inferred for the two disc components. In Section 4 we present the constraints on the chemical evolution and formation history of M 31. We address the gas content of the merging satellite in the M 31 disc in Section 5, and conclude in Section 6. Spectral range covers from 3685Å in the blue to approx 9200 Å in the red, with a spectral resolution in the range 850 − 1500. The oxygen and argon direct abundance measurements, via the detection of the temperature sensitive line [O iii] at 4363 Å , with their errors, for a magnitude limited sample of PNe in the M31 disc, are described in Paper IV which includes the catalogue with the measured quantities (their errors). We also refer to Paper IV for a comprehensive description of the PNe magnitude limited sample and the abundance gradients. In this work we extend the analysis of the oxygen and argon abundances to constrain the chemical enrichment and the star formation timescales of the discs of M31.

Data sample, abundance measurements and gradients
To this aim, we implement the identification of the thin/thicker discs in M31 of Bhattacharya et al. (2019b), based on the ages and dynamical properties of their PN populations, and then derive their star formation histories and chemical evolution. In the M31 magnitude-limited disc PNe sample, there are 75 high-extinction, 2.5 Gyr and younger PNe associated with the more rapidly rotating thin disc, and 130 low extinction, 4.5 Gyrs and older PNe, associated with the thicker disc, which has a larger asymmetric drift. The high extinction PNe are found at smaller radii (R GC < 22 kpc) than the low extinction PNe, which instead cover the entire extent of the disc (out to R GC = 30 kpc).
We note that the structural properties of thin/thick discs in nearby spiral galaxiess show thin discs embedded in thicker discs, with the latter having longer scale lengths than the thin component (Yoachim & Dalcanton 2006). This is similar to what is found for the M31 PNe 2.5 Gyr and younger sample. Furthermore the age determined for the thicker disk in M31 falls within the age range, from 4 to 10 Gyr, determined from Lick indices for the thick discs of local spirals (Yoachim & Dalcanton 2008).

The PN sample in M31 as tracers of the ISM chemical properties
As a first step in using the PN oxygen over argon ratio vs. argon abundance, 12 + log(Ar/H), distribution to constrain the ISM conditions, we explore the dependency of the PNe log(O/Ar) ratio on i) their circumstellar dust properties, ii) initial stellar masses and iii) ages. According to Ventura et al. (2017) theoretical models, the surface oxygen in the stellar atmospheres maybe modified during the AGB stage in two ways: oxygen in the stellar atmosphere can be enriched by the Third Dredge Up (TDU) which may occur in 2 Gyr old stars, with masses between 1.5 − 3M , and metallicity range [−0.5 : 0.0]. In the case of the M31 disc low extinction PNe, which are 4.5 Gyr and older, with low mass M * < 1.25M (see Paper II), their measured oxygen abundance are free of alterations from PAGB evolution, and TDU in particular. oxygen can be depleted by hot bottom burning (HBB) which occurs in M * > 3.5M progenitors during the AGB phase. In Paper IV (see their Appendix D), we found that the contribution to the M31 disc PN sample from very young (<300 Myr old) massive (M > 3.5M ) stars is negligible, therefore there is no evidence for oxygen depletion in the current M31 disc high extinction PN sample.
In summary, the current sample of low extinction PNe in the M31 disc are too old and of low masses for TDU to occur, while there is no evidence for the presence of very young (300 Myr) massive (> 3.5M ) PNe in the current high extinction PNe sample in M31, which may be affected by oxygen depletion. We refer to Appendix A for more in depth discussion of oxygen and argon abundances in PNe. Differently from oxygen, argon is known to be invariant during the AGB evolution (Delgado-Inglada et al. 2014;García-Hernández et al. 2016;Ventura et al. 2017). Given the absence of oxygen modifications due to PAGB effects in the current M31 disc PNe sample, in what follow we then proceed to use the PN oxygen and argon abundances to study the chemical properties of the ISM at the time the PN stellar progeny were formed.

Oxygen and argon radial gradients for thin and thicker disc in M31
In Paper IV, we measured two distinct abundance distributions for the oxygen and argon, for the thin and thicker disc in M31.
The mean value of the oxygen abundance for the thin disc, < 12 + (O/H) > high−ext = 8.57 ± 0.03, is higher than that of the thicker disc, < 12 + (O/H) > low−ext. = 8.48 ± 0.02, although both distributions have large standard deviation values. Same trend for the mean values of the argon abundance: the argon abundance of the thin disc, < 12 + (Ar/H) > high−ext = 6.32 ± 0.03, is higher than that of the thicker disc, < 12 + (O/H) > low−ext = 6.25 ± 0.02. When the two abundance distributions are compared in pairs, the two-sample Anderson-Darling test rejects the null hypothesis that the two distributions of each element are drawn from the same, underlying, distribution. Hence the two discs in M31 are chemically distinct in oxygen abundance distribution and argon abundance distribution.
Regarding the abundance gradient with radius, in Paper IV we find a steeper negative radial gradient for the oxygen abundance for the thin disc, (∆(O/H)/∆R) high−ext = −0.013 ± 0.006 dex/kpc, which is consistent with that measured for the HII regions (Zurita & Bresolin 2012). Paper IV also measured a nearflat and slightly positive radial gradient for the oxygen abundance of the thicker disc, (∆(O/H)/∆R) low−ext = 0.06 ± 0.003 dex/kpc. The measured radial gradients for the argon abundance are (∆(Ar/H)/∆R) high−ext = −0.018 ± 0.006 dex/kpc and (∆(O/H)/∆R) high−ext = −0.05 ± 0.003 dex/kpc respectively. The results of Paper IV are consistent with results of previous studies (Sanders et al. 2012;Kwitter et al. 2012;Peña & Flores-Durán 2019), whose oxygen gradient measurements were dominated by the more numerous low-extinction PNe associated with the thicker disc. We refer to Paper IV for a more extensive comparison of these radial gradients with those of the MW and other spirals.

Radial gradients of the log(O/Ar) values
Oxygen and argon abundances are reliably measured in the M31 PNe (see Section 2.2, Appendix A and Paper IV). Since PNe evolve from parent stellar populations covering a range of ages, their log(O/Ar) ratios probe the ISM conditions at the differ- The bestfitting radial log(O/Ar) gradient is shown for HII regions (black), high-(blue) and low-extinction (red) PNe. The middle panel displays three independent linear fits for three radial ranges: the solid line is the linear fit to the entire data; for within 14 kpc and beyond 18 kpc, the linear fits are shown with dashed lines. ent epochs of their birth, and thus provide important constraints for the chemical evolution models of the ISM (see review by Nomoto et al. 2013 and work by Kobayashi et al. 2020a).
In M31 we identified three population of tracers in three different age ranges. The HII regions are tracing the chemistry of the ∼ 300 Myr young stellar population, while the high extinction PNe, which trace the thin disc of M 31, have ages about 2.5 Gyrs or younger. The low-extinction PNe that are associated with the thicker disc of M 31, are 4.5 Gyrs or older; see Paper II and Paper IV for further details. Figure 1 shows the galactocentric radial distribution of the log(O/Ar) values for HII regions 1 (upper panel), high-extinction PNe (middle panel) and low-extinction (lower panel) PN samples, in the R GC = 2-30 kpc radial range. Parameters of linear fits are also noted in Table 1. We find a negative radial gradient for the HII regions, with relatively large error bars. For the low-extinction PNe, we determine a slightly positive radial gradient, with a mixture of high and low log(O/Ar) values at any radius over the 2 − 20 kpc radial range of the disc. For the high-extinction PNe, the figure indicates distinct log(O/Ar) distributions in three radial regions. Within R GC < 14 kpc, < log(O/Ar) >= 2.28 ± 0.05 with a slightly negative radial gradient. In the R GC > 18 kpc outer region < log(O/Ar) >= 2.62 ± 0.28 with a negative radial gradient as for the inner region. In the intermediate radial range 14−18 kpc, which includes the most active star forming region in the M31 disc (Kang et al. 2009), the log(O/Ar) values have a wide spread, including the largest inferred values (> 2.5) for the PNe in our sample. We discuss the radial gradients and the constraints from the log(O/Ar) distributions of the different age ranges on the chemical evolu-   tion models for the thin and thicker disc in M31 in the following sections.

Oxygen vs argon as tracers of the enrichment history in the M 31 discs
The star-formation history of a galaxy leaves chemical imprints in its ISM through enrichment with different elements. The information of the ISM chemical conditions at the time of birth of a star is encoded in its element abundances. While both argon and oxygen are produced from core-collapse supernovae, argon is additionally produced by Type-Ia supernovae (Kobayashi et al. 2020a,b). Thus even though both argon and oxygen are α-elements, they do not have lockstep behavior and the star-formation histories of parent stellar populations will imprint information on the log(O/Ar) vs. argon abundance, 12 + log(Ar/H), distribution of stars with different ages.

M 31 PNe distribution of the in the log(O/Ar) vs. 12 log(Ar/H) plane
The left panel of Figure 2 shows the oxygen-to-argon abundance ratio, log(O/Ar), vs. 12+log(Ar/H) for all the M 31 PNe. In order to identify the general trend in the log(O/Ar) vs. argon abundance, we adopt the following procedure. We divide the high-and low-extinction PNe in bins of argon abundance, The number of PNe in each bin is lower for the high-extinction PNe as they are fewer in total 2 . We then calculate the mean log(O/Ar) values of PNe in each 12+log(Ar/H) bin after removing the 2 × σ outliers from the mean 3 . In Table 2 and Table 3 we provide the mean values of log(O/Ar) in the 12+log(Ar/H) bins, their standard deviation in bins, the error on the mean value and mean measurement error in each bin, for the low and high extinction PNe samples. We discuss them in turn.
In Table 2, we compare the standard deviation with the mean measurement error in each bin for the low extinction PNe: the scatter is of order √ 2× the mean measurement error for most bins. We thus infer that there is an intrinsic scatter of the measured log(O/Ar) for the low extinction PN sample, and this scatter is of a similar order of magnitude to the mean error in the bins. We thus adopt as the error for the log(O/Ar) average value, in each bin, the error on the mean log(O/Ar) value added in quadrature to the mean observation error of the PNe, in each bin.
For the high extinction PNe in Table 3, we notice that the standard deviation is smaller than the average measured error for about half of the bins. We still adopt as the error for the log(O/Ar) average value, in each bin, the error on the mean log(O/Ar) value added in quadrature to the mean observational error in each bin, keeping in mind that the dispersion of values is smaller for 12+log(Ar/H) range 6.21 -6.32, and for the larger The distribution in the right panel of Figure 2 is clearly a function of argon abundance, with the highest log(O/Ar) values ( 2.5) at the lowest measured argon abundance (12+log(Ar/H) 6.0) and the lowest mean log(O/Ar) measured value ( 2.08) for the most argon abundant (12+log(Ar/H) 7.5) bin, for the high extinction PNe. We also show the clipped outlier PNe individually. They are few in number and do not affect the general trend.

Tracing chemical enrichment in the old disc with PN abundances
We now focus on the distribution of, and trend in, the mean log(O/Ar) values vs. 12+log(Ar/H) for the low-extinction PNe, which trace the early evolution. In Figure 3, their values are compared with the evolution tracks for oxygen and argon for the MW thick and thin discs (Kobayashi et al. 2020a). According to these models, argon is made in both SN II and SN Ia, with a greater fraction in SN II than for Fe and a smaller fraction than for oxygen. Therefore the decrease in the log(O/Ar) vs. 12+log(Ar/H) plane would signal the increase in the SN Ia contribution, which is indeed shown by the models. In the solar neighborhood, 34% of Ar is produced from SN Ia by present, see Kobayashi et al. (2020a) Hayden et al. (2015), while the knee at [Ar/H] = −1.2 is caused by the dependency of argon production on the metallicity of core collapsed Type II SNe. We first compare the log(O/Ar) vs. 12+log(Ar/H) binned distribution for low-extinction PNe with the chemical evolution models. The O and Ar abundance distributions of the M31 lowextinction disc PNe show a possible slightly positive and null (within the errors) radial gradient, respectively, see Paper IV. In A&A proofs: manuscript no. main_03082022 turn, the log(O/Ar) distribution has a near constant average values and similar large scatter over the entire 2-30 kpc range (see Section 2.4, Fig. 1 and Table 1 ). The measured null radial abundance gradients for the low-extinction PNe and the large and constant scatter across the entire radial range supports the view that their distribution in the log(O/Ar) vs. 12+log(Ar/H) plane is the result of an effective radial mixing of different chemical enrichment tracks over the 2 − 30 kpc radial range. Significant radial mixing of the old stellar population over the entire PHAT area was also reported by Williams et al. (2017).
In Figure 3, the comparison of the log(O/Ar) vs. 12+log(Ar/H) binned distribution for low-extinction PNe with the Kobayashi et al. (2020a) models shows that the highest log(O/Ar) at low 12 + log(Ar/H) values for M 31 disc lowextinction PNe are representative of the stellar population that forms soon after α elements are produced in core collapse SNs, from short-lived massive stars. After Type Ia supernovae start producing additional argon relative to oxygen, the ISM is enriched and stars subsequently form with decreasing log(O/Ar) values at increasing argon abundance (12 + log(Ar/H)). The same chemical evolution models, following the aforementioned processes, also show the decreasing trends in the analogous [α/Fe] vs. [Fe/H] plot for stellar abundances, as seen for the MW (see e.g. Hayden et al. 2015).
In the MW (also recently in UGC 10738; Scott et al. 2021), it was found that the thick disc is more metal-poor and α-enriched compared to the thin disc, occupying distinct regions of the [α/Fe] vs. [Fe/H] plot (see e.g. Hayden et al. 2015;Kobayashi et al. 2020a). The explanation for the α enhancement in the MW thick disc is that it has had a faster chemical enrichment timescale than the MW thin disc, with high star-formation efficiency in the thick disc at early times, while the extended star formation in the latter generates higher iron abundance values than those in the MW thick disc. Figure 3 shows that correspondingly the MW thin disc also reaches higher Ar abundances than the thick disc. On the contrary, differently from the MW thick disc, the 4.5 Gyr and older thicker disc in M31 reaches higher argon abundances than the MW thick disc and as high as those reached by the extended star formation in the Solar neighborhood.

Fresh infall of gas: oxygen and argon abundances of HII regions in the M31 disc
A useful insight on the chemical evolution of the M31 discs is given by the properties of the younger generation of stars with respect to the older low-extinction PNe in the M31 thicker disc. We thus explore the log(O/Ar) vs. 12 + log(Ar/H) distribution ( Figure 4) measured for the HII regions in the M 31 thin disc (see Section 2.4). Similar to the M31 PNe, the HII regions follow a decreasing trend in log(O/Ar) with 12 + log(Ar/H) with a significant scatter. When we compare the HII region log(O/Ar) distribution with the MW thick/thin chemical enrichment models, we see that the HII region values do not cluster near the end of the track for the MW thin disc, at high [Ar/H] values. At a given argon abundance, their large scatter supports the hypothesis that their distribution in the log(O/Ar) vs. 12 + log(Ar/H) plane is due to differences in present day abundances. The range of log(O/Ar) values, between 2.15 to 2.6, at argon abundance 12+log(Ar/H) 6.0 indicates spatial variations of the abundances in the M31 gas phase of the thin disc, due to different degree of mixing of metal poor gas with the enriched ISM in the disc, and also with the starburst in the 10 − 17 kpc region (Kang et al. 2009). The spatial variations are plausibly linked with the recent production of oxygen in the starburst and rainfall of gas over the M31 disc, which may be related to the extra planar HI gas (Westmeier et al. 2008).

Constraints on the chemical enrichment of the thin disc in M31 with PN abundances
In this section we investigate the differences in the chemical enrichment of the high-extinction PNe in M 31, which are tracing the kinematically thin disc of M 31 (Paper II), with respect to the low-extinction older PNe. We note that high-extinction PNe occupy regions of the log(O/Ar) vs. 12+log(Ar/H) plane that, even if with some overlaps, are somewhat different from the lowextinction PNe, see Figure 2 right plot. We find that the decrease for the mean log(O/Ar) vs. 12+log(Ar/H) is linear for the lowextinction PNe. The mean log(O/Ar) vs. 12+log(Ar/H) for the high-extinction PNe are larger than those for the low extinction PNe in most bins, but for the 12+log(Ar/H) = 6.1 -6.3 range (3 bins), where they are lower.
We compare the linear fit to the mean log(O/Ar) values in the 12+log(Ar/H) = 6 -6.7 range, for the high-and low-extinction PNe in the M 31 disc, and validate them statistically. The highand low-extinction PNe can be fitted by linear functions with a Bayesian Information Criterion 4 (BIC; Schwarz 1978) values of -50.42 and -48.99, respectively. The linear fit to the lowextinction PNe (slope = −0.54 ± 0.09; intercept = 5.7 ± 0.58) is the best-fit, see Figure 5, with higher order functions not producing better fits (giving higher BIC values instead). A third order polynomial does however provide a better fit to the highextinction PNe (BIC=-56.17) as shown in Figure 5, thus validating the deviation from simple linear decrease.
If both the high-and low-extinction PNe were following the same enrichment history, the high-extinction PNe, which have an age of ∼ 2.5 Gyr or younger (Paper II), should populate low log(O/Ar)-higher argon abundances region of the evolutionary track, see Section 3.2 and 3.3. The existence of high-extinction younger PNe with high log(O/Ar) values at relative lower argon abundances support a secondary infall event with less chemically evolved gas. 4 It is a criterion for model selection among a finite set of models based on bayesian statistics. Overfitting may result from adding parameters to increase the likelihood of a function. BIC introduces a penalty term for the number of parameters in the model in order to avoid overfitting.

Constraints on chemical evolution and the formation history of the M 31 disc
As illustrated in Section 3, galactic chemical evolution models of the MW thick disc and solar neighbourhood (Kobayashi et al. 2020a) have a near linear decrease of log(Ar/O) values in the 12+log(Ar/H) = 6 -6.7 range (see Figure 3). Deviations from a simple linearly decreasing chemical evolution track may be related to either quenching or a secondary infall of gas onto a galaxy which then causes a modification of the ISM chemical abundances (see Section 4.2 for details and Matteucci 2021 for a review).

Regions of homogeneous chemical evolution in M31 discs
The chemical enrichment models do not include the effects of radial metallicity gradients. Therefore, as the high-extinction PNe have i) a steeper argon abundance gradient (Paper IV) and ii) have distinct behaviour of the log(O/Ar) values in three radial ranges, as described in Section 2.4, the steps prior to the chemical modeling include the identification of the disc regions with nearly flat log(O/Ar) radial gradients. In Figure 6 we show the distribution of the high-extinction PNe in the three radial ranges, inner disc (R GC < 14 kpc) [Left], starburst region (14 ≤ R GC ≤ 18 kpc) [Middle], and outer disc (R GC > 18 kpc) [Right], with respect to the low extinction PN values for the entire M31 disc. While in Section 3.4 the high extinction PNe have larger log(O/Ar) values than the low extinction PNe, we find that the high extinction, younger, PNe in the inner disc cluster towards the end, and below, the low extinction PNe distribution at relatively high [Ar/H], see left panel of Fig. 6. Considering their younger age, they are consistent with the chemical evolution only if also there is a dilution of the ISM by infall of metal poorer gas. In the middle panel of Figure 6, the high extinction PNe show a large range in argon abundances for a 2.5 Gyr and younger evolution. Those points on the top of the low-extinction PNe distribution may be consistent with the strong starburst located in this disc region (Kang et al. 2009) and a range of initial gas mixing. In the right panel of Figure 6, the high log(O/Ar) values at relatively low argon abundances for the younger high extinction PNe, with respect to the low extinction PNe distribution, indicate a higher star formation efficiency in those outer regions, i.e. shorter timescales in the thin disc with respect to the thicker disc, with the progenitor of these younger PNe being formed mostly out of the less enriched gas.

Imprint of a secondary gas infall on the oxygen and
argon abundances in the inner (R GC < 14 kpc) M31 thin disc from a chemical evolution model  Grisoni et al. (2017); Spitoni et al. (2019Spitoni et al. ( , 2021. According to these models, the MW thick disc stars form rapidly with high star-formation efficiency which leads to high [α/Fe] values for its stars. Eventually star-formation is suppressed as gas is depleted in the MW thick disc. After some time (t max ), further infall of gas, either pristine or pre-enriched, occurs which dilutes the available ISM and then triggers star-formation, at a lower efficiency. This may . Thin disc stars in the MW are formed initially at a lower metallicity which then increases with time. In the the MW thin disc, stars have lower [α/Fe] values than thick disc stars because of the diminished star formation efficiency. The time interval t max , the amount of accreted gas, and whether the latter is pristine or pre-enriched, drives the extent of these loops in the galactic chemical evolution tracks. Larger loops can be present for larger t max values (see Figure 14 in Spitoni et al. 2019) and/or larger amount of pristine gas. The loop also becomes narrower if the infalling gas is pre-enriched (Spitoni et al. 2021). As illustrated in Section 3 and 4.1, the oxygen and argon abundance distribution and ratios for young high extinction PNe support a secondary event of gas accretion. Hence we construct new, galactic chemical evolution models for the inner disc regions (R GC < 14 kpc) with two gas inflow phases for the M31 disc, using the same chemical evolution code by Kobayashi et al. (2020a). The models are constrained using the star-formation history and metallicity distribution function measured for the M 31 disc within the PHAT footprint (R GC 18 kpc) using isochrone fitting to the observed RGB stars (Williams et al. 2017). See Appendix B for a detailed illustration of these models.
In these chemical evolution models, the thicker disc older stars are formed by the first gas inflow, with the infall timescale of 1 Gyr. The SFR of the first burst differs in the models though, with a weaker first burst by a factor ≈ 2.5 for the second model, which is still within the uncertainties of the total star formation rate vs. age for the entire PHAT footprint, see Williams et al. (2017) and Appendix B. The star formation rate ceases after 2 Gyr and 3 Gyr, in the fiducial and weaker first burst model, respectively. The thin disc stars instead are formed by the second inflow, started at 9 Gyr with the infall timescale of 2 Gyr. In the two models, the SFR for the second burst differs at the 20% level, with the fiducial model having a relatively weaker second burst. The chemical composition of both gas inflows is set to be primordial, which results in the loop in the log(O/Ar) vs. 12+log(Ar/H) plane. The star formation timescales are 10 Gyr and 2 Gyr for the thicker and thin disc stars respectively, which means that the star formation efficiency is five times higher during the thin disc formation than for the thicker disc. A small amount of gas outflow (with the timescale of 5 Gyr) is also included. These new models are shown in Figure 7 together with the high extinction PN oxygen and argon measurements for the M31 ISM within 14 kpc, in the log(O/Ar) vs. 12+log(Ar/H) plane, and the low extinction PNe log(O/Ar) vs. 12+log(Ar/H) measurements over the entire disc. The list of the relevant parameters for the chemical evolution models is given in Table B.1.
We discuss the two chemical enrichment tracks in turn. In Figure 7 [Left], the fiducial model described in Appendix B reaches 12+log(Ar/H) ∼7.25 at lookback time ∼ 5 Gyr ago, following which the star formation is suppressed in the thick disc. Infall of primordial gas occurs ∼ 5 Gyr ago; it rapidly reduces the mean 12+log(Ar/H) value of the ISM and leads to a second burst of star-formation. At the start of the second burst of star-formation, prior to supernovae Type-Ia eruption, the log(O/Ar) value at corresponding lower metallicity increases, thereby producing the rising part of the loop (Figure 7-[Left]). The loop turns over once the supernovae Type-Ia eruption kicks in, leading to a decrease of the log(O/Ar) value with increasing 12+log(Ar/H), over the past ∼4 Gyr. The track reaches 12+log(Ar/H) values 7.2 at present times. The high argon abundance high extinction PNe are better reproduced by this fiducial model . In Figure 7 [Right], because of the weaker first star formation, the second model reaches 12+log(Ar/H) ∼6.7 at lookback time ∼ 5 Gyr ago. With the dilution and loop following a similar pattern as for the model in the left panel, the track reaches 12+log(Ar/H) values 6.7 at present times. This chemical evolution model of the ISM in the M 31 thin/thicker disc is in good agreement with the values traced by the low-extinction PNe and the loop covers the distribution of the high extinction PNe log(O/Ar) and 12+log(Ar/H) values ar R GC < 14 kpc. However, super-solar metallicity values traced by the high-extinction PNe are not reached in this model, and the model seems to underestimate the amount of low-metallicity gas added. Such remaining discrepancies may be addressed in a future investigation with more extended chemo-dynamical modelling.
In Figure 8 we show the chemical evolution enrichment track computed with reference to a burst of star formation from a gas infall with primordial composition and a timescale of 2 Gyr to reproduce the chemical enrichment of the thin outer disc, at R GC > 18 kpc. The lookback time shows the time from the start of the gas infall. The adopted timescale is even shorter than that of the MW thick disc, and it is as short as for nearby early-type galaxies It clearly illustrates the shorter star formation timescale and the lower argon abundance of the outer thin disc in M31.

Possible gas content of the merging satellite in the M 31 disc
From the age-velocity dispersion relation in the 14-17 kpc and 17-20 kpc, in Paper II we constrained the baryonic mass of the satellite that merged with the M31 disc. The estimated merger mass ratio from the age-velocity dispersion relation is 1 : 5 (Hopkins et al. 2008(Hopkins et al. , 2009, hence the baryonic mass of the satellite M sat = 1.4 × 10 10 M . Adopting the best-fit galaxy gasfraction to galaxy stellar mass relation (Díaz-García & Knapen 2020), the estimated gas fraction for the satellite would then be ∼19% , with a total mass in gas of M HI,sat = 2.8 × 10 9 M . We can compare the estimated gas mass brought in from the satellite with the i) mass in stars formed during the last 4.5 Gyr and ii) the amount of neutral HI still detected in the M31 disc. Williams et al. (2017) found that the total stellar mass produced over the last ∼ 2 Gyr in the M 31 disc is M stellar ∼ 4.5 × 10 9 M , and the present day total mass of HI gas in the M 31 disc is M HI = 4.23 × 10 9 M (Chemin et al. 2009) . Thus, prior to the burst of star-formation ∼2.5 Gyr ago (Bernard et al. 2015;Williams et al. 2017), a gas mass of M gas,2 Gyr 9×10 9 M must have been available in the M 31 disc. This amount is much larger than the estimated amount brought in by the merging satellite. Thus residual gas had to be still present in the pre-merger M 31 disc, which was then diluted by ∆[Ar/H] −0.5 dex as observed in the argon abundance distribution. The high extinction PNe were formed out of this diluted ISM within 14 kpc, as illustrated in Figure 7, and mostly from satellite gas at R GC > 18 kpc, see Figure 8 .
From the mass-metallicity relation (e.g. Zahid et al. 2017), the gas from the merging satellite galaxy would be pre-enriched, but at lower metallicity values than the ISM of the M31 premerger disc. A galaxy of mass equal to M sat is expected to have a metallicity ∼ −0.5 dex. The metallicity distribution of the RGB stars in the GSS from Conn et al. (2016), andCohen et al. (2018) provide evidence for a metal poor population of stars with [Fe/H] in the range [−1.0 : −0.2] dex. Then a second infall of gas related to the wet merger of the satellite onto the M31 disc would dilute the previously enriched ISM in the M31 pre-merger disc. The star burst triggered by the newly acquired gas would reproduce the lower [Ar/H] values and the log(O/Ar) distribution of the high-extinction PNe, as discussed in Section 4.2. The lower metallicity of the gas in the satellite, estimated either from the mass-metallicity relation or the observation of the GSS stars, is consistent with extent of the "loop" towards lower [Ar/H] abundances in the chemical evolution models shown in Figure 7.

Conclusions
We use the largest sample of PNe in the M 31 disc with oxygen and argon measurements from Paper IV, as well as archival measurements of these elements for HII regions. We compare their distributions with chemical evolution models for the MW thick disc, the Solar neighborhood and chemical enrichment models for the M31 discs with secondary infall phases. Our main con-clusions are that the log(O/Ar) vs. argon abundance distributions for the nebular emissions of the stellar evolution phases (e.g. HII regions and PNe) are valid alternatives to the [α/Fe] vs. [Fe/H] diagrams that have been widely used to constrain the star formation efficiency in the MW. By using this new approach, we determined a different chemical evolution history for the thin and thicker discs in the Andromeda galaxy, with the thin disc in M31 having been affected by a wet merger. In detail, our conclusions can be further illustrated as follows: -The oxygen-to-argon ratio, log(O/Ar), values of PNe are high at low argon abundances and monotonically decrease to lower values at higher argon abundances, consistent with predictions of the ISM properties from galactic chemical evolution models. This is the first vivid example that galactic chemical enrichment is imprinted in the log(O/Ar) vs. Ar plane. -The log(O/Ar) vs. argon abundance distribution for the lowextinction, 4.5 Gyr and older, PNe associated with the M 31 thicker disc supports a chemical evolution where the thicker disc forms from infalling gas with primordial composition at early times and reaches a high argon abundance similar to that of the Solar neighborhood in the MW (12+log(Ar/H) 6.7). -The log(O/Ar) vs. 12+log(Ar/H) distribution for the highextinction, 2.5 Gyr and younger, PNe associated with the thin disc of M31 overlaps with that of the low extinction PNe, with quantifiable differences. Three regions are identified in the M31 disc: inner disc, starburst region and outer disc. -In the inner disc (R GC < 14 kpc), the high-extinction, 2.5 Gyr and younger, PNe form after a secondary event with infall of metal poor gas which mixed with the pre-enriched ISM in the pre-merger M31 disc. This is in agreement with the twoinfall galactic chemical evolution model described in Section 4.2 and Appendix B. In the starburst region (14 ≤ R GC ≤ 18 kpc), the high-extinction, 2.5 Gyr and younger, PNe display a range of initial gas mixing and high log(O/Ar) values, consistent with the strong starburst located in this region of the M31 disc (see Kang et al. 2009). In the R GC > 18 outer disc, the high extinction PNe progeny are formed mostly in a burst of star formation 2 Gyr ago from satellite less chemically enriched gas. -The log(O/Ar) vs. 12+log(Ar/H) distribution for the HII regions in M31 has a large scatter at fixed argon abundance 12+log(Ar/H) 6.2. This spread is consistent with spatial variations and different degree of mixing, which is reasonably linked to current rainfall of metal poor gas from extra planar HI, part of which may have also originated from the satellite. -In contrast to the chemical enrichment history and structure of the thin disc in the MW, the thin disc in M 31 is less radially extended, formed stars more recently and at a higher star-formation efficiency, and had a faster chemical enrichment timescale than the thicker disc in M31.
The next steps of this investigation will explore the kinematical and chemical abundance properties of the PNe associated with the substructures in the outer disc and inner halo of M31.
A&A proofs: manuscript no. main_03082022 In Paper IV (see their Appendix D for a detailed discussion), using a larger sample of 101 MW PNe with abundance measurements and dust characterisation compiled by Ventura et al. (2017) compared to the chemical evolution tracks by Kobayashi et al. (2020a), we found that both CRDs and ORDs (as well as PNe with featureless dust) follow the MW stellar evolution tracks with no oxygen enrichment or depletion relative to argon abundances. The same can be seen in Figure A.1. However, as also noted Paper IV, many of the MW PNe with mixed chemistry dust (MCDs) that are metal-rich (12 + log(Ar/H) > 6.3) preferentially have log(O/Ar) values below the model tracks, indicating lower oxygen or oxygen depletion. These MCDs are the youngest (<300 Myr) most metal-rich PNe in the MW sample (García-Hernández & Górny 2014).
As such young PNe are expected to exhibit enhanced nitrogen abundances based on the AGB evolution models of their progenitors (García-Hernández & Górny 2014; Ventura et al. 2017), we can check whether any M 31 PN has enhanced nitrogen as further indication of a very young age and oxygen depletion, if any.
We measure nitrogen ionic abundances from the [N ii] 6548,6584 Å line fluxes for a sub-sample of 87 M 31 disc PNe which were observed in 2018 and 2019 with Hectospec at the MMT. The nitrogen elemental abundance is obtained by using the ICF from Delgado-Inglada et al. (2014) (the nitrogen abundances for the full sample will be presented in a future publication). The nitrogen abundances for 7 PNe in the M 31 disc are also measured by Kwitter et al. (2012) using Cloudy photoionisation models, independent of ICF. Figure A.2 [Top] shows that the nitrogen abundance measurements of these 7 PNe derived here are consistent with those of Kwitter et al. (2012), indicating that the Nitrogen abundance measurements are not dependent strongly on the adopted ICF. Figure A.2 [Bottom] shows the log(O/Ar) vs. 12 + log(Ar/H) distribution of these PNe coloured by their nitrogen-to-argon abundance ratio, log(N/Ar). Young PNe, possibly affected by HBB, would show higher values of log(N/Ar) than the general trend, owing to their expected enhanced nitrogen abundances Table B.1. Timescales of the two infall models for the M31 2 ≤ R GC ≤ 14 kpc region. t 1 t 2 t 3 τ i,1 = τ i,2 τ s,1 τ s,2 τ o,1 = τ o,2 fiducial model 3 9 12 5 1 1 4 weaker 1st SF 2 9 -5 5 1 3 outer disc ---1 1 -- Notes. This table summarised the timescales of the two infall models. The units are in Gyr. The timescales are: the timescales of the infall (τ i ), that of star formation (τ s ), and outflow (if there is, τ o ). The table lists three epochs that identify the truncation of the first infall (at t 1 ), the onset of the second infall (at t 2 ), and the truncation of the second infall (at t 3 ), also. See Appendix B for more details.
as well as low log(O/Ar) values due to oxygen depletion. To obtain the general trend of log(O/Ar) values as a function of 12+log(Ar/H), as in Section 3.1, we divide the PNe in bins of 12+log(Ar/H) such that there are 10 PNe in each bin (the bins with the smallest and largest argon abundance have fewer measurements, five and two respectively). We then calculate the mean and standard deviation of the log(O/Ar) values of PNe in each 12+log(Ar/H) bin, shown as diamonds in Figure A.2 [Bottom]. The PNe with the highest 25-percentile log(N/Ar) values have a mean offset of 0.02 dex (σ offset = 0.11 dex) from the general trend, while those with the lowest 75-percentile log(N/Ar) values have a mean offset of 0.01 dex (σ offset = 0.14 dex). Therefore, the PNe with higher log(N/Ar) values do not preferentially occupy lower log(O/Ar) values below the general trend. We can thus state that the number of PNe affected by HBB in our M 31 subsample is very small and does not drive the measured log(O/Ar) trends. This is consistent with the M 31 lowand high-extinction PNe having average ages ∼ 4.5 Gyr and ∼ 2.5 Gyr respectively with the bulk of the latter having likely formed in a burst of star formation ∼2 Gyr ago (Paper II), implying therefore a very small number of PNe with very young massive progenitors (affected by HBB).
To summarise, we find no conclusive evidence of AGB evolution effects with modification of the oxygen abundance in the nebula to be driving the trends in log(O/Ar) vs. 12 + log(Ar/H) for the M 31 disc PNe studied in this work. Any such effect is within the measurement errors. We thus conclude the oxygen abundance measurements for M 31 PNe reflect their birth ISM chemical abundances, within the errors.