1 Introduction

The evolution of the chemical properties of stellar populations and of the interstellar and intergalactic medium across the cosmic epochs provides unique information on the evolutionary processes driving the formation and evolution of galaxies. Theory and cosmological simulations give a relatively simple scenario on how dark matter evolve from the primeval perturbations, forming dark matter halos and large scale structures that accrete within their gravitational potential (e.g., Springel et al. 2018). However, the evolution of the baryonic component is much more complex as baryons interact with radiation and are subject to dissipative processes. The evolution of the baryonic component and how this results into the formation of stars and in the properties of galaxies as we see them in the local universe and across the cosmic epochs, has been subject to numerous models and cosmological simulations, which use different prescriptions and assumptions. The investigation of the evolution of the content of chemical elements provides tight constraints on such models. Indeed, the content of metals gives a measure of not only the integrated star formation in galaxies, but also on the fraction of metals lost through outflows and stripping. The metallicity, i.e., the content of metals relative to hydrogen and helium, is also sensitive to dilution resulting from inflow of pristine gas. Therefore, the investigation of the metal content and of the metallicity in galaxies provides truly crucial information on the key mechanisms involved in the evolution of galaxies. In addition, different chemical elements are enriched on different timescales by different populations of stars; therefore, the relative abundance of elements enables us to obtain unique constraints on the star formation history and on the late stages of the evolutions of single and multiple stars, stages which dominate the production of heavy elements.

The analysis of the metallicity and chemical abundances on spatially resolved scales (gradients) gives additional information on the processes that have regulated the growth and assembly of galaxies (e.g., inside-out or outside-in formation and/or quenching), as well as information on other internal phenomena such as galactic fountains, stellar migration and radial gas inflows.

Major advances of several observational techniques and the advent of new major observing facilities have recently enabled astronomers to probe the chemical evolution in the early cosmic epochs, by directly probing the early enrichment process of galaxies and of the intergalactic/circumgalactic medium, hence setting tight constraints on the models of early galaxy formation.

In this review we primarily provide an overview of the chemical properties of galaxies in the local universe and at high-redshift from an observational perspective, but by also discussing how such properties give important constraints on the galaxy evolutionary processes, including extensive comparisons with theoretical models and numerical simulations.

We will generally discuss the overall statistical properties of galaxies rather then focusing on individual galaxies, although we will unavoidably open some themes by quickly discussing the Milky Way (MW) or some nearby well-studied targets.

In the first part we also give a more technical overview of the methods, definitions and techniques adopted to measure the metallicity and chemical enrichment in stellar populations and in the interstellar medium, also discussing the strengths and weaknesses of the various methods. A large fraction of the review will be dedicated to the metallicity scaling relations, i.e., the trends of the stellar/gas metallicity with other galaxy properties, such as mass, star formation rate, gas content, environment, etc. We will then discuss how these properties on resolved scales, i.e. the investigation of metallicity gradients, which has recently been a steadily growing field. The relative chemical abundances of various elements is another major topic that, as mentioned above, provide key information, and to which we dedicate an entire major section; however, we cannot realistically review all chemical elements; hence we will mostly focus on some specific abundances ratios that are particularly useful to constrain the star formation history and galaxy evolution, and which have been measured across large samples of galaxies. In all of these major topics we discuss the finding both in the local universe and the evolution of these properties at high-redshift. We will dedicate a section to the current constraints of the metallicity in the host galaxies of Active Galactic Nuclei (AGN), while we only provide very brief discussions about resolved stellar populations and absorption systems associated with the intergalactic medium (IGM). Finally, we give an overview of the current understanding of the global metal content of galaxies and discuss their metal budget.

1.1 Expressing metallicity and chemical abundances

Different definitions are adopted to measure the abundance of metals and of the individual chemical elements. In contrast to other scientific disciplines (Agricola 1556), in astrophysics the term “metals” refers to the all elements heavier than helium. The “metallicity” Z indicates the mass of all metals relative to the total mass of baryons (dominated by hydrogen and helium):

$$\begin{aligned} {Z} \equiv {M}_{\mathrm{metals}}/{M}_{\mathrm{baryons}}. \end{aligned}$$
(1)

The relative abundance of two generic chemical elements X and Y is generally expressed in terms of relative number densities N, relative to the solar value, with the following notation:

$$\begin{aligned} {[X/Y]} \equiv \log {({N}_{{X}}/{N}_{{Y}})} - \log {({N}_{{X}}/{N}_{{Y}})_\odot }. \end{aligned}$$
(2)

When expressing the abundance of chemical elements relative to hydrogen, the following expression is also often used:

$$\begin{aligned} {{12}}+\log {({X}/{H})} \equiv 12+ \log {({N}_{{ X}}/{N}_{{H}})} , \end{aligned}$$
(3)

where the value 12 was introduced so that any element, even the most rare, has a Solar positive value of the expression 3.

Since oxygen is generally the most abundant heavy element in mass, often the “metallicity” is expressed in terms of oxygen abundance, generally under the implicit assumption that the abundance of all other chemical elements scale proportionally maintaining the solar abundance ratios. Therefore, often the metallicity is indicated as

$$\begin{aligned} \mathrm{12}+\log {(\hbox {O}/\hbox {H})} \equiv 12+ \log {({N}_{\mathrm{O}}/{N}_{\mathrm{H}})}. \end{aligned}$$
(4)

The Solar (photospheric) reference is still not completely settled. The review by Asplund et al. (2009) gives \({ 12} + \log {(\hbox {O}/\hbox {H})_\odot } = 8.69 \pm 0.05\) based on three-dimensional (3D), time-dependent hydrodynamical model of the solar atmosphere. However, for instance, Delahaye and Pinsonneault (2006) give \(\mathrm{12}+\log {(\hbox {O}/\hbox {H})_\odot } = 8.86 \pm 0.05\) based on helioseismology. However, it should be clear that in the context of this review the Solar metallicity (or Solar chemical abundances) does not have a particular meaning, at most being representative of the chemical enrichment of the (pre-)Solar Neighbourhood, i.e., a specific region of the Milky Way disc. The Solar metallicity/abundances should only be considered as reference values. What is important is that when comparing different studies they should be scaled to the same Solar reference value.

However, it should be clear that using the O/H abundance is only an approximation of the real metallicity of the gas, as the relative abundance of the chemical elements can vary in a drastic way with respect to the solar value.

1.2 The origin of the elements

While primordial nucleosynthesis accounts for the origin of hydrogen, deuterium, the majority helium and a small fraction of lithium (e.g., Cyburt et al. 2016), all other elements are produced by stellar nucleosynthesis or by the explosive burning and photodisintegration associated with the late stages of stellar evolution. Boron, beryllium and a small fraction of lithium are exceptions because they are produced by cosmic rays spallation of heavier elements.

Extensive reviews have been published on the origin of chemical elements through stellar processes (e.g., Rauscher and Patkós 2011; Matteucci 2012; Nomoto et al. 2013); therefore, in this section we only provide a quick overview of the basic processes associated with the production of heavy elements and a short summary of the primary sources of some key element, as well as the associated production timescales.

Stars on the main sequence burn hydrogen atoms, producing \(^{4}\hbox {He}\) atoms, either through the proton–proton nuclear reaction chain (pp-chain), or through the CNO-cycle; the relative role of these two processes depends on the stellar mass (the latter dominates at masses larger than about \(1.3~\hbox {M}_{\odot }\)) and, of course, metallicity.

At later stages, when hydrogen is exhausted in the core and this becomes hotter, helium burning starts, which produces \(^{12}\hbox {C}\) through the triple-alpha reaction, and \(^{16}\hbox {O}\) is also produced through the capture of an additional helium nucleus. During this phase hydrogen burning continues around the helium burning core.

If the initial stellar mass is less than \(8~\hbox {M}_{\odot }\) no additional burning phases take place.

If the stellar mass is less than about \(2~\hbox {M}_{\odot }\), the star undergoes the so-called core helium flash, a runaway process which results in an explosive expansion of the outer core. This results in a white dwarf surrounded by a planetary nebula.

Stars between 2 and \(8~\hbox {M}_{\odot }\) after leaving the main sequence enter the red giant phase and, once the helium burning core is exhausted, they go trough the so-called Asymptotic Giant Branch (AGB), characterized by both a hydrogen- and a helium-burning shell. Thermonuclear runaway of the latter results in a sequence of several He-shell flashes, with increased mass losses and strong stellar winds, which eject a large fraction of the previous burning products (in particular carbon and nitrogen). Such flashes are also responsible for producing large convection zones (which mix the product of nuclear reactions) and for the production of s-process nuclei. Eventually these stars also end their life as white dwarfs surrounded by planetary nebulae.

It should be noted that during the CNO-cycle, hence for stars more massive than about \(1.3~\hbox {M}_{\odot }\), while the primary reaction chain leaves unaffected the “catalytic” nuclei, secondary branches of the reaction produce \(^{14}\hbox {N}\) at expenses of \(^{16}\hbox {O}\) and \(^{12}\hbox {C}\). The result is a strong enrichment of \(^{14}\hbox {N}\) and also a drastic reduction of the \(^{12}\hbox {C}/^{13}\hbox {C}\) ratio. These effects are obviously strongly dependent on the stellar metallicity, which boosts the role of the CNO cycle and explains the “secondary” nature of nitrogen, i.e., whose production is greatly enhanced in metal-rich environments.

In stars more massive than \(8~\hbox {M}_{\odot }\), after helium burning, the core shrinks and increases the temperature to ignite carbon burning into \(^{20}\hbox {Ne}\) or \(^{23}\hbox {Na}\). Then heavier elements are subsequently burned as the lighter elements in the core are exhausted and, consequently, the core temperature increases. This process in particular results in the production and burning of elements such as oxygen, magnesium and silicon. During this sequential process occurring in the stellar core, the burning of lighter elements occurs in stratified shells around the core that, in the meantime, have reached the adequate temperature for igniting the associated nuclear reactions. The process ends when the core is primarily composed of iron and nickel. At this stage no more energy is gained by further nuclear reactions; the energy released can no longer sustain the hydrostatic equilibrium with the weight of the outer layer and collapse begins, yielding to a core-collapse supernova (type II, or Ib–Ic). The outward propagating shock wave results in the stellar material undergoing shock heating and explosive nucleosynthesis. The latter affects the abundance pattern of the outer layers, and especially the distribution of \(\alpha \) elements and the production of Fe-peak elements. Supernovae also produce elements heavier than iron through neutron rapid capture of iron-seed nuclei (r-process), followed by (slower) \(\beta \)-decay. The resulting yield of the various chemical elements implies complex calculations (e.g., Nomoto et al. 2013, and references therein) and, in particular, depends on the location of the “mass-cut”, i.e., the boundary between the core remnant that retain metals and the envelope that is expelled.

Type Ia supernovae are an additional important source of elements. SNIa are thermonuclear explosions of C–O white dwarfs that accrete mass either as a consequence of mass exchange in a close binary system with a non-degenerate star (“single degenerate scenario”) or as a consequence of the merging of two white dwarfs (“double degenerate scenario”) (see Maoz et al. 2014, for a review). SNIa explosions are thought to arise from the ignition of carbon burning in the C–O core, which results in the total disruption of the white dwarf (as the nuclear energy released by the dwarf is higher than its gravitational binding energy, Leibundgut 2001). The resulting nucleosynthesis produces elements primarily not only around the iron peak, but also silicon, argon, sulfur, and calcium. Clearly, since SNe Ia require first the formation a white dwarf from a low-mass (less than 8 M\(_\odot \)) star, and then a significant mass transfer from a companion star, the production of the SNIa is delayed with respect to the onset of star formation. Specifically, the first SNe requires a minimum timescale of 30 Myr (Greggio and Renzini 1983), although the bulk of SNIa explodes on longer timescales, due to the longer timescales associated with lower-mass stars and mass transfer, see Maoz and Mannucci (2012) and Fig. 1.

Finally, the merging of binary neutron stars is an additional source of elements beyond the iron peak, driven by the r-process.

Fig. 1
figure 1

Timescales of production of various elements after a single episode of star formation (a single stellar population, SSP) of solar metallicity, based on the model by Vincenzo et al. (in prep.), see text for details. The upper panel shows the production rate in M\(_\odot \) Gyr\(^{-1}\) normalized to 1 M\(_\odot \) of formed stars. The lower panel shows the cumulative mass produced, normalized to the amount after 1 Hubble time. Oxygen (red line) is mainly produced by CC SNe and, therefore, has the shortest formation timescales. Iron (blue line) is dominated by type Ia SNe, Carbon (black) has contributions from both kinds of SNe and from AGB stars. The production of Nitrogen (green) is dominated by AGB stars. In this plot, the production of elements before 30 Myr is due to CC SNe, type Ia SNe are described by a power-law \(t^{-1}\) after 40 Myr and up to the Hubble time, and AGB stars give additional contributions above this power-law at intermediate ages of \(\sim \) 0.04 to 5 Gyr

Clearly, each of these processes, and the associated release of chemical elements into the interstellar medium, is linked to the terminal stages of stars with a specific stellar mass on the main sequence and, therefore, with different timescales.

The amount of metals injected into the interstellar medium (ISM) by each star at the end of its lifetime is quantified by the so-called stellar yield, \(p_i\), which is defined as the fractional mass of newly formed ith chemical element, injected into the ISM, relative to the mass of the stellar progenitor on the main sequence. The computation of the stellar yields is quite complex and often subject to large uncertainties, as they also depend on the assumed mass loss and stellar rotation. The yields also depend on the progenitor metallicity, and in some cases the dependence is very strong (such as in the case of nitrogen, as discussed above). Compilations of stellar yields are given in Romano et al. (2010) and Nomoto et al. (2013). Here we only qualitatively summarize the main production channels of some of the key elements.

Figure 1 summarizes the production timescales of a few critical elements (O, C, N, and Fe) after a single episode of star formation, i.e., for a single population of stars (SSP) created together along a given initial mass function (IMF). Most of the \(\alpha \) elements (e.g., O, Ne, Mg) are thought to be produced by stars more massive than \(8~\hbox {M}_{\odot }\) and by the associated core-collapse SNe and, therefore, released into the ISM promptly, soon after the beginning of the star formation. Iron peak elements (e.g., Fe, Ni) are partly produced by massive stars, but most of them are produced by type Ia SNe and hence are released into the ISM with a delay ranging from \(\sim \) 40–50 Myr to a few Gyr, depending on the stellar initial mass function (IMF) and on the star formation history. It should be noted that zinc is often referred to as an iron-peak element, but it also seems to have an \(\alpha \)-like enrichment component and not always closely follow iron (e.g., Berg et al. 2016b). Elements such as carbon and nitrogen are partly produced by massive stars, but most of the production is due to intermediate mass stars (2 \(\hbox {M}_{\odot }<{M}_{\mathrm{star}}<8\,\hbox {M}_{\odot }\)), primarily through their AGN winds (or Wolf Rayet stars).], and, therefore, these elements are also subject to a delayed enrichment. The results in Fig. 1 are obtained by Vincenzo et al. (in prep.) for stars of solar metallicity, with the IMF from Kroupa et al. (1993), stellar lifetimes from Kobayashi (2004), stellar yields from Nomoto et al. (2013), SNe Ia yields from Iwamoto et al. (1999), and a \(t^{-1}\) delay-time distribution for the type Ia SNe (Maoz and Mannucci 2012). In Sect. 7, we will see that the different enrichment mechanisms and timescales of these different elements provide a powerful tool to constrain the evolution of the star formation history in galaxies.

2 Measuring metallicities of stellar populations

UV, optical and infrared spectra of galaxies contain a wealth of information about their stellar populations. Except for a number of galaxies in the local group where single stars can be resolved, galaxy spectra consist of the integrated light of the stellar population which is virtually always a composition of different generations. The broad-band colours and the continuum of the spectra are dominated by the distribution of stellar ages and metallicity, modulo dust reddening. The emission lines reflect the ionizing properties of the most-massive stars (and of the AGN, when present), the absorption features reveal the properties of stellar photospheres, stellar winds, and interstellar gas.

Several methods have been developed during the years to derive the chemical abundances from galaxy spectra. Two are the components that are usually considered, stars and the interstellar medium (ISM). This section deals with the former component, the next section with the latter one.

2.1 Rest-frame optical spectra

Extracting information about the metallicity of stellar populations is subject to an ongoing effort. A complete discussion of this topic is far beyond the scope of this review. Here we summarize the two main techniques that are often used to “invert” the spectra and derive the physical properties of the stellar population. Both are based on comparing observed with synthetic spectra.

The first tool to be developed was a set of standardized indices, calibrated on model spectra, aimed at maximizing the sensitivity to some parameters (e.g., age, or metallicity) while minimizing the dependence on the other parameters. These indices, including the “Lick indexes” (see, e.g., Fig. 2), were introduced in the 1970s and later refined and re-calibrated, and are still widely used, especially when only low-resolution spectra are available (e.g., Faber 1973; Worthey et al. 1994; Worthey and Ottaviani 1997; Trager et al. 1998, 2000b; Thomas et al. 2003, 2011; Schiavon 2007). Each index is defined either by a central bandpass whose flux is compared with those of two adjacent wavelength ranges intended to estimate a “pseudo-continuum”, or by the flux ratio in two nearby bandpasses to measure a break. The Lick indices use relatively large bandpasses, of the order of 50Å; this is useful to increase the S/N when needed but it is not optimal to derive individual abundances. Some indices, such as the depth of the 4000 Å break \(D_n(4000)\) (Hamilton 1985; Balogh et al. 1999) and the depth of the Balmer absorption lines, are mainly sensitive to age and to the fraction of young stars relative to old ones. Other indices are defined to be sensitive to abundances of particular elements, such as Fe and Mg in case of the [MgFe]\(^\prime \) and [Mg\(_2\)Fe] indices defined by Thomas et al. (2003) and Bruzual and Charlot (2003). The results often depend on the particular set on indices available or used, but these indices are at the base of most studies about unresolved stellar populations in nearby galaxies (see Sect. 5.1.1).

Fig. 2
figure 2

Image reproduced with permission from Onodera et al. (2015), copyright by AAS

The wavelength ranges covered by the set of Lick indices (grey rectangles) used by Onodera et al. (2015) are overplotted on a stacked spectrum of a sample of quenched galaxies at \(z\sim \,1.6\). The red line is the stacked spectrum, with orange lines showing the \(\pm \,1\sigma \) uncertainties. The green line is the best-fitting model stellar spectrum

More recently, spectro-phometric models of the stellar population have been developed aimed at reproducing the full observed spectra with a combination of input stellar populations with different properties (e.g., Bruzual and Charlot 2003; Conroy et al. 2009; Leitherer 2014; Chevallard and Charlot 2016; Wilkinson et al. 2017; Cappellari 2017, among many others), see Conroy (2013) for a recent review, and http://www.sedfitting.org by Budavari et al. for a complete and updated list of the available tools. These methods are in principle very powerful because they use all the information contained in the spectra, often also including broad-band photometry to simultaneously derive chemical abundance, age distribution (i.e., the star-formation history), dust extinction, and possibly other parameters such as the IMF. The weak point is that often the solution is not unique and there are strong degeneracies among the parameters. In particular, the problem is strongly non-linear in stellar mass and age, with the youngest and more massive stars often completely outshining the older, less massive stars that constitute most of the mass (see, e.g., Maraston et al. 2010). Uncertainties in the spectra of the input stellar populations, for example due to the presence of stellar rotation and binary stars (e.g., Levesque et al. 2012; Leitherer et al. 2014; Stanway et al. 2016; Choi et al. 2017) also affect the results. Despite these uncertainties, these methods are becoming the standard tool to analyse composite stellar population, in particular when spectra with enough spectral resolution and S/N are available.

2.2 UV spectra

At high-redshift (\(z>1\)), the UV part of the spectrum becomes more easily accessible at optical wavelengths and often constitutes the most important piece of information about the properties of stellar populations. At these wavelengths spectra are dominated by the photospheric properties of young, massive stars, allowing us to derive the properties of the on-going star formation activity.

A metallicity dependence of the depth of many UV absorption features is expected on the base of theoretical models (e.g., Leitherer and Heckman 1995; Leitherer et al. 2010; Eldridge and Stanway 2012) and was actually observed in the spectra of local starburst galaxies (e.g., Heckman et al. 1998, 2005; Leitherer et al. 2011).

In the UV, correlations between metallicity and equivalent widths (EWs) of many absorption features are a consequence of several physical causes. First, the contribution from the photospheres of the O and B stars can be seen, and these spectra depend on metallicity. Second, the strong, metal-dependent winds from hot stars dominate the spectra close to high-ionization lines, such as CIV and SiIV (e.g., Walborn et al. 1995). Third, the ISM produces the absorption of interstellar lines whose column density is related to metallicity. The interstellar lines are usually heavily saturated (e.g., Pettini et al. 2002b, see Savage and Sembach 1996 for a review) and the observed EWs are more sensitive to velocity dispersion than column density (González Delgado et al. 1998; Leitherer et al. 2011).

A complete list of features with their physical origin can be found in Table 1 of Leitherer et al. (2011). Here we only mention that the features that have usually attracted the highest interest are stellar wind lines such as NV\(\lambda \)1238,1242, SiIV\(\lambda \)1393,1402, CIV\(\lambda \)1548,1550, and HeII\(\lambda \)1640, interstellar lines such as SiII\(\lambda \)1190,1260,1304,1526, SiIII\(\lambda \)1206, OI\(\lambda \)1302, CII\(\lambda \)1334, SiIV\(\lambda \)1393,1402, FeII\(\lambda \)1608, AlIII\(\lambda \)1670,1854,1862, NiII\(\lambda \)1710,1741,1751, and MgII\(\lambda \)2796,2803, and photospheric lines such as FeV\(\lambda \)1360-1380, OV\(\lambda \)1371, SiIII\(\lambda \)1417, CIII\(\lambda \)1427, FeV\(\lambda \)1430, SiII\(\lambda \)1533, and CIII\(\lambda \)2300 (Leitherer et al. 2001; Pettini et al. 2002b; Leitherer 2014). Some of these are shown in the rest-frame UV composite spectrum of distant galaxies in Fig. 3. Weaker interstellar lines, which are not saturated, are sometimes also detected (Pettini et al. 2000, 2002b), specifically: FeII\(\lambda \)1144, SII\(\lambda \)1250,1259, NiII\(\lambda \)1317,1370,1703,1709,1741,1751, and SiII\(\lambda \)1808.

These features usually have complex dependences on age, metallicity and IMF. For these reasons, similar to the optical case, several authors have identified a number of indices optimized to depend only or mainly on metallicity (Heckman et al. 1998; Leitherer et al. 2001; Rix et al. 2004; Maraston et al. 2009; Leitherer et al. 2011; Sommariva et al. 2012; Zetterlund et al. 2015; Faisst et al. 2016; Byler et al. 2018). These works are based either on empirical spectral libraries or on theoretical stellar spectra. Both approaches have strengths and weaknesses. Empirical libraries automatically take into account the uncertain effects of stellar winds, but usually include only a narrow range of metallicities, limited by the metallicity distribution of young stars in the local universe, and are affected by the unknown contribution of interstellar lines (e.g., Kudritzki et al. 2016). In contrast, theoretical libraries can be built for a large number of different metallicities over a large range, but are affected by model uncertainties and do not take into account several effects present in the real spectra. The resulting indices depend on a number of assumptions, such as the age of the star-forming episode and the evolution of the SFR, usually assumed to be constant, instantaneous, or exponentially declining.

Fig. 3
figure 3

Image reproduced with permission from Steidel et al. (2016), copyright by AAS

Stacked rest-frame UV spectrum of star-forming galaxies at \(\langle z \rangle \sim 2.4\). The most important absorption and emission features are shown, colour-coded according to their nature: red, stellar absorption features; green, interstellar absorption features; blue, nebular emission lines; violet, fine structure lines

The IMF of the stellar population also affects the results. While this introduces yet another source of uncertainty, it also opens a way to test the dependence of the IMF on galaxy type, environment, star-formation activity and cosmic age, as discussed, for instance, by Lemasle et al. (2014), La Barbera et al. (2016), Fontanot et al. (2017), Sarzi et al. (2018) and De Masi et al. (2018).

3 Measuring metallicity and chemical abundances of the gaseous phase

In this section we provide a brief overview of the methods used to infer the metallicity and chemical abundances of the gaseous phase focusing on the interstellar medium (ISM), but also discussing the techniques exploited for the circum-galactic medium (CGM) and intergalactic medium (IGM), as well as for the intra-cluster medium (ICM). For each method we critically assess advantages and limitations. A more detailed, recent review can be found in Peimbert et al. (2017).

3.1 Direct method based on electron temperature

The spectra of ionized gas in astrophysical conditions are usually very rich of collisionally excited bf emission lines (CEL). The flux of each metal line is given by the abundance of the element (specifically the observed ionic species) times its emissivity. If the latter can be measured, then the abundance can be constrained accurately (Aller 1984).

The emissivity of these lines depends on both electron temperature \(T_{\text {e}}\) and on electron density \(n_{\text {e}}\). Once these two parameters are measured, ionic abundances follow from relations based only on atomic physics. For at two levels-ion, the rate of collisional de-excitation (per unit volume) of the transition 2 → 1 is given by n2n0C21, where n2 is the density of ions whose level 2 is populated, n0 is the density of colliding particles (typically electrons), and C21 is the collisional de-excitation coefficient given by

$$\begin{aligned} C_{21} = \left( \frac{2\pi }{kT_\text {e}}\right) ^{0.5} \frac{\hbar ^2}{m^{3/2}} \frac{\varOmega (1,2)}{\omega _2}, \end{aligned}$$
(5)

where \(\varOmega (1,2)\) is the “collision strength” of the transition, \(\omega _2\) is the statistical weight of the upper level 2, and there is a mild dependence on \(T_{\text {e}}\) as \(\sqrt{T_\text {e}}\).

The rate of collisional excitation is similarly given by n1 n0 C21, where

$$\begin{aligned} C_{12}=\frac{\omega _2}{\omega _1} \text {e}^{-E_{21}/kT_\text {e}}C_{21}, \end{aligned}$$
(6)

that depends exponentially on \(T_{\text {e}}\).

A CEL is produced by the collisional excitation of the upper level followed by a radiative de-excitation. Neglecting stimulated emission (usually not important in diffuse nebulae) and absorption, the population \(n_2\) of the upper level is given by

$$\begin{aligned} \frac{\text {d}n_2}{\text {d}t} = -n_2(A_{21}+n_{\text {e}}C_{21}) + n_1n_\text {e}C_{12}, \end{aligned}$$
(7)

where \(A_{21}\) is the Einstein coefficient.

At equilibrium (\(\text {d}N_2/\text {d}t=0\)):

$$\begin{aligned} \frac{N_2}{N_1} = \frac{N_\text {e} C_{12}}{A_{21}+n_\text {e} C_{21}} \end{aligned}$$
(8)

The critical density \(N_\text {c}\) is given by

$$\begin{aligned} N_\text {c} = A_{21}/C_{21}, \end{aligned}$$
(9)

and, therefore,

$$\begin{aligned} \frac{n_2}{n_1} = \frac{n_\text {e}/n_\text {c}\exp ^{-E/kT_\text {e}}}{1+n_\text {e}/n_\text {c}}. \end{aligned}$$
(10)

When the medium has a density much lower than the critical density of the transition (\(n_{\text {e}}\ll n_{\text {c}}\)), \(n_2\) depends exponentially on \(T_{\text {e}}\) and \(N_1\sim N_X\), where \(N_X\) is the density of the ion X. The volumetric emissivity \(J_{21}\) of a line is, therefore, given by

$$\begin{aligned} J_{21} = h\nu _{21}\frac{N_2A_{21}}{4\pi } \sim n_{\text {e}} N_x \text {e}^{-E_{21}/kT_\text {e}}. \end{aligned}$$
(11)

Once \(n_{\text {e}}\) and \(T_{\text {e}}\) are measured, ionic abundance can be obtained by comparing the flux of the CEL to the hydrogen recombination line. Adding up the abundances of the observed ions, and assuming an ionization correction for the unobserved ones, the total elemental abundance is derived (e.g., Aller 1954; Dinerstein 1990; Pilyugin and Thuan 2005; Pilyugin et al. 2006b, 2009, 2010a; Bresolin et al. 2009a; Perez-Martinez 2014; Pérez-Montero 2014; Pérez et al. 2016).

Electron density are usually measured by density-sensitive doublets, i.e., doublets that have critical densities not far from the gas densities, hence whose emission ratio depend strongly on the density around this regime, such as [OII]\(\lambda \)3726,3729 and [SII]\(\lambda \)6717,6731. Electron temperatures can be measured through “auroral” lines, i.e., lines coming from high quantum levels, whose excitation is very sensitive to temperature. The most commonly used of such lines are OIII]1661,1666, [OIII]\(\lambda \)4363, [OII]\(\lambda \)7320,7330, [SII]\(\lambda \)4069, [NII]\(\lambda \)5755, [SIII]\(\lambda \)6312 (e.g., Castellanos et al. 2002), whose fluxes is compared to other brighter lines from the same species but from a very different energy level. A more complete list of the line ratios used can be found in Pérez-Montero (2017). The most accurate results are only obtained when several of these lines are used, because they trace different regions of the emitting nebula, in particular regions of different ionization levels. The [SIII]\(\lambda \)6312 is an interesting line, although not often used. In fact, it can be observed to higher metallicity than, for example, [OIII]\(\lambda \)4363, but the corresponding nebular lines needed to measure the abundance are [SIII]\(\lambda \)9069,9532, which are in a part of the spectrum which is often not observed. Specialized routines are now available that derive temperature from line ratios, such as PYNEB (Luridiana et al. 2012, 2015), a PYTHON version of the STSDAS NEBULAR routines. A list of sources for atomic data can be found, for example, in García-Rojas et al. (2009) and Fang et al. (2015). A tutorial on the use on this method was recently published by Pérez-Montero (2017).

The practical use of this method is often limited by the intrinsic faintness of the auroral lines, which typically are 10–100 times fainter than the corresponding Balmer lines. While the auroral lines are routinely observed in local or low redshift, metal-poor, star-forming galaxies (Kennicutt et al. 2003; Izotov et al. 2006a, b, 2011, 2012, 2018a; Izotov and Thuan 2007; Kreckel et al. 2015; Haurberg et al. 2015; Amorín et al. 2015; Pilyugin et al. 2015; Lagos et al. 2016b; Ly et al. 2016; Sánchez Almeida et al. 2016), at higher redshift (\(z>1\)) only a few detection claims exist. While some of these measurements are based on UV lines such as OIII]\(\lambda \)1661,1666 (Villar-Martín et al. 2004; Erb et al. 2010), most of the claims are based on optical lines observed in the near-IR and are often affected by very low signal-to-noise ratios (often \(S/N\,<\) 2), are based on uncertain identifications (due, for example, to the nearby H\(\gamma \)), and would produce unrealistically high [OIII]\(\lambda \)4363/[OIII]\(\lambda \)5007 ratios (Yuan and Kewley 2009; Christensen et al. 2012; Stark et al. 2013; James et al. 2014b; Maseda et al. 2014; Sanders et al. 2016b). In contrast, auroral lines are often detected in stacked spectra of both local and distant galaxies (Andrews and Martini 2013; Trainor et al. 2016; Curti et al. 2017; Bian et al. 2018). Average relations between the temperatures of regions of different ionizations in HII regions and galaxies are also computed (Pérez-Montero and Díaz 2003; Izotov et al. 2006b; Pilyugin 2007; Pilyugin et al. 2006b, 2009, 2010a; Curti et al. 2017). These relations can be used when some auroral lines are not observed. Direct relations between “strong” lines and auroral line fluxes, the so-called ff relations, have been proposed by Pilyugin and Thuan (2005) and Pilyugin et al. (2006a, 2009). Recently, Curti et al. (2017) and Curti (2019) have extended these relations and applied them to SDSS galaxies, obtaining very tight relations between auroral and strong line fluxes.

A source of uncertainty of this method is due to the thermal and density structure of the emitting nebulae. The exponential dependence on \(T_{\text {e}}\) means that emission is dominated by the regions of higher \(T_{\text {e}}\) and the derived values of the parameters can be biased (e.g., Peimbert 1967; Stasińska 2002; Liu 2002, 2003; Esteban et al. 2004; García-Rojas and Esteban 2007; Peimbert et al. 2007). This problem is typically addressed by introducing a temperature-fluctuation parameter \(t^2\) that can be estimated by comparing different measurements of \(T_{\text {e}}\) (e.g., Peimbert 1967; Peimbert et al. 2004, 2012, 2017).

The problem of temperature fluctuations can be alleviated by exploiting coronal lines from various ionization stages of the same elements. These lines are emitted by different parts of the HII regions, including the partially ionized regions, and the thermal structure of the HII region can be better sampled (e.g., Campbell et al. 1986; Garnett 1992; Kennicutt et al. 2003; Bresolin et al. 2005; Pilyugin et al. 2009; Curti et al. 2017). This point is further discussed in Sect. 3.4.

3.2 Abundances from metal recombination lines

The permitted, recombination lines (RL) of metal ions are in principle the most direct way to derive the chemical abundances. In the usual conditions on HII nebulae (no stimulated emission), the volumetric emissivity of a permitted line of the ion X due to a transition between the levels i and j is

$$\begin{aligned} J_{ij} = \frac{h\nu _{ij}}{4\pi } N\left( X^{+i+1}\right) n_{\text {e}} \alpha ^{\mathrm{eff}}_{ij} \left( X^{+i},T_\text {e}\right) , \end{aligned}$$
(12)

where \(\alpha ^{\text {eff}}_{ij}\) has only a slow (about linear) dependence on \(T_{\text {e}}\). The ionic abundance is computed by comparison with hydrogen recombination lines, which have the same dependence on density, hence the estimated abundances are nearly insensitive to the gas density. The total element abundance is measured after assuming an ionization correction for the unobserved ions (Peimbert 2003; Tsamis et al. 2003; Esteban et al. 2004, 2014; López-Sánchez et al. 2007; Peimbert et al. 2007, 2014; Peimbert and Peimbert 2014; Toribio San Cipriano et al. 2017).

The RLs most commonly used to measure metallicity are OI8446,8447, OII4639,4642,4649, OIII3265, OIV4631, NII4237,4242, NIII4379, CII4267, CIII4647, and CIV4657.

The mild dependence on density and temperature reduces the impact of clumping and temperature fluctuation that can affect the CEL-\(T_\text {e}\) method described above. As the line emissivity in Eq. (12) is proportional to the abundance of each element, recombination lines from metallic species are very faint when compared to the H recombination lines, of the order of \(10^{-3}\)\(10^{-4}\) with respect to Balmer lines even for the most abundant elements like C and O. The detection of RL from metals is practically limited only to bright HII regions, planetary nebulae (PNe) and supernova remnants (SNRs), with spectra of high-resolution and high signal-to-noise ratio (e.g., Peimbert et al. 2004; García-Rojas and Esteban 2007).

3.3 Photoionization models

The widely adopted alternative to the direct method consists in using photoionization models to predict or interpret the relative strength of some of the main nebular lines to constrain the gas-phase metallicity. This approach has high potential, but also a number of limitations, as unavoidably only a small number of parameters involved in the photoionization calculations can be realistically explored, with simplified geometrical configuration, generally not properly reflecting the complexity and distribution of real HII regions. However, the advantage is that in principle there is no limit on the possible properties of star-forming regions that can be explored, especially in terms of metallicity range and properties of the ionizing spectrum. This flexibility enables also the exploration potential of systems at high-redshift, even in extremely metal poor environment or extreme ionizing spectra, which do not have local counterparts (e.g., Schaerer 2003; Kewley et al. 2013a; Jaskot and Ravindranath 2016; Xiao and Stanway 2018). The additional advantage is that such models can constrain, together with metallicity, other properties of the ionized gas, such as the ionization parameter. Finally, as discussed in the next sections, photoionization models enable the calibration strong-line diagnostic associated with lines that are not possible to calibrate empirically through the direct method because the required data are not available.

Generally, the classical approach is to use detailed photoionization codes, such as CLOUDY (Ferland et al. 2013) or Mappings (Binette 1985; Sutherland and Dopita 1993; Dopita et al. 2013) to generate a grid of models out of which a number of line ratios are extracted and proposed as diagnostic of the gas metallicity (e.g., Kewley and Dopita 2002; Nagao et al. 2011; Dopita et al. 2013, 2016; Jaskot and Ravindranath 2016; Gutkin et al. 2016; Chevallard and Charlot 2016; Feltre et al. 2016). Generally, with the exception of a few cases (Jaskot and Ravindranath 2016), models assume a simple plane-parallel geometry. Most of them assume that chemical abundances scale proportionally to solar, except generally for nitrogen whose abundance is assumed to scale with the global metallicity assuming a fixed relationship (see, e.g., the discussion in Nicholls et al. 2017). The effect of dust is generally included, both in terms of dust extinction and in terms of dust depletion of chemical elements, and the assumptions on dust distribution affect the resulting structure of the HII region (e.g., Stasińska and Szczerba 2001); dust depletion is generally inferred from Galactic studies and assuming that the dust-to-metal ratio remains constant with metallicity, which, however, may not be the case at low metallicities (De Cia et al. 2016; De Cia 2018). Another important assumption is that most ionized clouds are ionization bounded (i.e., the ionized zone is not truncated by the dimension of the cloud), but this assumption may not apply in number of galaxies, especially in some young, strongly star-forming systems (Nakajima and Ouchi 2014); some models have incorporated this possibility to investigate metallicity diagnostics, although restricted to the UV (Jaskot and Ravindranath 2016). The primary quantities that are varied in the grid of parameters are: metallicity and ionization parameter, defined as the dimensionless ratio of the incoming photon flux density and gas density at the cloud surface, normalized by the speed of light \(U=q/c=Q_{\text {ion}}/(4\pi r^2n_{\text {e}}c\)) where \(Q_{\text {ion}}\) is the number of ionizing photon emitted per unit time from the source, and r is the distance to the emitting cloud . Metallicity and ionization parameter are, for a given shape of the ionizing flux, the two parameters most important in affecting the flux ratios of the main nebular lines, and are often subject to degeneracies, in the sense that most emission line ratios depend on both parameters.

Some authors, especially in early models, adopted a constant gas density. Later Dopita et al. (2014) have pointed out that radiation pressure on dust is probably the primary physical mechanism regulating the physical properties of gaseous clouds and that, therefore, it is physically more sensible to adopt a constant pressure and derive the density distribution accordingly. In either case the assumed density or pressure can be another parameter that is varied to construct the grid of models.

While original models assumed only a single representative shape of the ionizing spectrum, more recent models typically explore a broad range of ionizing stellar continua from stellar population synthesis models, lately even expanding the range to include the contribution of stellar binarity (especially at low luminosity) and rotation (Levesque et al. 2012; Leitherer et al. 2014; Stanway et al. 2016; Choi et al. 2017; Xiao and Stanway 2018). The possible presence of turbulence introduces further uncertainties (Gray and Scannapieco 2017).

While such photoionization models are widely used, we warn about some important caveats and issues. Despite extensive efforts, these photoionization models are still quite simplistic, and are still far from capturing the complexity of the HII regions and their distribution in galaxies. An indication of this is that, while models tend to statistically reproduces fairly well many diagnostic diagrams, when individual systems are considered it is often very difficult to find a single model that simultaneously reproduce all observed nebular line ratios. The basic issue is likely to be that both within each HII region and among the several HII regions typically sampled by the large projected aperture of extragalactic surveys, the ionization parameter is not characterized by a single value, but by a broad distribution. The same issue applies to other parameters such as density, temperature and metallicity. Future generation models will hopefully incorporate this additional feature, although it is likely demanding to implement.

Another issue is that assuming a fixed relationship between N/H and O/H, typically of the steep form \(\hbox {N/H}\propto (\hbox {O/H})^2\) in the intermediate/high-metallicity regime, makes the flux of the nitrogen nebular lines (and in particular [NII]\(\lambda \)6584) hypersensitive to metallicity; however, one should take into account that nitrogen nebular lines are probing directly the abundance of nitrogen; if a system deviates from the assumed relation (which, for instance seems to be the case in galaxies of different masses and other sub-galactic regions), then the nitrogen emission line may provide deceiving information. A similar problem applies to carbon nebular lines when using UV spectra.

More recent models are introducing a Bayesian, or multi-features fit, approach in which the information coming from multiple nebular lines (and sometimes also from the continuum) is combined to identify the best model among those provided by the grid (e.g., Tremonti et al. 2004;Blanc et al. 2015; Chevallard and Charlot 2016; Pérez-Montero et al. 2016; Vale Asari et al. 2016; Pérez-Montero and Amorín 2017). These codes are certainly a step forward. However, some of them still make the assumption of a fixed relationship between N/H, C/H and O/H given a priori to the code. As mentioned, this risks the determination of the metallicity to be dominated by the flux of a single nitrogen or carbon line and, in addition, it is not really possible to distinguish between effect of differential chemical abundances and global metallicity. However, some of the most recent codes do consistently derive the global metal abundance and, separately, the abundance of nitrogen and carbon, without making assumptions on their abundance scaling with metallicity (Blanc et al. 2015; Pérez-Montero et al. 2016; Vale Asari et al. 2016; Pérez-Montero and Amorín 2017).

As will be discussed in Sect. 3.4, photoionization models tend to overestimate gas metallicity from \(\sim 0.2\) up to even \(\sim 0.6\) dex. The discrepancy is particularly strong at high metallicities. The origin of the discrepancy is not totally clear. Dust depletion is certainly a factor contributing for at least 0.1 dex; indeed, photoionization models do account for dust depletion, while the direct methods simply give the actual gas metallicity. Another potential problem is how nitrogen is included in the photoionization models; its assumed quadratic dependence on the metallicity may artificially boost the inferred metallicity. Alternatively, \({T}_{\mathrm{e}}\)-based metallicities may be biased low as a consequence of low-metallicity regions being characterized by brighter auroral lines. An interesting additional possibility is that the basic assumption that electrons in HII regions follow a Maxwell–Boltzmann distribution may be incorrect and that instead they follow a distribution similar to that observed in astrophysical plasmas in and beyond the solar system, i.e., a so-called \(\kappa \)-distribution (a generalized Lorenzian distribution), which is characterized by a more extended tail towards high energies than the Maxwell–Boltzmann distribution. This possibility has been investigated recently by Binette et al. (2012) and Nicholls et al. (2012) and Dopita et al. (2013), who have pointed out that, very interestingly, the introduction of a \(\kappa \)-distribution in models of HII regions solves some of the outstanding problems, such as the discrepancies in the temperatures measured by some auroral lines (see Sect. 3.4), and makes the photo-ionization derived metallicities in better agreement with those inferred from the auroral lines. This is certainly an area of study worth developing further. A revision of some of the nebular line diagnostics by including the effects of a \(\kappa \) distribution has already been presented by Dopita et al. (2013). It is also possible that the introduction of the \(\kappa \)-distribution is actually a way to reproduce the effect of gradients in the properties of the HII regions or of a collection of HII regions with different conditions, for example different temperatures.

3.4 Comparison among the different methods

In general the three methods described above give different results. The ratio of the abundances obtained from the more direct methods, i.e., from RLs and \(T_{\text {e}}\)-CELs, is often referred to as “Abundance Discrepancy” (AD, Tsamis et al. 2003; García-Rojas and Esteban 2007; García-Rojas et al. 2009, 2013; Peimbert et al. 2017; Toribio San Cipriano et al. 2017) and can be very high, up to a factor of \(\sim 5\) (Tsamis et al. 2003). This discrepancy is a serious limitation to our understanding of ionized nebulae and its origin is debated. Temperature fluctuations in HII regions are definitely present (García-Rojas and Esteban 2007) but their real effect on abundance determination is not clear and there also are indications, based on the comparison with fine-structure infrared lines, that they are not the source of the AD (Tsamis et al. 2003). The presence of small-scale, chemical inhomogeneities due to a clumpy, not well mixed gas distribution, possibly related to the presence of SN remnants, has been proposed by Tsamis and Péquignot (2005) and Stasińska et al. (2007) as a way to explain the different results from the RLs (dominated by the metal-rich clumps) and the CELs (dominated by the more diffuse medium), and by Binette et al. (2012) to explain the difference in temperature between [OIII] and [SIII]. Flourescence excitations via continuum or resonance lines have also been proposed as a way to provide unrealistic RL flux ratios and, therefore, affecting RL-derived metallicities (Grandi 1976; Liu et al. 2001; Tsamis et al. 2003).

Even higher differences exist between the results of CELs and RLs and those of the photoionization models; systematic discrepancies of 0.6–0.7 dex are not uncommon (e.g., Kewley and Ellison 2008; Moustakas et al. 2010; López-Sánchez et al. 2012), although some photoionization models tend to be in better agreement with the \(T_\text {e}\) method (Pérez-Montero 2014). Generally, direct \(T_{\text {e}}\) metallicities tend to be lower than those derived from photoionization models, with RLs usually in between these two values.

When making these comparisons one should also take into account that the direct \(T_{\text {e}}\) and RL methods estimate the metallicity of the gas phase only, while photoionization models generally take into account dust depletion and provide results in terms of total metallicity of the ISM (i.e., gas and dust). Therefore, \(T_{\text {e}}\) and RL metallicities should be increased by about 0.1 dex in order to have a fair comparison with the photoionization models.

The problem of understanding these differences and obtain the best estimates of the gas-phase metallicities can be addressed by comparing the values obtained from CELs, RLs, and photoionization methods with those of the young stars formed in the same region. The use of young stars is needed to guarantee that the stellar metallicity is similar to gas-phase metallicity, as young stars have been recently created from the same ISM that is observed in the HII regions, with the possible exception of some CNO reprocessing.

Independent estimates of metallicities can be obtained for blue supergiants (BSG) and red supergiants (RSG). Spectral type A and B BSG are massive (12–40 M\(_\odot \)), young (less than 10 Myr) stars in the short evolutionary stage leading to RSG. They are very bright and therefore their spectra can be obtained with the high-S/N needed to accurately measure chemical abundances. Their spectra show absorption lines from many elements (e.g., C, N, O, Mg, Al, S, Si, Ti and Fe). High-resolution spectra can provide abundances with uncertainties of less then 0.05 dex (e.g., Kudritzki et al. 2012), while lower resolution spectra can also be used to obtain uncertainties of \(\sim \) 0.1 dex (Kudritzki et al. 2008, 2016).

RSGs have initial masses between 8 and 35 M\(_\odot \), high-bolometric luminosities peaking in the near-IR, and young ages (less than 50 Myr). Abundances of RSG can also be reliably measured if good spectroscopy is available (Cunha et al. 2007; Davies et al. 2010, 2017a; Patrick et al. 2016).

The metallicity estimates derived for RSG and BSG are largely independent because the physical conditions of the atmosphere of the two types of stars are very different. The BSG method exploits hot atmospheres, and obtain metallicities from optical lines of singly and doubly ionized ions, while for RSG much colder atmospheres are used and near-IR lines from neutral metals are measured. The agreement between these two methods is excellent (e.g., Gazak et al. 2015; Davies et al. 2017a), and Zahid et al. (2017) recently showed that the derived values are also in excellent agreement with those derived by spectral fitting of stacked spectra of many galaxies.

The stellar values can be compared to those derived for the ISM with various methods. This comparison has been subject of a number of recent studies (Bresolin 2007; Kudritzki et al. 2012; Hosek et al. 2014; Lardo et al. 2015; Gazak et al. 2015; Bresolin et al. 2016; Davies et al. 2017a; Toribio San Cipriano et al. 2017) in which metallicity gradients of individual, local galaxies with different metallicities are derived with all the available methods to allow a comparison on local scales (see Fig. 4). In most cases there is a good agreement between stellar metallicity and the results of the \(T_{\text {e}}\)-CEL method for the ISM. The status is summarized by Bresolin et al. (2016) (see Fig. 5; Davies et al. 2017a). Photoionization models (Kewley and Dopita 2002; Tremonti et al. 2004; Maiolino et al. 2008) usually provide higher metallicities, with differences of the order of 0.2–0.3 dex.

Fig. 4
figure 4

Image reproduced with permission from Davies et al. (2015), copyright by ESO

The metallicity gradient in M300 measured by three different methods: BSG (blue triangles), RSG (red dots), HII with “\(T_{\text {e}}\)” metallicities. The agreement is remarkable

It is still not clear whether the stellar metallicities agree better with RLs or \(T_{\text {e}}\)-CEL derived values. According to the summary in Bresolin et al. (2016) (Fig. 5), a good agreement is found with the \(T_{\text {e}}\)-CEL at any metallicity, while the RL method seems to overestimate metallicity with respect to the stellar values for sub-solar metallicities. The situation is not clear yet and better agreement with RL metallicities are also found in the Orion nebula (Simón-Díaz 2010; Simón-Díaz and Stasińska 2011).

Fig. 5
figure 5

Image reproduced with permission from Bresolin et al. (2016), copyright by AAS

Difference between stellar and gaseous metallicity as inferred “directly” from the \(T_{\text {e}}\) method (circles) and from the recombination lines (orange squares) in a sample of local galaxies and star-forming regions, as a function of stellar metallicity. The vertical line shows the adopted Solar value

Summarizing, \(T_{\text {e}}\) direct gas-phase metallicities are in better agreement with the metallicities of young stars and therefore are now considered to be a more solid base for the “strong line” methods described in the next sectionFootnote 1 (Denicoló et al. 2002; Pettini and Pagel 2004; Andrews and Martini 2013; Curti et al. 2017).

3.5 Strong line calibrations

The emission lines required to directly measure the gas metallicity are very weak and challenging even with large telescopes. This has generally confined the use of the direct method to a few tens/hundreds local galaxies and HII regions, or resorting to the use of stacking of large number of spectra (e.g., Andrews and Martini 2013; Curti et al. 2017). At high-redshift the auroral lines needed to apply the \(T_\text {e}\) method have been detected (or marginally detected) only in a handful of sources. This issue has prompted astronomers to calibrate alternative diagnostics of the metallicity that exploit relatively strong nebular emission lines, which can be detected more easily, even in low S/N spectra. This technique is often referred to as the “strong line method”.

It is important to note that the strong line method is not a primary technique to derive metallicity but is a way to allow for an easier, albeit less precise, application of one of the primary methods introduced above.

These strong line calibrations have been performed empirically, through the direct methods (e.g., Pettini and Pagel 2004; Pilyugin and Thuan 2005; Pilyugin et al. 2010b; Pilyugin and Grebel 2016; Curti et al. 2017), through photoionization models (e.g., Zaritsky et al. 1994; McGaugh 1991; Kewley and Dopita 2002; Kobulnicky and Kewley 2004; Tremonti et al. 2004; Nagao et al. 2011; Dopita et al. 2016), or a combination of the two (e.g., Denicoló et al. 2002; Nagao et al. 2006a; Maiolino et al. 2008).

An important warning is that a number of the “strong line” metallicity diagnostics are highly degenerate with other parameters (ionization parameter, density, pressure,...) or even exploit indirect correlations between metallicity and other parameters, such as the correlation between metallicity and ionization parameter (Nagao et al. 2006a; Dopita et al. 2006) and the correlation between metallicity and nitrogen abundance. It is important to be aware of these issues and using multiple diagnostics is strongly encouraged. Within this context, for what concerns the direct calibrations, Pilyugin et al. (2010b) and Pilyugin and Grebel (2016) have developed methods which use strong line ratios to consistently provide metallicity, nitrogen abundance and mitigate the effects of ionization parameter.

In the context of photoionization models, the calibration of strong line ratios is rapidly making the way, as mentioned above, to codes which simultaneously fit multiple line ratios (Blanc et al. 2015; Chevallard and Charlot 2016; Pérez-Montero et al. 2016; Pérez-Montero and Amorín 2017), and, since most of these codes are publicly available, having calibrations for individual “strong line” diagnostics is becoming less compelling. However, the use of individual strong line diagnostics is still popular and handy, especially for high-redshift studies, where generally only a few nebular lines are measured.

The use of hybrid calibrations (i.e., both empirical, through the direct methods, and exploiting photoionization models) was in the past needed to properly cover the metallicity range (Nagao et al. 2006a; Maiolino et al. 2008). Indeed, while the low-metallicity range is fairly populated in terms of galaxies and HII regions with auroral lines detections for the \(T_\text {e}\) method, at high metallicities the weakness of the auroral lines has resulted in poor statistics and sparse sampling; this, together with concerns of biases of the auroral-line detected samples in the high-metallicity range, has led to efforts to complement the empirical calibrations (at low metallicities) with photoionization models-based calibrations at high metallicities. However, more recently extensive stacking of SDSS galaxies (Curti et al. 2017) has mitigated these issues and fully empirical (\(T_\text {e}\)-based) calibrations are available up to high- metallicities.

We finally point out that all these strong line calibrations, either empirical or through photoionization models, are based on HII regions and/or star-forming galaxies. This is an important caveat for various reasons. The presence of diffuse ionized gas (DIG, see Sect. 3.5.4) or of contamination by other excitation mechanisms, such as photoionization by shocks or harder sources, such as AGNs or evolved post-AGB stars, can strongly affect the nebular line ratios in a way that is independent of the metallicity and can vary from galaxy to galaxy. The selection of HII regions/star-forming galaxies is typically based on the so-called “Baldwin–Phillips–Terlevich” (BPT) diagnostic diagrams which attempt to isolate HII regions from regions excited by other mechanisms and sources (see Sect. 3.5.4). However, the simple demarcation reported by some authors (Kauffmann et al. 2003; Kewley et al. 2006b) is identified in an empirical way or by using considerations based on theoretical models. In reality the transition between different excitation mechanisms is certainly much smoother and mixed. It is therefore very likely that the datasets used for calibrations include contamination from AGN, shocks and post-AGB stars and, vice-versa, samples of star-forming regions are missed from the calibrations, especially in the high-metallicity regime. The same contamination issue is certainly present when applying these calibrations to galactic regions which are marginally resolved, or to the integrated spectra of galaxies in which different contributions are likely mixed. We will discuss these issues in the following.

3.5.1 Strong line calibrations: optical lines

Optical nebular line ratios are among the most widely exploited to constrain the metallicity in galaxies, both because some of the strongest nebular lines are in this wavelength range, and because this spectral range is easily accessible from ground with a variety of facilities and huge amount of data have been delivered by several surveys.

Table 1 provides a list of the main strong line diagnostic, and their definition, that have been proposed and calibrated either empirically or theoretically by several teams (McCall et al. 1985; Skillman 1989; McGaugh 1991; Zaritsky et al. 1994; Denicoló et al. 2002; Kewley and Dopita 2002; Kobulnicky and Kewley 2004; Tremonti et al. 2004; Pettini and Pagel 2004; Pilyugin and Thuan 2005; Nagao et al. 2006a, 2011; Maiolino et al. 2008; Pilyugin et al. 2010b; Pilyugin and Grebel 2016; Brown et al. 2016; Dopita et al. 2016; Curti et al. 2017). In some of these studies authors have suggested a combination of them to account for secondary (or primary!) dependences, such as the dependence on ionization parameter or nitrogen abundance (e.g., Kobulnicky and Kewley 2004; Pilyugin and Thuan 2005; Pilyugin et al. 2010b; Pilyugin and Grebel 2016; Curti et al. 2017).

In the third column we provide the references for some of the calibrations that have been proposed. We strongly suggest to use the empirical calibrations (\(T_\text {e}\)-based) obtained by Curti et al. (2017) and Curti (2019), and shown in Fig. 6, as they are based both on individual HII galaxies and stacks of SDSS galaxies in which also low-ionization coronal lines are detected and which, therefore, enable the empirical calibration to extend to high metallicities, while also mitigating potential biases.

As mentioned, given that each of these line ratios is generally degenerate with other parameters, it is advised to combine multiple of these diagrams to disentangle dependences (a publicly available routine to combine these calibrations can be found at http://www.arcetri.astro.it/metallicity/). However, if only some of these diagnostics are available it is important to be aware of a number of potential caveats, strengths and weaknesses, which are summarized in column 4 and 5 of Table 1, and discussed in the following.

R23 (\(\equiv \) log(([OII]\(\lambda \)3727 + [OIII]\(\lambda \)5007,4958)/H\(\beta \))) has been one of the parameters most widely adopted, as it involves emission lines of both the main ionization stages of oxygen, O\(^+\) and O\(^{2+}\), hence it is less affected by the ionization stage of the HII regions. However, it is subject to a significant dependence on the ionization parameter, which some authors attempt to correct by including additional transitions such as O32 (see below), as discussed below (Kobulnicky and Kewley 2004; Pilyugin and Thuan 2005; Nagao et al. 2006a). The additional problem of this diagnostic is that it is double-branched (i.e., a R23 value can be associated with two very different metallicities), hence identifying which of the two branches applies requires the use of another diagnostics. Also, it has a very weak dependence on metallicity at the switching between the two branches (roughly between \(12+\log (\mathrm{O/H})=8.0\) and 8.5), therefore loosing its importance on a significant range of metallicities. Finally, the use of [OII]\(\lambda \)3727 line (far in wavelength from H\(\beta \) and [OIII]\(\lambda \)5007) implies that this diagnostic is significantly sensitive to dust reddening, hence requiring correction.

R2 (\(\equiv \) log([OII]\(\lambda \)3727/H\(\beta \))) and R3 (\(\equiv \) log([OIII]\(\lambda \)5007/H\(\beta \))) are individually strongly dependent on the ionization parameter and on the ionization state of the gas, hence are not really meant to be used in isolation, unless each of them is really the only indicator available. Nevertheless, Fig. 6 shows that, at least for the stacked spectra, a tight, well-behaved relation exist between R3 and metallicity. This is mainly driven by the relation between metallicity and ionization parameter, but seems to work well across a large range of metallicity. R3 is also insensitive to dust reddening. In contrast, R2 is affected by the same issues as R23 in terms of being double-valued and being affected by dust reddening.

N2 (\(\equiv \) log([NII]\(\lambda \)6584/H\(\alpha \))) is very convenient to use, especially at high-z, as the two lines are very close-by hence requiring a very small spectral coverage and being completely free from dust reddening (unless differential dust extinction is present inside the galaxies). However, this diagnostic is actually primarily tracing the abundance of nitrogen, hence if galaxies follow a different N/O–O/H relation than the average sample used for the calibration then this indicator can be deceiving. This diagnostic is also well known to be strongly sensitive to ionization parameter (e.g., Nagao et al. 2006a). This ratio is also one of the parameters used to select star-forming galaxies in the BPT diagram, hence its calibration and application is highly sensitive to the detailed BPT-boundary adopted to select HII regions, implying that the resulting metallicity distribution of galaxies is also quite subject on the assumed BPT boundary. Finally, given its quadratic dependence on metallicity, N2 can be very low in low metallicity galaxies, down to \(-\,2\)therefore the, and the [NII]\(\lambda \)6584 line can be vary faint. This can introduce severe selection biases through the undetected galaxies if this effect is not properly taken into account.

S2 (\(\equiv \) log([SII]\(\lambda \)6717,31/H\(\alpha \))) is very similar to N2, but it has the advantage of not being linked to the nitrogen abundance. Sulfur is an \(\alpha \)-element like O, therefore S/O is expected to evolve much less than N/O.

O32 (\(\equiv \) log([OIII]\(\lambda \)5007/[OII]\(\lambda \)3727)) is mostly a tracer of the ionization parameter (Díaz et al. 2000; Kewley and Dopita 2002; Nagao et al. 2006a). The strong dependence on metallicity (Fig. 6) is mostly secondary and is mostly due to the correlation between ionization parameter and metallicity. Therefore, this diagnostic should really be used only if no other tracers are available and always bearing in mind that it is a very indirect tracer, through the UZ relation (which may evolve with redshift, see Kewley et al. 2013b), and subject to large scatter. Both O32 and R3 can be used to distinguish between the two branches of R23.

Ne3O2 (\(\equiv \) log([NeIII]\(\lambda \)3870/[OII]\(\lambda \)3727)), proposed by Nagao et al. (2006a) and Pérez-Montero et al. (2007), is equivalent to O32, i.e., is primarily a tracer of the ionization parameter, however it has the advantage that [NeIII]\(\lambda \)3870 and [OII]\(\lambda \)3727 are very close in wavelength, hence they require a very small wavelength range to be observed simultaneously (quite convenient at high-redshift), and their ratio is essentially unaffected by dust reddening. Also, it is the only ratio available in the near-IR bands above \(z\sim 3.6\) when H\(\beta \) and [OIII]\(\lambda \)5007 exit the near-IR K band. The drawback is that the [NeIII]\(\lambda \)3870 line is typically ten times fainter than [OIII]\(\lambda \)5007.

Fig. 6
figure 6

Strong line diagnostics as a function of oxygen abundance obtained by Curti et al. (2017) and Curti (2019). For the definition of the strong line diagnostics on the Y-axes please refer to Table 1. Green symbols are individual galaxies and circles are measurements from stacks of galaxies from the SDSS, colour coded by the number of galaxies in each stack (top colour bar). These calibrations are based on the most recent realization of \(T_{\text {e}}\) direct method

O3N2 (\(\equiv \) log([OIII]\(\lambda \)5007/H\(\beta \)) − log([NII]\(\lambda \)6584/H\(\alpha \))) is used quite extensively as it is a monotonic function of metallicity and unaffected by dust extinction (unless significant differential extinction is present). However, it is very sensitive to the ionization parameter and to the N/O abundance ratio, hence with the same major caveats discussed above for the other diagnostics affected by the same issues.

O3S2 (\(\equiv \) log([OIII]\(\lambda \)5007/H\(\beta \)+ [SII]\(\lambda \)6717,31/H\(\alpha \))) is a very promising diagnostics, proposed only recently by Curti (2019) and Kumari (in prep.). It is essentially equivalent to R23, but where [OII]\(\lambda \)3727/H\(\beta \) is replaced by [SII]\(\lambda \)6717,31/H\(\alpha \). Hence, as R23, it probes both the high- and low-ionization stages of the gas, but it has the additional advantage of being nearly insensitive to dust reddening. However, it shares some of the same issues as R23, in the sense that it has a secondary dependence on the ionization parameter and is double valued, although the inversion point is conveniently shifted to lower metallicities relative to R23.

N2S2H\(\alpha \) (\(\equiv \) N2/S2 + 0.264\(\cdot \)N2) was developed by Dopita et al. (2016) and expected to be less sensitive to the ionization parameter with respect to N2 and S2, and share with them the same requirement for the small wavelength range and the feature of being little affected by extinction. The calibration provided by (Dopita et al. 2016) through photoionization models is strongly dependent on the assumed N/O–O/H relation. This diagnostic has been re-calibrated empirically by Curti (2019), but the inclusion of the nitrogen line still preserves its dependence on the nitrogen abundance.

N2O2 (\(\equiv \) log([NII]\(\lambda \)6584/[OII]\(\lambda \)3727)) and N2S2 (\(\equiv \) log([NII]\(\lambda \)6584/[SII]\(\lambda \)6717,31)) are both very good tracers of the abundance of nitrogen relative to \(\alpha \) elements (e.g., Pérez-Montero and Contini 2009); however, they are often used also as tracers of metallicity thanks to the correlation of the nitrogen to \(\alpha \)-elements abundance ratio (N/\(\alpha \), such as N/O and N/S) with O/H at high metallicities. Yet, one should bear in mind that the dependence on metallicity is therefore indirect and that these diagnostics can be deceiving for systems deviating from the N/\(\alpha \)–O/H average relation. Both diagnostics fail to be a sensitive metallicity diagnostics at low metallicities where the N/\(\alpha \)–O/H flattens. N2S2 has the advantage, relative to N2O2, of being much less sensitive to dust reddening and requiring a small wavelength range.

Even more indicators, based a different combinations of several lines or on other elements such as Argon, have been proposed, among others, by Pérez-Montero and Díaz (2005), Stasińska (2006), Pérez-Montero et al. (2007), Kobulnicky et al. (1999), Dopita et al. (2016) and Pilyugin and Grebel (2016). Most importantly, some of these authors provides interlaced combination of multiple line diagnostics that attempt to explicitly take into account (and therefore mitigate) the effects of other physical parameters, such as the ionization parameter.

It is important to recall that, when attempting to investigate the nitrogen abundance N/O as a function of global metallicity O/H, obviously none of the strong line metallicity diagnostics involving the nitrogen lines should be used to measure the global metallicity O/H, as in this case both axes would essentially measure more or less directly the nitrogen abundance; for this kind of studies the global metallicity O/H should be measured by using only nitrogen-free diagnostics. Calibrations giving N/O (or N/S) in terms of N2O2 and N2S2 are given, for instance, by Pérez-Montero and Contini (2009) (based on empirical calibration) and Strom et al. (2017a) (through photoionization modelling).

Summarizing, the results of strong line methods based on optical lines critically depend on the method used to calibrate their line ratios, and calibrations based on the \(T_{\text {e}}\) method should be preferred. When lines of other non-alpha elements, such as Nitrogen, are used to measure Oxygen abundance, variations of the relative element abundances add another important source of uncertainty.

Table 1 List of optical strong line ratio metallicity diagnostics that have been proposed in the literature

3.5.2 Strong line calibrations: UV lines

Although extensive models and codes have been developed to infer the metallicity from the UV lines (Fosbury et al. 2003; Jaskot and Ravindranath 2016; Gutkin et al. 2016; Feltre et al. 2016; Chevallard and Charlot 2016; Pérez-Montero and Amorín 2017; Byler et al. 2018; Nakajima et al. 2018), and some works also measure metallicities directly using UV auroral lines (e.g., Erb et al. 2010; Berg et al. 2018), no clear, simple calibration has been proposed to derive the gas metallicity from ratios of UV strong lines. Ratios between CIII]1908, CIV1549 and HeII1640 (often some of the brightest line detected in the UV spectra of star-forming galaxies) are all strongly sensitive to the ionization parameter and shape of the ionizing continuum, and they are also produced/contributed by different physical processes: CIII]1908 is a collisional excited nebular line, CIV1549 is a resonance line which is generally blended with interstellar absorption and stellar P-Cygni profile, HeII 1640 is an interstellar recombination line in the highly ionized part of HII regions, but is also produced by the atmospheres of Wolf–Rayet stars. Pérez-Montero and Amorín (2017) also investigated the ratio of UV carbon lines with the optical hydrogen recombination lines; however, besides the problem of these ratios being extremely sensitive to dust reddening, no clear trend was found with metallicity, likely also because of the strong dependence on the ionization parameter.

However, the UV range contains the optimal nebular lines for the measurement of the C/O abundance ratio. Indeed, based on empirical calibrations, Pérez-Montero and Amorín (2017) have shown that the flux ratio

$$\begin{aligned} \mathrm{C3O3}\equiv \log {\left( \frac{\hbox {CIII}]1908+\hbox {CIV}1549}{\hbox {OIII}]1664}\right) }, \end{aligned}$$

is primarily sensitive only to the C/O ratio, and they provide an empirical relation in the form \(\log {(\hbox {C/O})} = -1.069+0.796~\hbox {C3O3}\). One however has to be aware of the caveats on the different production mechanisms of these lines discussed above, especially for what concerns CIV1549. If information on the electron temperature is available from auroral lines, then a more accurate determination of the C/O abundance (or at least of \(\hbox {C}^{2+}/\hbox {O}^{2+}\)) can be obtained from the CIII]1908/OIII]1664 ratio (Garnett et al. 1995, 1997), hence not having to rely on CIV1549 and the issues associated with this transition.

3.5.3 Strong line calibrations: far-infrared lines

The advent of IR space observatories such as Spitzer (Werner et al. 2004) and Herschel (Pilbratt et al. 2010), the prospect of new major spectroscopic IR surveys with the next generation satellites (SPICA, Roelfsema et al. 2018), and the possibility of observing the rest-frame far-IR wavelength range in distant galaxies with ground based sub-mm/mm observatories, such as ALMA and NOEMA, has fostered the investigation of mid/far-IR transitions as metallicity tracers.

Far-IR, fine-structure lines are expected to become a main coolant of HII regions at moderately high metallicities and low temperatures (e.g., Stasińska 2002). One problem of these transitions is that the most useful of them, within this context, are sparsely distributed across the broad IR wavelength range, often posing observational challenges. The additional problem is that these transitions tend to have low critical densities, making the gas density an additional important parameter affecting these diagnostics. Moreover, some of these transitions come from multiple gas components, making their modelling difficult and subject to large uncertainties and assumptions; for instance [CII]158\(\upmu \)m is emitted partly in the HII regions and in the Photo-Dissociation-Regions (PDR).

However, despite these caveats, IR transitions offer the possibility of investigating the metallicity of galaxies virtually without suffering from any dust extinction effects, therefore are worth exploring and using, whenever possible (Moorwood et al. 1980a, b; Lester et al. 1987; Rubin et al. 1988; Tsamis et al. 2003). For example, the metallicities of the central, obscured regions of starburst galaxies are only accessible via these far-IR lines, while metallicities derived from optical lines are likely to apply only to the outer, less dust-extincted part of these galaxies (e.g., Puglisi et al. 2017; Calabrò et al. 2018).

Most efforts have used photoionization modelling, and primarily using nebular lines coming only from HII regions, hence reducing model uncertainties and assumptions. After investigating various possible line ratios, Nagao et al. (2011) identified the line flux ratio ([OIII]\(\lambda \)52\(\upmu \)m+88\(\upmu \)m)/[NIII]57\(\upmu \)m as a good tracer of the gas metallicity, which is little dependent on gas density, ionization parameter and source of ionizing continuum. This has been confirmed later on by Pereira-Santaella et al. (2017) by adding a correction factor that minimizes the effects of some scaling relations not directly associated with metallicity (Fig. 7, left). Note however, that the metallicity sensitivity of this diagnostic is primarily due to the relation N/O–O/H assumed in the photoionization models adopted to calibrate it (hence the steep dependence at high metallicities and flattening at low metallicities). So one should be aware that this diagnostic is primarily tracing N/O.

Fig. 7
figure 7

Image reproduced with permission from Fernández-Ontiveros et al. (2017), copyright by ASA

Two metallicity calibrations based on far-IR lines. Left: (2.2\(\times \)[OIII]\(\lambda \)88\(\upmu \)m+[OIII]\(\lambda \)52\(\upmu \)m)/[NIII]57\(\upmu \)m flux ratio as a function of the gas metallicity, for different ionization parameters and sources of ionization, according to the photoionization models developed by Pereira-Santaella et al. (2017). Right: ([NeIIII]16\(\upmu \)m+[NeII]13\(\upmu \)m)/([SIV]11\(\upmu \)m+[SIII]16\(\upmu \)m) as a function of metallicity along with a range of photoionization models

In absence of one of the two [OIII] transitions, the [OIII]\(\lambda \)52\(\upmu \)m/[NIII]\(\lambda \)57\(\upmu \)m and [OIII]\(\lambda \)88\(\upmu \)m)/[NIII]\(\lambda \)57\(\upmu \)m ratios can still be used individually, but they are much more sensitive to the gas density. In fact, the line ratio [OIII]\(\lambda \)88\(\upmu \)m/52\(\upmu \)m in sensitive to density and can be used to measure it, together (Lester et al. 1983).

Fernández-Ontiveros et al. (2017) have also proposed the ratio ([NeIIII]16\(\upmu \)m+[NeII]13\(\upmu \)m)/([SIV]11\(\upmu \)m+[SIII]16\(\upmu \)m) as a good metallicity diagnostics, by computing both an empirical calibration and a range of photoionization models (Fig. 7, right), (although strongly affected by dust extinction) and Croxall et al. (2013) have proposed a calibration based on the ratios [OII]\(\lambda \)88\(\upmu \)m/H\(\alpha \) and [NeII]\(\lambda \)12.8\(\upmu \)m/[NeIII]\(\lambda \)15.5\(\upmu \)m.

Finally, Nagao et al. (2012) proposed that [NII]\(\lambda \)205\(\upmu \)m/[CII]\(\lambda \)158\(\upmu \)m could be a viable metallicity tracer for observations of high redshift galaxies with ALMA (at \(z>4\), where both lines are detectable from ground), again exploiting the N/O–O/H, although subject to significant model uncertainties.

3.5.4 Excitation source and BPT diagrams

Understanding the source of excitation of the nebular emission lines is important as the calibration of the metallicity diagnostics is typically restricted to a class of ionizing/excitation sources, generally young hot stars in star forming regions. This section gives a quick overview of some of the diagnostics used to classify galaxies and galactic regions, as well as of some of the main issues.

Diagrams comparing two or more emission line flux ratios contain a wealth of information about the conditions of the emission line nebulae. Different pieces of information can be obtained from different line ratios. Usually the influence of dust extinction is minimized either by using pairs of lines close in wavelength, or by plotting extinction-corrected line ratios.

The BPT diagrams, originally proposed by Baldwin et al. (1981), are among the most widely used. In these diagrams the [OIII]\(\lambda \)5007/H\(\beta \) ratio is used for its dependence on ionization parameter, and compared either to [NII]\(\lambda \)6584/H\(\alpha \) (N2-BPT), [SII]/H\(\alpha \) (S2-BPT), or [OI]/H\(\alpha \) (OI-BPT) (Veilleux and Osterbrock 1987). Diagnostic diagrams similar to the BPT plots have been discussed also considering lines at wavelengths other than optical, such as mid- and far-IR (Dopita et al. 2013) and UV (Byler et al. 2018).

In general, line flux ratios depend on many quantities such as total metallicity, chemical abundance ratios, ionization parameter and hardness of the ionizing continuum. For this reason BPT diagrams are widely used to compare observational and theoretical results on emission lines (e.g., Lehnert and Heckman 1996; Kewley et al. 2001a; Kewley and Dopita 2002; Dopita et al. 2006). These are the first diagrams that any photoionization model must reproduce, at least at the statistical level. Diagrams have also been introduced that are optimized to reveal the effect of single parameters. For example, the distribution of galaxies in the diagram O32 vs. R23 is sensitive to both the ionization parameter and metallicity, and can be used to break the degeneracies between these two quantities both in the local universe and at high redshift (Lilly et al. 2003; Hainline et al. 2009; Richard et al. 2011; Nakajima et al. 2013).

The presence of diffuse ionized gas (DIG) in local and distant galaxies (e.g., Oey et al. 2007) can also affect the position of galaxies on the BPT diagrams and bias the measurements of metallicity especially, in galaxies with low sSFR and low surface brightness of H\(\alpha \) emission (Sanders et al. 2017; Zhang et al. 2017; Lacerda et al. 2018).

Note that, since [NeIII]\(\lambda \)3870 and [OIII]\(\lambda \)5007 are emitted by essentially the same region and scale nearly proportionally to each other, the [NeIII]/[OII] and O32 have the same meaning and both depend mainly on ionization parameter. This is the reason why the [NeIII]/[OII] vs. O32 plots show small scatters (Steidel et al. 2016). The two lines can therefore be used equivalently.

Local galaxies show a typical bimodal distribution in the classic N2-BPT ([OIII]/H\(\beta \) vs [NII]/H\(\alpha \)), S2-BPT ([OIII]/H\(\beta \) vs [SII]/H\(\alpha \)) and O1-BPT ([OIII]/H\(\beta \) vs [OI]/H\(\alpha \)) diagrams, which are not sensitive to dust extinction. HII regions and galaxies whose emission is dominated by star formation, i.e., whose ionizing continuum is due to young, hot stars, follow a well defined sequence, where high values of [OIII]/H\(\beta \) correspond to low values in [NII]/H\(\alpha \), [SII]/H\(\alpha \) and [OI]/H\(\alpha \), as shown in Fig. 8 for the former two (Veilleux and Osterbrock 1987; Kewley et al. 2001a, b; Kauffmann et al. 2003; Bamford et al. 2008; Lamareille 2010). The position of the galaxies along the star-forming sequence is dominated by the luminosity-weighted ionization parameter. This is linked to metallicity (Nagao et al. 2006a; Dopita et al. 2006), therefore this is also a metallicity sequence, increasing toward higher values of [NII]/H\(\alpha \) (e.g., Andrews and Martini 2013; Curti et al. 2017).

Nebulae excited by a harder, AGN continuum are located in a different part of the diagram, with high values of all these diagnostic ratios: [NII]/H\(\alpha \), [SII]/H\(\alpha \), [OI]/H\(\alpha \) and [OIII]/H\(\beta \); therefore, the BPT diagrams are routinely used as a classification tool to separate the two populations of star-forming and AGN-dominated galaxies. The reason why AGN are located in these regions of the BPT diagrams is primarily the combination of two effects: high energy photons produced by AGN increase the fraction of highly ionized oxygen and also penetrate deeper into the gaseous clouds producing extended, partially ionized regions where singly ionized nitrogen and sulfur, as well as neutral oxygen, are abundant and where their collisionally excited transitions are the main coolant; moreover, the intense radiation field produced by AGNs also often results into high ionization parameter, which further boosts the [OIII]5007 emission.

Fig. 8
figure 8

Image reproduced with permission from Lamareille (2010), copyright by ESO

The three classical diagnostic diagrams N2-BPT, S2-BPT, and R3 vs R2 for SDSS galaxies. Galaxies are colour-coded according to the N2-BPT (blue and purple) and S2-BPT (green and cyan) diagrams. Blue are star-forming galaxies, purple are composite galaxies, green are Seyfert nuclei and cyan are LINERs/LIERs.

Composite galaxies with contributions from both star formation and AGN occupy an intermediated position.

The region of the BPT diagrams with high values of [NII]/H\(\alpha \), [SII]/H\(\alpha \)and [OI]/H\(\alpha \) , and relatively low values of [OIII]/H\(\beta \) are occupied by the so-called “Low Ionization Nuclear Emission Line Regions”, LINERS (Heckman 1980). The bimodality of the population in this region is even more evident in the S2-BPT diagram (Kewley et al. 2006a). Recent works have shown that this kind of emission is not confined to the nuclear region and is actually extended on kilo-parsec scales (Sarzi et al. 2010; Singh et al. 2013; Belfiore et al. 2016b), it has therefore been suggested to rename these regions as “LIER” (i.e., dropping the “N” that stands for “Nuclei” in the original acronym). It has been shown that LINER/LIER are a composite population that can be ionized by weak AGN, by post-AGB stars (associated with evolved stellar populations) or shocks (e.g., Shull and McKee 1979; Chevalier et al. 1980; Lehnert and Heckman 1996; Allen et al. 2008; Rich et al. 2010, Binette et al. 2012; Newman et al. 2014, see Sect. 3.5).

The most commonly used separation boundaries between these classes of galaxies are from Kewley et al. (2001a, b, 2006a) and Kauffmann et al. (2003), but other classifications have been proposed by Veilleux and Osterbrock (1987), Tresse et al. (1996), Dopita et al. (2000), Stasińska et al. (2006), Cid Fernandes et al. (2010) and Lamareille (2010). In all cases, galaxy classification based on different BPT diagrams (either N2, S2, or OI) can be very different (see Fig. 8). To obtain more robust and stable results Vogt et al. (2014) defined multiparameteric, 3D BPT diagrams, simultaneously using more than two line ratios.

Recently, a number of large size, integral field unit (IFU) spectrographs have been used to map a number of emission line galaxies. These data have been used to obtain spatially-resolved BPT diagrams studying the ionization properties of galactic sub-regions, both in local galaxies (e.g., Belfiore et al. 2015, 2016b; Davies et al. 2016, 2017b) and at high redshift (e.g., Newman et al. 2014; Curti 2018). In particular, the use of the VLT IFU spectrograph MUSE (Bacon et al. 2010) on nearby galaxies has allowed astronomers to obtain spatially resolved spectroscopy with very high spatial resolution on entire galaxies (e.g., Cresci et al. 2017; Venturi et al. 2017, 2018) distinguishing different ionisation conditions in different parts of galaxies.

The new, near-IR, multi-object spectrographs on 8-m class telescopes have allowed the observation of the near-IR spectra of a large number of intermediate redshift galaxies, obtaining high S/N ratios and reliable fluxes for the rest-frame optical lines. These studies have confirmed and clarified earlier results, based either on few galaxies or stacked spectra, that had shown the presence of a shift between the BPT distribution of high-z galaxies with respect to local forming galaxies (Shapley et al. 2005, 2015; Erb et al. 2006a, 2010; Kriek et al. 2007; Liu et al. 2008; Hainline et al. 2009; Finkelstein et al. 2010; Yabe et al. 2012, 2014, 2015a; Cullen et al. 2014; Newman et al. 2014; Masters et al. 2014; Hayashi et al. 2015; Salim et al. 2015; Sanders et al. 2016c; Steidel et al. 2016; Trainor et al. 2016; Strom et al. 2017a, b; Kashino et al. 2017, among others). High redshift galaxies appear to have higher [OIII]/H\(\beta \) for a given [NII]/H\(\alpha \), or higher [NII]/H\(\alpha \) for a given [OIII]/H\(\beta \), or both, see Fig. 9. Studies comparing the redshift evolution of only one line ratio with respect to other properties of the galaxies such as mass or sSFR also find evidences for evolution, but in this case the results are more ambiguous because it is not easy to distinguish an evolution of the ionizing properties of the gas from a simple metallicity evolution (Cullen et al. 2016; Holden et al. 2016)

Fig. 9
figure 9

Image reproduced with permission from Shapley et al. (2015), copyright by AAS

Three classical diagnostic diagrams for high-z galaxies. The N2-BPT, S2-BPT, and O32 vs. R23 diagrams are shown. Grey dots are the local SDSS galaxies, while green dots are star-forming galaxies at \(z\sim 2.3\) observed with MOSDEF at Keck. High-z galaxies show a clear offset from the local HII-sequence on the N2-BPT diagram, but not offset on S2-BPT diagram (although this claim has been debated, Strom et al. 2017b) nor on the local O32 vs R23 diagram

Many possible explanations have been proposed. Kewley et al. (2013a) ascribe this shift to an evolution of the ISM conditions in high-z starbursts, such as higher densities, harder ionizing radiation, or higher ionization parameter (see Fig. 10). Several authors (Brinchmann et al. 2008; Hayashi et al. 2015; Shimakawa et al. 2015; Kewley et al. 2015; Cullen et al. 2016; Kashino et al. 2017; Hirschmann et al. 2017; Kaasinen et al. 2018) explain the effect as the consequence of a high ionization parameter due either to higher SFR densities, or to a top-heavy IMF. Higher densities or higher ISM pressures are suggested by Shirazi et al. (2014), but this possibility is not supported by Strom et al. (2017a) who found no correlation between density and offset from the local BPT. Higher densities and higher ionization parameters are also seen in the rare local galaxies that fall into the same, offseted part of the N2-BPT diagram (Bian et al. 2016, 2017, 2018). These galaxies have low mass and high sSFR, but other SDSS galaxies with similar mass and sSFR do not necessarily have the same properties of the ISM, pointing towards other explanations for their peculiar properties. Also Poetrodjojo et al. (2018) do not find a relation between ionization parameter and sSFR in their sample of local, face-on spiral galaxies. Higher N/O abundance ratio has been proposed, as an explanation of the offset, by Pérez-Montero and Contini (2009), Masters et al. (2014, 2016), Shapley et al. (2015) and Yabe et al. (2015a), based on the evolution in the N2-BPT and not in the S2-BPT, on the absence of offset in the O32 vs. R23 diagrams (see Fig. 9), and on some N/O ratios measured at high redshift (Kojima et al. 2017). Jones et al. (2015a) reached the same conclusion on the base of a metallicity calibration at \(z=0.8\) based on the “direct” method. The difference in N/O could be related to the presence of hot (\(T_\text {e} \sim \) 80 000 K) Wolf–Rayet stars enriching the ISM with N, or to selection effects at high redshift (Masters et al. 2016). Cowie et al. (2016) proposed a combined effect of increasing N/O and increasing ionization parameters. However, a variable N/O has been questioned by Kashino et al. (2017) and Strom et al. (2017a), as these authors find evolution also in the S2-BPT diagram, and by the results of the models proposed by Hirschmann et al. (2017). Steidel et al. (2014), Nakajima et al. (2016) and Strom et al. (2017a) proposed the effect of harder stellar ionizing radiation as the dominant contribution to the BPT evolution, but this is not supported by the lack of correlation between the ionizing photon production efficiencies and the offset from the local BPT diagram found by Shivaei et al. (2018). In addition to the hardness of the spectra, Nakajima et al. (2016) suggested that matter-bounded HII regions could boost the value of the O32 ratio. Weak AGNs and shocks (Wright et al. 2010; Newman et al. 2014) could also have a role, albeit Steidel et al. (2014) found no evidence for high ionization lines from AGNs, and possibly selection effects could be present (Juneau et al. 2014).

Fig. 10
figure 10

Image adapted with permission from Kewley et al. (2013a), copyright by AAS

Illustration of the effect of the variation of several parameters on the offset from the HII sequence on the N2-BPT diagram. It shows how the SDSS star-forming galaxy sequence (in red) is modified by a harder ionizing field into the orange line, by an increase of electron density into the green line) by an increase of the ionization parameter q into the blue line) and by an increase of N/O into the purple line. The effects of most of these parameters are highly degenerate

In conclusions, it is not clear if there is one dominant effect and, actually, it is likely that many of them contribute to the observed BPT “evolution” with redshift. This large range of possible explanations ultimately derives from the intrinsic degeneracies of the photoionization models with respect to many parameters, as shown in Fig. 10.

Most importantly, in the context of this review, it is not clear how the strong-line calibrations derived for local galaxies are affected by this evolution of the excitation mechanisms/processes or of the relative abundances. Definitely the shift from the local BPT diagram results in an internal inconsistencies of the different methods (e.g., Kewley and Dopita 2002); the effect on metallicity can be severe when using N-based indicators (Newman et al. 2014; Cullen et al. 2016) or relatively minor (Salim et al. 2015), especially when using many line ratios (Brinchmann et al. 2008). The problem of evolution of metallicity diagnostics could be solved if direct detection of [OIII]\(\lambda \)4363 and other auroral lines in a significant number of high redshift galaxies could be obtained. Using observations at \(z\sim 0.8\), Jones et al. (2015a) obtained that the calibrations remain stable within 0.01 dex with respect to the local ones, revealing a small effect of the shift in the BPT diagram on metallicity determinations. Similarly, by exploiting the currently limited number of high redshift (\(z\sim 2\)) galaxies for which “direct” \(T_{\text {e}}\) metallicity measurements are available (mostly lensed galaxies), Patrício et al. (2018) have shown that diagnostics such as R\(_{23}\) and \(O_3\) appear to provide still reliable metallicity measurements, consistent with the local calibrations, while diagnostics involving nitrogen seem to be affected by a systemic offset, suggesting that a redshift evolution of the N/O abundance ratio, or selection effects, are probably affecting the reliability of these diagnostics. Diagnostics that are primarily tracing the ionization parameter and tracing metallicity only indirectly (or through a secondary dependence), such as O32 and Ne3O2 are those that show the largest dispersion and deviations in terms of ‘true’ metallicity with respect to that inferred from applying the local calibrations, confirming that enhanced ionization parameter in distant galaxies is heavily affecting these diagnostics. Certainly more direct \(T_{\text {e}}\) based measurements at high redshift are needed to further assess these issues.

Summarizing, the various flavours of the BPT and other diagnostic diagrams are critical tools to study the conditions of the ionized ISM in local and distant galaxies and determine the main contribution to the ionization, allowing for a classification of the galaxy. Clear signs of evolution are seen between local and high-redshift (\(z\sim 2\)) galaxies, and the origin of this evolution is debated. It is likely that the values of metallicity derived for distant galaxies with strong line methods are affected by this evolution, but the importance of this effect is still not clear.

3.6 Interstellar and intergalactic absorption lines

The ultraviolet spectral region is rich in resonant transitions of metal lines, associated with different kind of ions. If a clump of gas is observed against a source of radiation, the absorption introduced by such transitions provide a measure of the column density of the associated ionic species (through the curve of growth, see Savage and Sembach 1996, for a review) . More specifically the depth of the absorption produced by the transition between levels m and n of an ion i along the line of sight is given by \(\hbox {e}^{-\tau (\nu )}\), where the optical depth is defined as

$$\begin{aligned} \tau (\nu ) = \frac{\pi \hbox {e}^2}{{m}_{\mathrm{e}} {c}} {f} {N}^{{i}} \phi (\nu ,{b}), \end{aligned}$$

where \({N}^{{i}}\) is the column density of the ion i, f is the oscillator strength and \(\phi (\nu ,{b})\) describes the line profile as a function of frequency, which is also dependent on the Doppler parameter b.

If the column of hydrogen is also measured through one of the Lyman absorption lines, and constraints on the ionization status of the gas can be obtained, then accurate measurements of the cloud metallicity, or of its relative chemical abundances, can be inferred. Provided that the spectrum has high enough signal-to-noise, and adequate resolution, this techniques provides some of the most accurate measurements of gas metallicity and chemical abundances. Indeed, by measuring directly the column of metals/ions the technique is often (nearly) model independent and subject to little uncertainties.

Metal absorption lines have been extensively used to investigate the chemical abundances of the circumgalactic/intergalactic medium and in the interstellar medium of galaxies (e.g., Le Brun et al. 1997; Rauch 1998; Pettini et al. 2002a; Noterdaeme et al. 2008; Lehnert et al. 2009; Steidel et al. 2010; Rafelski et al. 2012; De Cia et al. 2018). In most cases quasars were used as background light, but sometimes also other bright sources such as Supernovae SNe and Gamma-ray Bursts (GRBs) constituted the beacon (e.g., Berger et al. 2012; Prochaska et al. 2007; Ledoux et al. 2009; Vreeswijk et al. 2014). The advent of sensitive spectrometers in space has enabled to extend these studies to local systems where stellar clusters or star-forming regions are used as background sources (e.g., James et al. 2014a; Tumlinson et al. 2011, 2017; Werk et al. 2014).

The fundamental difference with respect to luminosity-selected galaxies is that absorption lines, sensitive to column-density cross-section, probe the most numerous objects which are expected to be low-mass galaxies with low metallicity, a population whose emission is not easily accessible at high-redshift. Continuum or line emission metallicity studies instead tend to probe preferentially more massive (hence more metal rich) systems. Moreover, while metallicities measured through emission lines or stellar photospheric features are luminosity weighted (hence provide a view of the metal content in galaxies biased towards the most active regions), absorption systems provide a totally unbiased view of the metal content from this point of view.

Damped Lyman Alpha systems (DLAs) are by far the class of systems most widely exploited to trace metal enrichment and chemical abundances through absorption features at cosmological distances. These are defined as systems detected in absorption with column densities of neutral hydrogen \(\hbox {N(HI)}>2\times 10^{22}~\hbox {cm}^{-2}\). The several metal transitions observed in absorption together with Ly\(\alpha \) often enable a detailed characterization of their metallicity and abundance pattern. Good overviews of DLAs can be found in Pettini (2004) and Wolfe et al. (2005), although much work has obviously been done since these reviews.

An important aspect to bear in mind about DLAs and other absorption systems, especially in the context of this review, is that it is not really clear which kind of systems they are probing and what is their relation to galaxies. Even in those cases in which the optical counterpart of the DLA is identified, it is not clear whether the absorption system is tracing the outer parts of the galactic disc, or gas in the circumgalactic medium that has been ejected or is in an accretion phase, or a satellite galaxy/clump, or simply clumps in the nearby ICM/CGM (e.g., Fynbo et al. 2010; Krogager et al. 2012; Fumagalli et al. 2015), possibly characterized by a complex multiphase status and a poor chemical mixing (e.g., Zahedy et al. 2018). The fact that DLA are not directly tracing the bulk of the active and luminous part of galaxies has been discussed for instance by Pettini et al. (1999) and Pettini (2004), based on the weak redshift evolution of the metallicity in DLAs.

Prochaska et al. (2003) found evidence for a significant redshift evolution of the DLA metallicities, subsequently confirmed by various authors (e.g., Rafelski et al. 2012, 2014; De Cia et al. 2018; Poudel et al. 2018), however the extrapolation of these evolutions to \(z=0\) is around \(\log (Z/\hbox {Z}_{\odot })\sim 0.7\)–0.8, i.e., well below the metallicity of the bulk of local galaxies (at least within their effective radius). For the same reasons, and in particular due to the fact that absorption systems probe different regions than the bright parts of galaxies, it is difficult to compare the metallicity inferred from absorption systems with those inferred from emission lines or stellar features in galaxies. While early studies claimed a strong discrepancy (i.e. that absorption systems were systematically more metal poor than the metallicity inferred from emission line systems), these claims have been significantly revisited by studies that have compared absorption-derived metallicities with close-by HII regions and found that they are indeed consistent with each other (e.g., Bowen et al. 2005). Statistical studies have also found consistency between metallicities inferred from absorption systems and the luminous part of their counterpart once metallicity gradients are taken into account (e.g., Krogager et al. 2017).

Even with the difficulties in associating DLAs with the galactic component, absorption systems still provide the most accurate measurements of the relative abundances on an extremely wide range of chemical elements. By probing both refractive and volatile species absorption systems also provide some of the best determination of the relative depletion of chemical elements in dust grains in different environments and as a function of metallicity, as it will be discussed in Sect. 3.8.

DLAs will be recalled and discussed in several parts of this review. It is however important here to recall some caveats and difficulties of this technique, in particular the saturation of absorption lines, which in some cases only provides a lower limit on the column density of the ionic species, and the need to correct for the ionic fraction. The latter is often the main source of uncertainty as, unless multiple transitions from different ionization stages of the same element are observed, it may imply some modelling to account for the fraction of the unobserved ionic species. However, in many cases either several transitions are available to tightly constrain the ionization stage, or the column density is large enough that one can safely assume that the bulk of the cloud is self-shielded and therefore mostly neutral or in low ionization stages (Wolfe et al. 2005).

In the case of absorption lines used to probe the ISM in galaxies, the clump(s) covering factor is another potential source of uncertainty, but it can often be constrained either through the depth of saturated lines, or by using doublets whose relative equivalent widths is tied to the relative oscillator strengths.

The fact that transitions of different elements can be traced at different redshifts, depending on the available wavelength band, and at different column densities (depending on the oscillator strength of different transitions), may introduce difficulties in properly comparing the metallicities observed at different cosmic epochs and different systems (e.g., Rafelski et al. 2012, 2014; Wolfe et al. 2005, and reference therein).

At very high redshift (\(z>5\)) the Ly\(\alpha \) optical depth of the IGM is so high to make the continuum blueward of the Ly\(\alpha \) wavelength at the redshift of the background source nearly totally absorbed. Therefore, at these redshifts it generally becomes nearly impossible to obtain constraints on the column of hydrogen. However, at very high redshifts absorption systems still provide extremely precious information on the relative chemical abundances (modulo ionization corrections) (Becker et al. 2012; D’Odorico et al. 2013).

Summarizing, absorption lines probe a different galaxy population (or different galactic regions) with respect to that investigated by luminosity-selected samples, because the former are most sensitive to the column density distribution of gas. Due to their high column density, accurate measurements of the abundance of several elements can be obtained for the DLAs, allowing for a detailed study of the evolution of the abundance ratios and of dust depletion.

3.7 Chemical abundances from X-ray spectroscopy

X-ray spectroscopy has been extensively used to measure the chemical abundances in hot plasmas (\(T\sim 10^6\)\(10^8\) K), especially in the case of the hot ICM, (Mushotzky et al. 1996; Mushotzky and Loewenstein 1997; de Plaa et al. 2007; Sato et al. 2007; Mernier et al. 2016, 2018; Simionescu et al. 2018), but also in the galactic winds (e.g., Ranalli et al. 2008; Nardini et al. 2013; Veilleux et al. 2014).

Chemical abundances and metallicity are inferred through thermal collisional excitation and ionization models. This is relatively straightforward, as these plasmas are generally optically thin and in collisional ionization equilibrium, although accounting for the charge exchange X-ray emission has complicated the interpretation of some spectra (e.g., Liu et al. 2012; Zhang et al. 2014). The X-ray spectrum of hot plasma is rich of transitions from nearly all elements, from carbon to nickel, hence rich of information that can be used to infer chemical abundances. Of course, no transitions from hydrogen are detectable in the X-rays, hence the absolute metallicity can only be constrained by modelling the metal lines emission relative to the underlying free-free emission.

Until recently, many studies had mostly used CCD spectra resulting into very low spectral resolution that produces blending of many metal emission lines and, therefore, significant uncertainties in the resulting chemical abundances. The use of X-ray dispersion grating spectrometers has enabled astronomers to greatly improve the determination of the chemical abundances in plasmas by disentangling the transitions of several chemical elements (e.g., Ranalli et al. 2008; Mernier et al. 2016, 2018). However, the dispersive nature of this technique has often limited the use of this method to the systems with high surface brightness, such as galaxy cluster cores. The use of the innovative arrays of micro-calorimeters on board of the Hitomi space observatory have delivered a fantastic combination of very high spectral resolution and high sensitivity, therefore enabling an unprecedented, extremely detailed determination of the abundances in the Perseus cluster (Hitomi Collaboration 2017; Simionescu et al. 2018). Unfortunately, the limited lifetime of the Hitomi observatory has prevented the extension of such studies to larger samples.

Grating X-ray spectroscopy has been used also to detect highly ionized metal lines in Warm-Hot systems along the line of sight of bright background quasars (Buote et al. 2009; Zappacosta et al. 2010, 2012; Fang et al. 2010; Nicastro et al. 2018). These observations are very challenging, and do not really provide constraints on the gas metallicity (which is generally assumed or constrained from other tracers), however they can provide key information on the content of baryons in these intervening systems.

3.8 Dust depletion

In the interstellar medium a significant fraction of metals is locked into dust grains, and this is a major source of uncertainty for the determination of the metallicity and chemical abundances of the ISM. The “direct” methods only probe the gas phase metallicity, while photoionization models assume, a-priori a fixed depletion pattern of the various chemical elements, hence assume that the dust-to-metals ratio does not depend on metallicity or other environmental effects.

As mentioned in Sect. 3.6, absorption line studies are the most effective in tracing the depletion pattern, as they can both trace elements that are heavily depleted onto dust (e.g., iron) and elements of the same group that are little depleted.

Extensive studies and reviews on the dust depletion patterns and properties have been published (e.g., Savage and Sembach 1996; Draine 2003; Jenkins 2009; De Cia et al. 2016; Galliano et al. 2018) and an extensive discussion goes beyond this paper. However, in the following we provide some basic information that is particularly useful when investigating the metallicity of the ISM.

The amount of depletion is heavily variable and depends on environment (e.g., Jenkins 2014). Typically, of the order of 30–50% of the metals is locked in dust grains. The depletion fraction is quite different from element to element, depending on its “refractory” nature. Typically, in the ISM nearly all (\(\sim 90\)–99%) of iron is locked in dust, together with most of the silicon (\(\sim 30\)–97%, making silicate grains). About \(\sim 0\)–40% of oxygen and carbon are in grains, while other elements such as nitrogen, are almost totally in the gas phase (Jenkins 2009) even in the densest environments (e.g., Caselli et al. 2002). Zinc is often assumed to be free from dust depletion and therefore to be a good tracer of the abundance of the iron peak elements, but detailed studies have shown that it can suffers of significant depletion as well (Jenkins 2009; Berg et al. 2015b). The resulting structure of dust is complex, with cores, mantles and intrusions of the different elements. These are all important aspects when measuring the metallicities by using specific metal lines.

The assumption that the depletion pattern and, more generally, the dust-to-metal ratio is independent of metallicity seems to be reasonable only around solar metallicities, while recent studies (mostly based on the DLA observations) have shown that the dust-to-metal ratio decreases significantly with metallicity (e.g., about 50% lower at \({Z}\sim 0.1~{\hbox {Z}}_{\odot }\)) (Vladilo et al. 2011; Wiseman et al. 2017; De Cia et al. 2013, 2016), suggesting that grain-growth in the ISM is a dominant mechanism of dust formation. The variation of dust depletion with metallicity introduces a level of complexity that has not been incorporated yet in photoionization models.

4 The landscape of galaxy chemical evolution models

Each galaxy is subject to a number of processes acting together that determine its chemical evolution. Metal poor gas is accreted from the IGM; gas is used inside the galaxy to form stars, but stars can also be acquired via major or minor mergers; stellar evolution provides chemically enriched gas to the ISM via SN explosions and stellar winds; part of this enriched gas can leave the galaxy through galactic winds to enrich the CGM or the IGM; the gas inside galaxies is recycled several times and goes through a number of stellar generations; dust is created and destroyed, locking and releasing metals; the presence of a central AGN can affect the properties of the ISM and of the CGM, either by heating the gas or by removing it through AGN-driven winds, and both these effects can suppress star formation; dynamical interactions with nearby galaxies and with the ICM can alter the properties of the ISM and affect the star formation activity. All these effects can depend on several parameters such as dark matter halo mass, cosmic time, and environment. The study of the chemical enrichment of galaxies can give information on all these effects.

Different kinds of models have been developed to take into account all the effects discussed above in a consistent way and inside a cosmological framework. Extensive reviews of the methods used to develop galaxy formation models and of the current status of the subject can be found in Silk and Mamon (2012), Silk et al. (2014), Somerville and Davé (2015), Naab and Ostriker (2017) and Dayal and Ferrara (2018). In the context of galaxy formation models, Matteucci (2012) provides an extensive overview on the theoretical approaches to model and reproduce chemical enrichment of galaxies..

In this section we give a brief summary of the main classes of existing models of galaxy formation and evolution, which will serve as a reference for the discussions presented together with the observational results. As emphasized by Somerville and Davé (2015), all these models are based on a common set of basic physical processes and obtain similar, although not at all identical, results (see, e.g., Mitchell et al. 2018).

One critical problem that all models must solve is how to keep star formation efficiency low, as in all galaxies only a small fraction of baryons are converted into stars. Galaxies contain a smaller fraction of baryons with respect to the cosmic average, and this fraction increases with mass up to halo mass of the order of \(10^{12}\) M\(_\odot \) (e.g., Baldry et al. 2008; Papastergis et al. 2012). This requires the existence of mechanisms to remove baryons from galaxies, or prevent accretion. Star formation must also be limited. For this, the gas can be put into some state unsuitable for star formation rather than removed from the galaxy, or further accretion of gas can be prevented therefore depriving the galaxy of the necessary ingredient for a sustained and prolonged star formation. Sources of these kinds of feedback have been proposed to be, among others, SNe, radiation from young stars, AGN, and ram-pressure stripping due to a intracluster medium.

All models contain a few main, critical parameters that are described via mathematical formulae in the analytical and SAM models or with physical prescriptions in the numerical ones:

  • The star formation “efficiency”, defined as the number of stars formed per unit time per unit gas mass

    $$\begin{aligned} \epsilon = \hbox {SFR}/{M}_{\mathrm{gas}} , \end{aligned}$$
    (13)

    often referred to as the inverse of the gas “depletion time”, i.e., the time required for star formation to completely use up all gas if there is no further gas accretion and the star formation rate remains constant. Note, however, that in case of no accretion the SFR cannot remain constant, must decline exponentially, as simply obtained by solving the “closed-box” differential equation

    $$\begin{aligned} {{\mathrm{SFR}}} = -\hbox {d}M_{\mathrm{gas}}/\hbox {d}t=\epsilon {M}_{\mathrm{gas}}, \end{aligned}$$
    (14)

    hence the term “depletion time” is not totally appropriate, as it is actually linked to the e-folding time of the SFR. The star formation efficiency is often taken as constant, independent of the gas mass or surface density; this is equivalent of assuming that the Schmidt–Kennicutt relation is linear, in contrast to the original formulation, which has a slope of 1.4 (see Kennicutt and Evans 2012). Whether the slope of the relation is linear or superlinear is still a debated issue, which goes beyond the scope of this review. More detailed models tend to distinguish between the star formation efficiency associated with the sole molecular gas, i.e., the phase out of which stars form, and the global star formation efficiency, i.e., including the atomic component of the cold gas, which is typically inactive and mostly a reservoir on larger scales.

  • The outflow mass loading factor, i.e., the ratio between mass outflow rate and star formation rate,

    $$\begin{aligned} \eta = \frac{\dot{{M}}_{\mathrm{outfl}}}{\hbox {SFR}}. \end{aligned}$$
    (15)

    For star-forming galaxies, in which the outflow is primarily driven by SNe and radiation pressure from stars, the loading factor is typically observed to be around unity (e.g., Steidel et al. 2010; Heckman et al. 2015; Fluetsch et al. 2018), but it is expected to anti-correlate with the mass of the galaxy (as a consequence of the deeper gravitational potential well, Somerville et al. 2012), and it can be boosted by a large factor by the presence of an AGN (Cicone et al. 2014; Fluetsch et al. 2018). Often models leave the outflow loading factor as a free parameter that is constrained by fitting the observational data. When talking about outflow, it is often important to discriminate between gas that does not escape the galaxy (falling back, hence being recycled, on relatively short timescales) and gas that leaves the halo (or which is potentially re-accreted only very long timescale, close to the Hubble time). Some models introduce the concept of “effective” outflow loading factor by considering only the fraction of gas that is expelled and not recycled (Peng et al. 2015). Most models assume that the metallicity of the outflowing gas is the same as the average metallicity of the ISM in the galaxy (or of the galactic sub-region), while some models also consider the possibility of differential outflows in which the species produced by core collapse SNe are expelled more efficiently than other elements.

  • The gas inflow rate. This is one of the most critical parameters, as it regulates most of the galaxy evolution. Many models assume that the gas inflow rate is proportional to the SFR in the galaxy as this enables a simpler analytical solution of the differential equations that describe galaxy evolution and chemical enrichment. However, although mathematically convenient, we warn that this assumption is unphysical, as there is no physical reason why the inflowing gas should know about the SFR in the galaxy. Moreover, assuming that the accretion rate is proportional to the SFR, when combined with Eqs. (13)–(15), automatically implies that the gas mass and the SFR must decline exponentially with time, hence implying that the galaxy cannot follow the observed tight relation between mass and SFR dubbed Main Sequence of Star Formation (MSSF, e.g., Brinchmann et al. 2004; Daddi et al. 2007). The only case in which the inflow rate can be physically assumed to be proportional to the SFR is in the condition of perfect equilibrium in which the inflow rate is exactly compensated by the SFR and outflow rate. More physically sound models assume the inflow rate proportional to the halo mass. Generally the inflowing gas is assumed to be chemically pristine; this assumption may not be appropriate in several cases; indeed, as we will see in Sect. 5.1.6 and Sect. 8. the circumgalactic gas can be significantly enriched, also at high redshift.

It is important to emphasize that a meaningful comparison between the results of the models, especially those based on numerical techniques, and the observation must take into account all the selection effects on both sides. In fact, generally large populations of observed galaxies are compared to a large number of model galaxies derived from the simulated evolution of baryons inside a distribution of dark matter halos, and the selection of the objects in the models can be different to those affecting the actual observations. While the former are usually selected on stellar mass or halo mass, the latter are selected based on luminosity or gas column density. In other words, simulations must be “observed” using the same selections used for the telescope observations, producing “synthetic observations” of the simulation outputs. This is not always done, but the result of the comparison and the values of the derived physical quantities can depend critically on these effects (e.g., Governato et al. 2009; Scannapieco et al. 2010; Guidi et al. 2016, 2018).

4.1 Numerical simulations

These models start from N-body simulations of the effect of gravity on dark and baryonic matter and the hierarchical growing of structures via a continuous process of merging and accretion. Several implementations exist that can be based either on particles, or on a geometric mesh, or hybrid. Critical parameters of these simulations are the spatial and mass resolution, i.e., the finest detail that can be resolved and reproduced, and the total size of the simulation, i.e., the number of particles used and the size of the universe volume that is simulated. Being limited by the total computational time available, usually high resolutions corresponds to small volume sizes, and vice-versa. Often several simulations are made with the same physics and different choices of resolution/size (e.g., McAlpine et al. 2016).

In hydrodynamical simulations baryonic physics is included by solving the hydrodynamical equations albeit in some simplified form (e.g., Finlator and Davé 2008). Complex processes on small scales such as star formation and the effect of SN explosions on the ISM would require sub-parsec resolution, which is never achieved, and sub-grid recipes are still used. The impressive advances in computational power of the last decade, and a better understanding of the description of the basic physical processes have allowed these codes to reach the level where reproducing most of the observational constraints is possible. Several large simulations are now available and have released their databases to the public to allow for a better and more extensive analysis of the results. Among the best know, ILLUSTRIS (Vogelsberger et al. 2014; Genel et al. 2014; Genel 2016), ILLUSTRIS-TNG (Springel et al. 2018; Pillepich et al. 2018; Torrey et al. 2017), EAGLE (Schaye et al. 2015; Crain et al. 2015; McAlpine et al. 2016; De Rossi et al. 2017), HORIZON-AGN (Dubois et al. 2014, 2016), FIRE (Hopkins et al. 2014; Ma et al. 2016), and MUFASA (Davé et al. 2016, 2017). The review by Somerville and Davé (2015) provides an extensive description of these models. Some hydro codes are optimized to reach the best possible resolution, albeit on single galaxies (e.g., Pallottini et al. 2014, 2017; Katz et al. 2015; Costa et al. 2015).

4.2 Semi-analytic models (SAMs)

These models start from considering the merging trees, obtained from N-body dark matter simulations or Monte-Carlo techniques, that give rise to the formation of cosmic structures. The resolution is given by the mass of dark matter particles used. Galaxies are associated with dark matter halos, and evolve inside them by describing the important physics with analytic formulae, usually based on a number of free parameters whose values are varied to reproduce the observations, in a forward-modelling approach. As mentioned above, crucial quantities and processes are, among others, gas cooling time, star formation efficiency, feedback from SNe, feedback from AGN, morphological transformation, and metals distribution (e.g., Kauffmann et al. 1993; Lacey and Cole 1993; Cole et al. 2000; Kauffmann et al. 2004; Croton et al. 2006; De Lucia and Blaizot 2007; Benson 2012; Collacchioni et al. 2018).

The advantage of SAM models is that they are less computationally expensive and allow for a faster exploration of the effects of changing the description of the physical processes. The role of these processes can therefore be better understood. The weak points are that the physics is controlled “by hand”, i.e., the assumptions on the physical processes are approximate and not necessary realistic (which is however an issue also for many simulations), the evolution of baryons and dark matter may not be self-consistent because baryons evolve inside pre-determined dark matter halos, and a large freedom in the choice of the parameters is allowed so that it is not guaranteed that a unique solution can be obtained. Examples of SAM models, discussing in particular the chemical properties of galaxies, can be found in Thomas et al. (1999), De Lucia et al. (2004, 2017), Somerville et al. (2008), De Lucia (2010), Hirschmann et al. (2013), Fu et al. (2013), Yates and Kauffmann (2014), Porter et al. (2014), Cousin et al. (2016) and Zoldan et al. (2017).

4.3 Analytical models

For the MW and the galaxies in the Local Group with resolved stellar populations, a set of detailed chemical evolution models have been developed to reproduce the wealth of data available in terms of total metallicity, chemical abundance ratios, and radial gradients of these quantities. Besides the galaxy-wide processes listed above, the (generally unknown) star-formation history of a galaxy, its dependence on radius, and many physical effects related to stellar evolution (stellar evolutionary sequences for each initial mass, SN explosions timescale, chemical and dust yields, radial re-distribution of matter, among others) are to be taken into account. All these quantities are often expressed as a function of some parameters whose value is to be obtained by comparison with the observations (e.g., Matteucci and Greggio 1986; Matteucci and Francois 1989; Matteucci and Brocato 1990; Matteucci 1994, 2001; Chiappini and Gratton 1997; Chiappini et al. 2001; Naab and Ostriker 2006; Pipino et al. 2006; Prantzos 2008, 2009; Calura and Menci 2009; Cescutti et al. 2015; Ryde et al. 2016; Vincenzo et al. 2016; Grisoni et al. 2017, 2018).

The treatment of dust production and destruction is an important ingredient that can affect the results especially for the elements most depleted into dust, and change important aspects of the chemical evolution models (Dwek 1998; Calura et al. 2008; Valiante et al. 2009; Gioannini et al. 2017a, b).

A class of models for star-forming galaxies are named “equilibrium” models (or “gas-regulator” models) because they consider how the interplay of all the on-going processes discussed above affects the gas reservoir of the galaxies and, therefore, their star formation, giving a slowly-evolving quasi-steady state, in which gas inflow is compensated by star formation and outflows, yielding a nearly constant, or slowly evolving, gas content (Bouché et al. 2010; Peeples and Shankar 2011; Dayal et al. 2013; Lilly et al. 2013; Forbes et al. 2014; Pipino et al. 2014; Peng and Maiolino 2014a; Feldmann 2015; Yabe et al. 2015b; Harwit and Brisbin 2015; Kacprzak et al. 2016; Hunt et al. 2016b, see Fig. 11).

Fig. 11
figure 11

Image reproduced with permission from Lilly et al. (2013), copyright by AAS

Basic processes of the equilibrium model

From the point of view of chemical evolution, there are still many chemical elements whose evolution is not satisfactorily reproduced in the solar neighbourhood. This is very probably due to failures in the stellar yields suggesting that nucleosynthesis calculations should be revised. Recently, Matteucci et al. (2014) have shown that a possibility of reproducing, for instance, the solar abundance of the r-process element Europium, as well as [Eu/Fe] versus [Fe/H] in the galaxy, is to assume that this element is mainly produced during merging of neutron stars. Such an event has been witnessed for the first time in connection to the gravitational waves event GW170817 (Pian et al. 2017). The other channel of Eu production is represented by core-collapse SNe, although their Eu production is not enough to reproduce the solar Eu. New nucleosynthesis calculations are necessary also for elements such as Mn, Cr, K, Ti (see Romano et al. 2010). Better dust production and destruction prescriptions are expected to better reproduce chemical evolution in the presence of dust.

Summarizing, models of galaxy formation and evolution have been computed using very different analytical and numerical methods, each of them with strong and weak points. All these models use different techniques to study the effect of the same set of critical mechanisms, in particular the processes limiting star formation in galaxies. Total metallicity and element abundance ratios are very sensitive to these processes and to the timescales of star formation, therefore they are among the main observables used to constrain the models.

5 Metallicity scaling relations in galaxies

Metallicity, both of gas and stars, shows clear scaling relations with several integrated properties of the galaxies. These relations are present both in star-forming and quiescent galaxies, and constitute one of the most revealing pieces of information about the evolution of galaxies, as discussed in Sect. 4. The processes of galaxy formation and evolution depend critically on several parameters, at least halo mass and environment; metallicity, being a consequence of all the star formation history, of star formation, gas accretion, merging, and gas outflow, is critically dependent on these parameters in the form of scaling relations.

5.1 The mass–metallicity relation (MZR)

The primary scaling relation of metallicity is observed to be with galaxy stellar mass and, as a consequence, with the quantities that scale with mass, such as luminosity in bands dominated by relatively old stars. The MZR exists for both gas-phase and stellar metallicities.

5.1.1 Stellar metallicity

The stellar MZR was first discovered in local ellipticals by studying colour-magnitude diagrams and stellar spectroscopy (McClure and van den Bergh 1968; Sandage 1972; Mould et al. 1983; Buonanno et al. 1985) and soon interpreted as the result of chemically-enriched SN winds preferentially ejecting metals from low mass galaxies, due to their shallower gravitational potential well (Tinsley 1974, 1978; Mould 1984).

In more recent times, several authors have applied the methods explained in Sect. 2 to optical spectra of local galaxies, in particular to the extensive SDSS spectroscopic database, to derive the stellar MZR of local galaxies, both quiescent and star forming (Trager et al. 2000a; Kuntschner et al. 2001; Gallazzi et al. 2005, 2006, 2008; Thomas et al. 2005, 2010; Panter et al. 2008; Graves et al. 2009; Harrison et al. 2011; Petropoulou et al. 2011; Conroy et al. 2014; González Delgado et al. 2014; Fitzpatrick and Graves 2015; Sybilska et al. 2017; Lian et al. 2018b; Kirby et al. 2013; Zahid et al. 2017; Zhang et al. 2018a). A clear MZR is seen, with metallicity increasing with mass, as illustrated in Fig. 12.

Fig. 12
figure 12

Image reproduced with permission from Zahid et al. (2017), copyright by AAS

Stellar MZRs of local galaxies derived by various authors. Green squares are the median values of local SDSS galaxies with metallicities derived from Lick indices (Gallazzi et al. 2005); black dots are also SDSS galaxies but with metallicities derived from spectral fitting of stacked spectra (Zahid et al. 2017); Blue stars are nearby, single galaxies whose metallicity is based on spectroscopy of individual supergiant stars (Bresolin et al. 2016; Kudritzki et al. 2016; Davies et al. 2017b). Red triangle are local dwarf galaxies (Kirby et al. 2013)

As detailed in Sect. 2, the stellar populations of local and distant galaxies are studied using not only rest-frame optical spectra, but also the UV spectra, especially in starburst galaxies. The UV stellar MZR have been derived using samples of galaxies selected in different ways in the local universe (Heckman et al. 1998; González Delgado et al. 1998; Leitherer et al. 2011; Zetterlund et al. 2015). When UV spectra are used, only the young, massive stars are sampled, and the results are expected to be similar to those derived for the ISM out of which young stars have recently formed (see Sect. 5.1.2).

By deconvolving the integrated spectra it is also possible to derive stellar metallicity as a function of stellar age in the same galaxy (e.g., Panter et al. 2008). With this technique evidence for temporal evolution is found, but the actual significance of the effect is debated because the process of deconvolving spectra is often model-dependent and always associated with significant uncertainties.

Usually, the derived metallicities are luminosity-weighted and can be significantly different from the mass-weighted values provided by models of galaxy evolution (Zahid et al. 2017; Lian et al. 2018b). The observed scatter is often greater than the observational uncertainties, showing that other parameters besides the stellar mass affect the chemical evolution of galaxies.

As detailed in Sect. 7.2.1, the same authors cited above also find a systematic increase of the \(\alpha \)-elements-to-iron ratio with stellar mass. The \(\alpha \) alpha/Fe abundance ratio is a powerful tool to constrain the relative contribution of SNIa and core-collapse SNe, hence of the star formation history timescale, and an enhanced \(\alpha \)-elements/Fe abundance ratio is usually explained as a consequence of shorter formation timescales in massive galaxies, the so-called downsizing (e.g., Thomas et al. 2010; Onodera et al. 2015).

Splitting the SDSS sample into the two classes of quiescent and star-forming galaxies, Peng et al. (2015) found evidence of the central role of “strangulation” or, more generally, “starvation” to explain the stellar metallicity properties of galaxies (Fig. 13). Strangulation is the suppression of gas accretion due to dynamical or physical processes (Larson et al. 1980; Balogh et al. 2000; Kereš et al. 2005). When the infall is halted, the galaxy evolves as a closed box (Searle and Sargent 1972; Tinsley 1980), reducing gas fraction, enriching gas with the residual, decreasing activity of star formation, and producing stars with rapidly increasing metallicities, primarily as a consequence of the lack of inflowing gas that dilutes the metallicity. As observed by Peng et al. (2015), in such a scenario higher metallicities are expected in quiescent (starved) galaxies with respect to star-forming galaxies which are still experiencing infall (hence dilution) of metal-poor gas. A similar interpretation is given by Spitoni et al. (2017) in their detailed analytical models, where the metal enrichment of the passive population relative to the star-forming population is achieved through a faster decline (relative to normal discs) of the inflow rate, which is similar to the simple “starvation” scenario but expressed through a smoother, more physically plausible declining function. Also the recent results of the hydrodynamical code EAGLE are consistent with this scenario (De Rossi et al. 2018). Whatever is its origin, the dominant feedback effect in these galaxy does not remove gas but prevents further accretions.

Fig. 13
figure 13

Image reproduced with permission from Peng et al. (2015), copyright by Macmillan

Stellar mass–metallicity relation for quiescent (red) and star-forming (blue) galaxies in the local universe. The difference is ascribed to a sudden reduction of gas infall, producing a reduction of star formation and a rapid increase in metallicity during the “starvation” phase, due to the lack of fresh supply of external gas diluting the metallicity

The role of environment in shaping the stellar MZR has been the subject of considerable efforts because it is considered, together with mass, one of the main, independent variables driving galaxy evolution (Trager et al. 2000a; Kuntschner et al. 2001; Thomas et al. 2005, 2010; Sheth et al. 2006; Sánchez-Blázquez et al. 2006; Pasquali et al. 2010; Zhang et al. 2018a). In clusters, galaxies on the red sequence show metallicities that do not depend on age or morphology: old, quiescent ellipticals and younger lenticular follow the same MZR with similar dispersion (Nelan et al. 2005; Mouhcine et al. 2011). It appears that while environment has a strong effect in defining the morphology, age, and the overall level of activity of the galaxies (Pasquali et al. 2010; Peng et al. 2010; Peng et al. 2012), the direct (i.e., not mediated by mass) effect on the MZR is modest (Thomas et al. 2010; Mouhcine et al. 2011; Fitzpatrick and Graves 2015; Sybilska et al. 2017). The effect of environment is larger for dwarf satellite galaxies in high overdensities, in the sense that they tend to be more metal rich than in low density environment (Peng et al. 2015; Trussler et al. 2018), an effect ascribed to starvation as these systems plunge in the hot environment of massive overdensities. In the Illustris simulations Engler et al. (2018) find a dependence of stellar metallicity of dwarf galaxies in clusters on the time elapse from the infall into the cluster, and explain this effect with stellar stripping.

Summarizing, increasing amounts of data are available about stellar metallicities in nearby galaxies and their dependence on the other main parameters of galaxies. Clear MZRs are present in all galaxy samples, and models can be tested if they can reproduce these relations in various galaxies, for stars of different ages, and in different environments. The difference in MZR between ETG and star-forming galaxies has been interpreted as evidence of gas starvation.

5.1.2 Gas-phase metallicity

This dependence of ISM metallicity on mass was observed, for the first time, in galaxies of the Local Group (van den Bergh 1968; Peimbert and Spinrad 1970). The first MZR was inferred for a small sample of local, star-forming galaxies (irregulars and blue-compact dwarfs) in an early work by Lequeux et al. (1979) in the form of a dependence of chemical abundance on total (dynamical) mass, and the correlation was soon after confirmed by Talent (1981) and Kinman and Davidson (1981), see Pagel and Edmunds (1981) for an early review. The same authors found a clear anti-correlation of metallicity with gas fraction, which is mainly driven by the anti-correlation between mass and gas-fraction (e.g., Rodrigues et al. 2012; Peeples et al. 2014). Total and stellar mass are difficult to obtain for a large number of galaxies, and this was even more true in the 1980s. As a consequence luminosity was often used as a proxy for mass when studying these scaling relations (Garnett and Shields 1987; Skillman et al. 1989; Garnett 2002; Lamareille et al. 2004).

The quality of the MZR observed in the local universe improved significantly with the use of SDSS spectra that allowed astronomers to measure the flux ratios of the main optical emission lines for more than 100 000 galaxies (e.g., Tremonti et al. 2004; Mannucci et al. 2010; Pérez-Montero et al. 2013; Lian et al. 2015), as shown in Fig. 14. The observed scatter around the relation in SDSS galaxies is of the order of 0.1 dex (Tremonti et al. 2004; Mannucci et al. 2010), somewhat larger than the metallicity measurement uncertainties. The MZR was also extended towards low mass galaxies, rare in the SDSS sample, albeit with a larger scatter and possible biases due to selection effects (Skillman et al. 1988; Lee et al. 2006; van Zee and Haynes 2006; Haurberg et al. 2013; Pilyugin et al. 2013; Haurberg et al. 2015). The existence of the MZR is nowadays established from \(\sim 10^7\) M\(_\odot \) to \(\sim 10^{12}\) M\(_\odot \)), with a steep dependence at low masses, up to \({M}_*\sim 10^{10}\) M\(_\odot \), that then flattens out at higher masses.

Fig. 14
figure 14

Image reproduced with permission from Tremonti et al. (2004), copyright by AAS

Gas-phase MZR in the local universe. Black dots are the median of the distribution of the SDSS galaxies (shown as gray dots). The black lines enclose 68% and 95% of the distribution, and the red line is a polynomial fit to the data. The inset histogram is the distribution of residuals of the fit. This MZR is derived using a metallicity calibration based on photoionization models, and this is likely to produce abundances higher than the actual values, as discussed in Sect. 3.4

Most of these studies are based on metallicities derived with the strong-line method, the only technique that can applied to large numbers of galaxies. It should be noted that both the shape and, even more, the overall normalization of the gas-phase MZR depend on the calibration used. In particular, the MZR based on photoionization models, like Tremonti et al. (2004), and also Mannucci et al. (2010) at high metallicities, provide high normalizations that, for example, do not fit the positions of the MW, the LMC and the SMC which are all significant more metal-poor than these MZR. For this reasons, of particular interest are the MZRs derived by using “direct” metallicities (see Sect. 3.1) (Berg et al. 2012; Andrews and Martini 2013; Ly et al. 2016; Curti 2019). These studies confirm the overall shape of the relations found with the strong line methods but with a significant lower normalization, as shown in Fig. 15.

The contribution from the DIG to the line ratios and, therefore, to the computation of metallicities, can also significantly affect the measured shape of the MZR and therefore is another source of uncertainty (Sanders et al. 2017).

The comparison between stellar and gas-phase MZR can be used to obtain scientific insights on several aspects of galaxies. For example,the small amount of ISM present in early-type galaxies (ETG) shows metallicities similar to that of the old stellar population (Griffith et al. 2018). This reveals that the ISM in these galaxies is not due to infalling gas (expected to be less metal rich than the stars) but to internal production.

Fig. 15
figure 15

Comparison of the MZRs obtained in the local universe (SDSS data, \(\langle z \rangle \sim 0.08\)) using different metallicity calibrations, from Curti (2019). The grey shaded area gives the average metallicity of galaxies using the \(T_{\text {e}}\)-based calibration in Curti et al. (2017), in good agreement with Andrews and Martini (2013) and Pettini and Pagel (2004). Some of the differences, especially at low masses, are due to the different galaxy selection criteria and are linked to the dependence of metallicity on SFR (see below). The curves with higher normalizations are obtained when using calibrations based on photoionization models, while lower curves are based on “direct” \(T_{\text {e}}\)-based metallicities. Stars and crosses show the MZR inferred by measuring the metallicities through BSGs and RSGs (Davies et al. 2017c).

5.1.3 Interpreting the MZR

Various possible driving mechanisms have been proposed to explain the existence and the shape of the MZR.

First, MZR could be shaped by outflows produced by feedback (e.g., Garnett 2002; Brooks et al. 2007). Outflows mainly due to SNe are very common in starburst galaxies, both in the local universe and at high redshifts (e.g., Heckman 2002; Law et al. 2007; Weiner et al. 2009; Steidel et al. 2010; Martin et al. 2012; Heckman and Thompson 2017). These outflows are observed to have metallicities higher than the ISM of the parent galaxies (Chisholm et al. 2018), and are expected to be much more effective in small galaxies, where the potential well is shallower, hence removing a larger fraction of metal-enriched gas from low-mass systems towards the CGM and the IGM (Tremonti et al. 2004; Tumlinson et al. 2011; Chisholm et al. 2018).

Second, it is known that high-mass galaxies evolve more rapidly and at higher redshifts than low-mass ones, the so-called “downsizing” (e.g., Cowie et al. 1996; Gavazzi and Scodeggio 1996; Somerville and Davé 2015, see also Sect. 7.2.1), therefore at the present stage they are expected to have converted a larger fraction of their gas into stars and metals, reaching a higher metallicity (Maiolino et al. 2008; Zahid et al. 2011). Under this interpretation the MZR is a sequence of evolutionary stages.

Third, the earlier evolutionary stage of smaller galaxies and their larger gas fraction (e.g., Erb et al. 2006b; Rodrigues et al. 2012; Lagos et al. 2016a) could be linked to the on-going infall of metal-poor gas, which, once mixed with the existing ISM, contributes to reduce metallicity and to the build up of the stellar population through star formation.

Fourth, the shape of the high-mass end of the IMF could depend on galaxy mass, introducing a systematic change in the average stellar yields and in the rate of metal enrichment (Köppen et al. 2007; Trager et al. 2000a; Mollá et al. 2015; Vincenzo et al. 2016; Lian et al. 2018a).

Fifth, the metallicity of the accreted gas, recycled from previous episodes of star formation, may be larger for larger mass galaxies (Brook et al. 2014; Ma et al. 2016).

The equilibrium models introduced in Sect. 4 are explicitly built to reproduce the MZR and therefore explain this relation as the consequence of the interplay of many on-going processes. The simple “bathtub” models explain the MZR without having to invoke any direct effect of the gravitational potential in terms of capability of retaining metals (i.e., no mass-dependent outflow rate), although mass is indirectly included through the DM framework, which accelerates the evolution of more massive systems (e.g., Bouché et al. 2010; Lilly et al. 2013; Peng and Maiolino 2014a; Dekel et al. 2013; Dekel and Mandelker 2014).

SAM and hydro numerical codes (see Sect. 4) were also tuned to reproduce the MZR and its evolution (De Lucia et al. 2004; Croton et al. 2006; Finlator and Davé 2008; Oppenheimer and Davé 2008; Oppenheimer et al. 2010; Dutton et al. 2011; Davé et al. 2011, 2012; Somerville et al. 2012; Dayal et al. 2013; Forbes et al. 2014; Lu et al. 2015; Pipino et al. 2014; Torrey et al. 2014, 2017, 2018; Zahid et al. 2014c; Feldmann 2015; Harwit and Brisbin 2015; Kacprzak et al. 2016; Christensen et al. 2016; Hirschmann et al. 2016; Rodríguez-Puebla et al. 2016). These models often produce an anti-correlation between metallicity and gas fraction (Schaye et al. 2015; De Rossi et al. 2015; Lagos et al. 2016a; Segers et al. 2016; De Rossi et al. 2017), which is actually observed, see Sect. 5.2.3.

A critical role in shaping the MZR is played by the properties of stellar and AGN feedback, the chemical yields, the actual metallicity of the outflowing wind with respect to the parent ISM, the fraction of metals that are re-accreted, and the evolution of the star formation efficiency. Many of these parameters are degenerate, and in some cases measures of metallicity help to break these degeneracies. For example, by comparing the stellar and gas mass–metallicity relations, Lian et al. (2018b) concluded that only two scenarios can reproduce both relations as well as the MSSF: either strong outflows remove most of the metals, or a steeper IMF characterises the at early stages of galaxy formation. As it will be discussed in Sect. 7, the study of the abundance ratios between different elements is potentially capable of breaking degeneracies because they are sensitive to the timescales of star formation.

An important issue when comparing observations with theory is that of normalization which drives the total amount of metals observed in the ISM. As shown in Fig. 15 and discussed in Sect. 3.4, metallicity calibrations totally or partially based on photoionization models provide a significantly higher normalization of the MZR than the direct methods. The amount of metals that the models must produce and disperse into the IGM depends critically on this problem. The direct method is now considered to be more reliable and the low normalizations should be preferred.

Summarizing, the MZR could be a sequence of metals removal by outflows, enrichment, dilution, or evolutionary stage, and all these effect could be simultaneously present. Information on the relative importance of these effects can be obtained by studying the effective yields (see Sect. 5.1.5).

5.1.4 Redshift evolution of the MZR

Measuring the evolution of the MZR with redshift is not an easy task. A large number of spectra are needed, and measuring metallicities generally require high S/N ratios spectra, especially to measure stellar metallicity but also when measuring the abundances in the gas phase with the “strong-line” method.

The evolution of the stellar MZR has been studied with optical spectroscopy to sample the stellar populations dominating the total stellar mass up to intermediate redshifts in massive galaxies, both in clusters and in the field (Kelson et al. 2006; Ferreras et al. 2009; Choi et al. 2014; Gallazzi et al. 2014; Onodera et al. 2015; Leethochawalit et al. 2018), finding a modest evolution, consistent with the passive evolution of the stellar population. A significant number of rest-frame UV spectra of high redshifts galaxies have been obtained (Mehlert et al. 2002; Fosbury et al. 2003; Shapley et al. 2003; Steidel et al. 2004, 2016; Rix et al. 2004; Savaglio et al. 2004; Halliday et al. 2008; Quider et al. 2009; Dessauges-Zavadsky et al. 2010; Erb et al. 2010; Mouhcine et al. 2011; Sommariva et al. 2012; Faisst et al. 2016). In particular, a number of bright, usually lensed galaxies have been studied in considerable detail in the UV, obtaining a wealth of information about their level of metal enrichment and about the chemical abundances ratios (e.g., Pettini et al. 2001; Villar-Martín et al. 2004; Dessauges-Zavadsky et al. 2010). Despite large uncertainties, as expected the stellar MZR derived from rest-frame UV spectra (hence probing young stars) is similar to the gaseous MZR, has a significant redshift evolution (see Fig. 16), and shows an inverse relation between metallicity and SFR (Faisst et al. 2016).

Fig. 16
figure 16

Image reproduced with permission from Faisst et al. (2016), copyright by AAS

Stellar MZR at \(z\sim 5\) for averages of galaxies with (blue) and without (red) Ly\(\alpha \) emission. The stellar MZR in compared with the gas-phase MZR in Mannucci et al. (2009)

Significant observational efforts by many groups have provided a clear picture of the evolution of the gas-phase MZR up to \(z\sim 3.5\), i.e., out to the maximum redshift where the main optical lines are still in the near-IR bands. This MZR is found to evolve monotonically with redshift, with metallicity declining with redshift at a given mass. At low redshift the evolution is faster at lower mass (see Fig. 17), while high mass galaxies have already reached their current metallicity at \(z\sim 1\), in what appears to be the chemical version of downsizing.

Fig. 17
figure 17

Image reproduced with permission from Zahid et al. (2014a), copyright by AAS

Redshift evolution of the MZR up to \(z=1.55\). Blue, yellow, black and cyan points refer to \(z=0.08\), 0.29, 0.78, and 1.55, respectively. These metallicities have the “high” normalization related to the use of photoionization models, see Fig. 15

This results is based on many observational works at intermediate (\(z\le 1.5\)) redshifts, mainly using optical spectroscopy (e.g., Contini et al. 2002; Kobulnicky et al. 2003; Kobulnicky and Kewley 2004; Maier et al. 2004, 2005; Savaglio et al. 2005; Maier et al. 2006; Cowie and Barger 2008; Zahid et al. 2011, 2013; Moustakas et al. 2011; Cresci et al. 2012; Foster et al. 2012; Pérez-Montero et al. 2013; Yuan et al. 2013; Nakajima et al. 2013; Guo et al. 2016; Pérez et al. 2016; Suzuki et al. 2017; Hoyos et al. 2005). Tracing the evolution at higher redshifts requires near-IR spectroscopy. The much brighter sky background, the reduced atmospheric transmission and the technological limitations of near-IR spectrographs with respect to optical ones, and the fainter apparent brightness of high-redshift objects produce a lower number of useful spectra, usually with lower S/N ratio and for galaxies with higher SFR. Nevertheless a number of authors have produced significant databases at \(z\sim 2\) (e.g., Erb et al. 2006a, 2010; Finkelstein et al. 2011b; Wuyts et al. 2012, 2014a; Cullen et al. 2014; Zahid et al. 2014b; Steidel et al. 2014; Sanders et al. 2015, 2016c, 2018; Onodera et al. 2015; Bian et al. 2017) and \(z\sim 3\) (Maiolino et al. 2008; Mannucci et al. 2009; Belli et al. 2013; Maier et al. 2014; Troncoso et al. 2014; Onodera et al. 2016).

In many cases, the detection of the faintest lines needed to measure metallicities is only possible when many galaxy spectra are stacked together. This procedure is not free from uncertainties because it involves choices on which galaxies to stack, how to do the stacking, and is prone to non-linear effects. Lensed and single bright galaxies have been extensively observed to obtain higher S/N ratio on the spectra of single galaxies, especially in the low mass range (Teplitz et al. 2000; Kobulnicky and Koo 2000; Richard et al. 2011; Wuyts et al. 2014b; Jones et al. 2015b; Perna et al. 2018).

While the resulting big picture outlined above is clear, the details are debated and depend on the method used.

First, metallicity indicators play a critical role, especially because at high redshift each observing program typically has access to a restricted number of diagnostics, hence introducing scatter and systematics among different surveys (Kewley and Ellison 2008, see Sect. 3.5). This is even more important when taking into account that different results can be obtained also when using different line ratios of the same calibration (e.g., Brown et al. 2016). Also, the possibility of a significant evolution with redshift of the calibrations due to different conditions of the star-forming regions, for example in terms of ionization parameter, ionizing spectra, density, pressure, and N/O ratio is always present (see Sect. 3.5.4), and could also be differential for the various methods (Kewley et al. 2013b, a, 2015; Steidel et al. 2014; Shapley et al. 2015; Strom et al. 2017a; Kashino et al. 2017).

Second, the method used to measure stellar mass (and SFR) affects the results, as discussed, e.g., in Yates et al. (2012) and Cresci et al. (2018).

Third, how galaxies are selected affects the final results. For example, when galaxies are selected according to the flux of a metallicity-sensitive emission line, such as [OIII]\(\lambda \)5007 (Izotov et al. 2011, 2015; Xia et al. 2012) the results can be biased toward lower or higher metallicities.

Fourth, similar but more subtle effects are introduced by the need of putting a S/N threshold on the flux the line to be used to measure metallicity. If a relatively high (3–5) minimum S/N is used on all the lines (e.g., Yates et al. 2012), metallicity-dependent selection effects can be introduced near the detection threshold. For example, the flux of [OIII]\(\lambda \)5007 is about 1/10 of H\(\alpha \) at solar metallicities (and low extinction) and about as bright as H\(\alpha \) at \(Z\sim 0.2 {\hbox {Z}}_\odot \). If a threshold of S/N=3 is used for [OIII]\(\lambda \)5007, at each SFR (i.e., roughly at each H\(\alpha \) luminosity), a larger fraction of low metallicity galaxies are included, while high metallicity galaxies are preferentially excluded altering the cosmic average and the resulting MZR. Opposite biases are introduced by [NII]\(\lambda \)6584, whose flux increases with metallicities. The result is the introduction of biases that are difficult to trace (Salim et al. 2014; Cresci et al. 2018). The opposite approach consists in using a high S/N threshold only on the emission line more directly related to SFR and less dependent on metallicity, such as H\(\alpha \) and H\(\beta \), obtaining a more SFR-selected sample (e.g., Mannucci et al. 2010). In both cases the resulting sample is usually not mass- or volume-selected.

Fifth, often spectra are obtained inside a fixed aperture irrespective to galaxy distance. For example, SDSS spectra are obtained with a 3 arcsec circular fiber placed on the galaxy center. This, together with the existence of radial metallicity gradients, can introduce spurious correlations of metallicity with distance and galaxy size that are not easy to estimate and correct.

Finally, the presence of a dependence of metallicity on SFR and other galactic properties, such as size and surface density (see Sect. 5.2), means that the shape of the observed MZR and its redshift evolution depend critically on how the targets are selected in terms of luminosities and redshift range. In other words, the MZR is only defined for a given average SFR at each mass, using intrinsically fainter or brighter galaxies for a given mass and redshift affects the shape of the MZR. Marked differences between published results (e.g., Steidel et al. 2014 vs. Wuyts et al. 2014a) can be explained by this effect, see the discussion in Cresci et al. (2018).

Despite all these uncertainties, there is a general agreement on the fact that the observed MZR evolves with redshift (see Fig. 17), especially at \(z>1\), i.e., the observed metallicity at a given stellar mass decreases with redshift at a rate which depends on redshift and mass. The evolution of the MZR can be parametrized in different ways (e.g., Moustakas et al. 2011). Maiolino et al. (2008) and Mannucci et al. (2009) used a simple 3rd-order polynomial fit, while Zahid et al. (2014c) introduced a different parametrization named “Universal Metallicity Relation”. The variation of only one of the parameters used in this analytic formula, the mass at which the mass-dependence on metallicity starts to flatten out, is often enough to reproduce the observed evolution. Curti (2019) proposed a similar parametrization but with one more parameter to better match the observations.

Most galaxy evolution models cited in the previous section reproduce the evolution of the MZR. The steady decrease of metallicity with redshift at constant mass is ascribed to several reasons, including higher efficiency in ejecting gas and reduced stellar yields. Based on the IllustrisTNG simulations, Torrey et al. (2017) identify one of the main drivers of the MZR evolution in the increasing gas fraction with redshift, in agreement with the anti-correlation between metallicity and gas fraction and SFR observed in the local universe (see Sect. 5.2.3), while Yabe et al. (2015b) attribute the evolution of the MZR to higher infalls and outflows at high redshifts , and Lian et al. (2018a) propose either higher metal loading factors or a steeper IMF at high redshifts. As discussed in Sect. 4, a meaningful comparison with the models should consider all the selection- and observational- effects listed above, and this is not always done.

Summarizing, there is a clear evidence that the gas-phase MZR and the stellar MZR based on UV observations (i.e., related to young stars) are evolving with redshift, with lower metallicities at earlier cosmic times. The details of the evolution depend on a number of possible selection effects and issues regarding how metallicities are estimated. The increase of metallicity with time for a given stellar mass is a feature which is commonly reproduced by the models of galaxy formation.

5.1.5 Effective yields

The effective yields are a way to measure the influence of infall and outflows on the chemical evolution of a galaxy (Matteucci 2001; Garnett 2002; Dalcanton 2007). The differential equations describing a closed-box system, i.e., with no infalling gas and not outflowing winds, with the assumptions of instantaneous recycling, and instantaneous mixing, can be solved (Edmunds 1990) to express the gas metallicity Z as a function of the gas fraction \(f_{\mathrm{gas}}\) and of the true stellar yield y, i.e., the ratio between the amount of metals produced and returned to the ISM and the mass of stars:

$$\begin{aligned} Z = y\ \ln (1/f_{\text {gas}}). \end{aligned}$$
(16)

The true yields y can be expressed in terms of the solar yield \(y_\odot =0.0142\), i.e., the fractional contribution of metals to the solar mass (Asplund et al. 2009).

Using the observed quantities for Z and \(f_{\mathrm{gas}}\), the effective yield is defined as:

$$\begin{aligned} y_{\mathrm{eff}} =Z/\ln (1/f_{\mathrm{gas}}). \end{aligned}$$
(17)

In general the measured values of \(y_{\mathrm{eff}}\) differ from the true stellar yields y if the system is not a closed box. In other words, comparing metallicity with gas fraction gives information on the gas flow from and into the galaxies. An outflow removes gas, hence reduces \(f_{\mathrm{gas}}\) and therefore reduces also \(y_{\mathrm{eff}}\). If the outflow preferentially eject metals, then it also reduces the metallicity, further reducing \(y_{\mathrm{eff}}\). Inflow of metal-poor gas decreases the metallicity and, although it increases the gas fraction, it can be shown that the net effect is to reduce \(y_{\mathrm{eff}}\). Both effects therefore tend to reduce \(y_{\mathrm{eff}}\) with respect to y, and the difference is a measure of how much these gas flows affect the chemical evolution of the galaxies.

Estimating gas masses via the Schmidt-Kennicutt relation, Tremonti et al. (2004) found a significant dependence of yeff, with lower mass galaxies having lower yeff . At variance with this result, Yabe et al., 2015b used the gas masses directly measured by Peeples and Shankar (2011) and metallicities estimated from the MZR , obtaining yeff decreasing with mass, as shown in Fig. 19-left. This shows the dependence of the results on how gas fractions are measured, but in either cases it proves that galaxies are not closed boxes, that infall and/or outflows of gas have an important effect on metallicity, and that earlier evolutionary stages in low mass galaxies cannot be the only explanation for the MZR. Either outflows of metal-rich gas, or infall of metal-poor gas, or both, must have measurable consequences on metallicity.

Fig. 18
figure 18

Image reproduced with permission from Tremonti et al. (2004), copyright by AAS

Effective yields as a function of total baryonic mass for SDSS galaxies (gray) and other low-mass galaxies, from Tremonti et al. (2004). The black diamonds are the median of the data distribution. The dashed line indicates true yield y if no metals are lost. The pink line is a fit with a physically-motivated analytic formula, and the pink stars denote galaxies that have lost 50% and 90% of their metals.

At high redshift yeff is found increase with gas fraction and, as a consequence, to decrease with stellar mass. This is observed at \(z\sim 1\) (Rodrigues et al. 2012; Yabe et al. 2015b), at \(z\sim 2\) (Erb et al. 2006a; Erb 2008; Wuyts et al. 2012), and at \(z\sim 3\) (Mannucci et al. 2009; Troncoso et al. 2014), as shown in Fig. 19, and mass loading factors of infall and outflows significantly larger than 1 are often obtained.

Fig. 19
figure 19

Images reproduced with permission from [left] Yabe et al. (2015b), copyright by AAS; and [right] from Mannucci et al. (2009), copyright by the authors

Left: Metallicity as a function of \(f_{\mathrm{gas}}\) for three sample of galaxies at various redshifts, from Yabe et al. (2015b). The properties of each galaxy sample are reproduced by models with gas infalls and outflows proportional to the SFR, where \(f_\mathrm{{i}}\) and \(f_\mathrm{{o}}\) are the relative loading factors. Right: Metallicity (upper panel) and effective yields (lower panel) as a function of \(f_{\mathrm{gas}}\) in the sample of \(z\sim 2\) by Erb (2008) (magenta stars) and \(z\sim 3\) galaxies from Mannucci et al. (2009) (blue dots). The black dashed line are the expectation of a close-box model (Edmunds 1990; Matteucci 2001, 2008). Blue and red lines are two different fits to the blue data points and emphasize the different contribution of outflows (\(f_\mathrm{{o}}\)) and infalls (\(f_\mathrm{{i}}\)). In this model, infalls are more effective in changing \(y_{\mathrm{eff}}\) for the gas-poor galaxies, while outflows are more effective in gas-rich galaxies

5.1.6 The mass–metallicity relation of DLAs and GRB host galaxies

At even higher redshifts (\(z>3.5\)), all the main optical lines leave the near-IR bands and metallicity can only be obtained with different techniques and on galaxy samples selected in very different ways with respect to the common magnitude- or mass-based selections.

A MZR relation up to \(z\sim 5\) has been derived for gamma-ray burst (GRB) host galaxies by Laskar et al. (2011) by measuring the metallicity through ISM absorption lines. To what degree this is representative of the high-redshift universe is debated because it is still not clear if long-duration GRBs select an unbiased sample of star-forming galaxies. A discussion of this complex issue is beyond the scope of this review. Here we only mention that evidence for a strong metallicity bias has been proposed and negated several times, (see, e.g., Fynbo et al. 2006; Mannucci et al. 2011; Arabsalmani et al. 2015a, 2018; Piranomonte et al. 2015; Trenti et al. 2015; Krühler et al. 2015; Vergani et al. 2017, and references therein).

As explained in Sect. 3.6, DLAs (e.g., Wolfe et al. 1986) provide a unique opportunity to obtain accurate measure of the metallicity of the CGM (and the ISM of the outer disc) in high-z galaxies (see, e.g., Pettini et al. 2002a; Pettini 2006). Comparing this information to that derived from luminosity-selected galaxies is not straightforward (e.g., Fynbo et al. 2008; Christensen et al. 2014), for the various reasons already discussed in Sect. 3.6.

In particular, DLAs are selected on the gas cross-section, and the information derived only applies to the part of the galaxy projected on the background source, whose distance from the center is usually not known. As a consequence, the same galaxy can show very different metallicity properties if, for example, two lines of sight at different radial distances are studied. We recall that it is not even clear whether DLAs probe the extended discs of galaxies, or clumps in the circumgalactic medium, or both [see the discussion in De Cia et al. (2018) and Krogager et al. (2017)].

Moreover, different elements are used to measure the metallicity of the ISM in emission and absorption: typically O and N abundances are estimated for emission-selected galaxies (see Sect. 3), while in DLA various elements can be traced (Zn, S, Fe, Si, O, C, Mg, ...) depending on the observed band and column density, etc.

Finally, the integrated properties of the galaxy associated with the DLA, such as stellar mass and SFR, are often not measured as the galaxy itself is often not even detected, although a smallnumber of DLA galaxies are actually identified in emission (e.g., Rhodin et al. 2018; Kanekar et al. 2018; Krogager et al. 2017; Noterdaeme et al. 2012; Péroux et al. 2011). At \(z\sim 0.7\) these galaxies follow a MZR similar to the emission-selected galaxies when a somewhat uncertain correction for metallicity gradients is included (Rhodin et al. 2018). Moreover, Krogager et al. (2017) identified the optical counterparts of a small sample of DLAs at \(z\sim 2\) and, by combining these with some additional previous detections, suggest that the DLA host galaxies follow a luminosity-metallicity relation.

However, more generally, as a consequence of the issues discussed above, most studies investigate the mass–metallicity relation by using the velocity dispersion of the DLA as a proxy for the mass (Haehnelt et al. 1998; Ledoux et al. 2006; Pontzen et al. 2008). A planar correlation between velocity dispersion (mass), metallicity and redshift has been proposed by Neeleman et al. (2013) that is capable of reducing the scatter about the relations involving only two quantities. The correlation between velocity dispersion and metallicity is found to be in place already at \(z=4\) (Ledoux et al. 2006), and Møller et al. (2013) explored the relation out to \(z=5\) by using a sample of  100 DLA at \(0.1<z<5\) (see also Arabsalmani et al. 2015b, for a confirmation that the same relation applies to DLA along the line of sight of both QSOs and GRBs).

Clearly the absolute normalization of the MZR for DLA and luminosity-selected galaxies cannot be easily compared, as this implies translating the DLA velocity dispersion (at the location of the DLA impact parameter) into a host galaxy stellar mass. However, the relative redshift evolution of the two relations can be compared, and this can give information on the nature of the galaxies associated with the DLA. Figure 20 shows the evolution of the DLA MZR from \(z=5\) to \(z=0\) from Møller et al. (2013) compared to the evolution of the MZR of emission-selected galaxies (Maiolino et al. 2008; Mannucci et al. 2009; Troncoso et al. 2014). Interestingly, the DLA MZR follow the same rapid evolution as low mass galaxies. However, clear evolution is only seen out to \(z=2.6\), beyond this redshift the normalization of the DLA MZR remains constant. The cosmic epoch of this break (\(z\sim 2.6\)) is interesting for various reasons: it is close to the peak of galaxy formation (Madau and Dickinson 2014) and, intriguingly, this is also the redshift where the FMR (see Sect. 5.2) starts to show clear signs of evolution (Mannucci et al. 2010), and it is also the epoch when the comoving HI mass density of galaxies change from rapidly increasing (at \(z>2\)) to remaining constant (at \(z<2\)) (Prochaska and Wolfe 2009). This change has been attributed to several possible effects: a variation of the IMF (Bailin et al. 2010), variations in the properties of the IGM and of equilibrium between accretion from the IGM and gas processing in galaxies (Dixon and Furlanetto 2009; Prochaska and Wolfe 2009; Prochaska et al. 2010), or the switch from “in-situ’ star formation, due to infalling gas, to ”ex-situ” formation followed by accretion of stellar systems already formed (as expected at these redshifts by some cosmological simulations Oser et al. 2010). In some of these scenarios, some of the hypotheses on which the gas-equilibrium models are based (e.g., Davé et al. 2011; Lilly et al. 2013, see Sect. 4) are no longer valid at \(z>2.6\).

Fig. 20
figure 20

Image reproduced with permission from Møller et al. (2013), copyright by the authors

Redshift evolution of the normalization of the mass–metallicity relation for DLAs. Dots show average of DLAs, from Møller et al. (2013), compared to the evolution of luminosity-selected galaxies from Maiolino et al. (2008), Mannucci et al. (2009) and Troncoso et al. (2014). At low masses the MZR of these works is based on “direct” \(T_{\text {e}}\) metallicities; therefore, the comparison is not strongly affected by the use of more modern calibrations (see Sect. 3.4)

Summarizing, DLAs follow a mass–metallicity relation which is also evolving with redshift, sign of ongoing process of chemical enrichment in these systems from high redshifts. The link with the evolution of the MZR in emission-selected galaxies is not well established, but there is similarity between DLA and low-mass galaxies.

5.1.7 Other classes of galaxies

Special classes of galaxies have also been studied to understand the key processes affecting chemical evolution.

Merging and interacting galaxies generally show lower metallicities than the MZR (Kewley et al. 2006a; Michel-Dansac et al. 2008; Ellison et al. 2008b, 2013; Reichard et al. 2009; Rupke et al. 2010a; Morales-Luis et al. 2011; Mouhcine et al. 2011; Torres-Flores et al. 2014; Chung et al. 2013; Cortijo-Ferrero et al. 2017).

This difference is interpreted as the effect of the tidally-driven gas infalls that also produce the increase of SFR (Reichard et al. 2009; Rupke et al. 2010b; Perez et al. 2011; Torrey et al. 2012; Ellison et al. 2013). This results can be used to study the nature of peculiar galaxies. Fox example, the agreement with the local MZR, together with the study of galaxy dynamics and metallicity gradients, allowed Bournaud et al. (2008) to conclude that a chain-like, clumpy galaxy is not an on-going merger but is actually a clumpy disk.

Starburst galaxies, Ultra-Luminous Infrared Galaxies (ULIRGs) and low-redshift analogues of Lyman-Break Galaxies (LBG) are usually found to have metallicity lower than expected for their mass (Liang et al. 2004; Rupke et al. 2008; Roseboom et al. 2012; Lian et al. 2015), and this can be explained by the effect of higher SFR discussed in the next Sect. (5.2).

Green peas (GP) are compact, star-forming, galaxies selected for the presence of a bright [OIII]\(\lambda \)5007 line (Cardamone et al. 2009). As the [OIII]\(\lambda \)5007 flux increases with decreasing metallicity (see Sect. 3.5), these galaxies, as other line-selected samples, are preferentially selected to have low metallicities and, as a consequence, are usually below the MZR (Amorín et al. 2010; Izotov et al. 2011; Xia et al. 2012; Amorín et al. 2012; Ly et al. 2015; Lofthouse et al. 2017; Senchyna and Stark 2018). Similar to the GPs, the extremely metal poor (XMP) galaxies are rare objects selected to have extreme line ratios and very low metallicities, down to a few percent solar in the local universe. They usually have low masses, high sSFR, and disturbed morphologies (Izotov et al. 2006b; Izotov and Thuan 2007; Morales-Luis et al. 2011; Izotov et al. 2012; Sánchez Almeida et al. 2015, 2016; Izotov et al. 2018b, c). By definition they fall below the MZR.

Lyman-alpha galaxies (LAG), i.e., galaxies selected from their Ly\(\alpha \) emission, are found to have low metallicities. The brightness of the Ly\(\alpha \) line in these galaxies is therefore interpreted in terms of low column density of dust associated with the low chemical abundance, which makes it easier to the Ly\(\alpha \) photons to escape (Finkelstein et al. 2011a, b; Nakajima et al. 2013; Song et al. 2014; Trainor et al. 2016).

5.2 Metallicity dependence on environment

The gas-phase MZR shows also some dependence on environment. In the local universe gas-phase metallicities are observed to correlate weakly with environment, with higher values in denser environments (Shields et al. 1991; Skillman et al. 1996; Mouhcine et al. 2007; Cooper et al. 2008; Ellison et al. 2009; Petropoulou et al. 2011, 2012; Kulas et al. 2013; Hughes et al. 2013; Pilyugin et al. 2017; Wu et al. 2017). However, the environmental effects are clearly seen only if satellite and central galaxies are considered separately. Indeed, central galaxies do not show significant correlation with the environmental density, while satellite galaxies do show, at a given stellar mass, a significant metallicity dependence on environment, i.e., satellite galaxies in denser environments are characterized by higher metallicities (Pasquali et al. 2012; Peng and Maiolino 2014b; Peng et al. 2015). In galaxy clusters the metallicity of satellite galaxies seems to also correlate with the stage of accretion of galaxies into the cluster (Maier et al. 2016). These trends have been interpreted as a combination of “strangulation” (i.e., satellite galaxies being prevented from accreting cold, near-pristine gas as they plunge into the hot halo of massive environments, see Sect. 5.1.1), ram-pressure stripping, higher metallicities in the infalling gas, and external pressure reducing the amount of metals lost to the CGM due to gas re-accretion (Pasquali et al. 2010; Peng and Maiolino 2014b; Spitoni 2015; Segers et al. 2016; Pilyugin et al. 2017). The correlation is probably also mediated by the gas fraction, which is higher in galaxies in low-density environments (Wu et al. 2017). Based on the results of the Illustris simulation, Genel (2016) attributes the metallicity difference partly to the different SF history, partly to the smaller size of forming disks, biasing star formation toward the inner, more metal rich parts. In contrast, Bahé et al. (2017) see the effect in the EAGLE simulation and identify gas stripping and suppression of metal-poor infalls as the main drivers of the effect.

At high redshift the situation if even more uncertain, with contrasting results and, in general, little evidence for the existence of any environmental dependence (Magrini et al. 2012; Kulas et al. 2013; Williams et al. 2014; Shimakawa et al. 2015; Kacprzak et al. 2015; Valentino et al. 2015). However, the statistics in these high-z studies are still poor, while it is well known that large statistics are required to identify the role of the environment and disentangle it from mass segregation effects.

Summarizing, the galaxy metallicity (and the metallicity scaling relations) does not depend significantly on environment for central galaxies, while satellite galaxies tend to have higher metallicities in denser environments. The latter phenomenon is possibly associated with the “strangulation/starvation” of galaxies as they plunge into hot halo of massive environments, but other phenomena may also play a role.

5.3 The dependence of metallicity on SFR and gas fraction

As soon as a sufficient precision was reached, evidences for further dependences of metallicity on other galactic properties started to appear. Tremonti et al. (2004) reported a correlation of the metallicity residuals from the MZR with galaxy colour, ellipticity, and central mass density. Soon after, Hoopes et al. (2007) pointed out a dependence of the shape of the MZR on galaxy size, and Ellison et al. (2008a) was the first to explicit state the dependence of metallicity on SFR for a given mass. Mannucci et al. (2010) introduced a 3D relation between mass, metallicity and SFR, named Fundamental Metallicity Relation (FMR), so that in the local universe the residual metallicity scatter across the median relation is reduced and becomes very low, \(\sim \) 0.05 dex (i.e., \(\sim 12\%\)), consistent with the uncertainties of the measurements. For a given stellar mass, metallicity decreases with SFR and sSFR, i.e., more actively star-forming galaxies have lower metallicities than more quiescent galaxies (Fig. 21). As introduced in Sect. 5.1.1, recently a dependence on SFR was also found for stellar metallicities by Faisst et al. (2016).

Fig. 21
figure 21

Image reproduced with permission from Mannucci et al. (2010), copyright by the authors

The fundamental metallicity relation. Left: dependence of the gas metallicity on mass, in bins of SFR. The gray shaded areas contains 68% and 95% of the full, unbinned galaxy sample. Right: dependence of metallicity on SFR in bins of stellar mass

Apparently at odds with the clear evolution of the MZR, Mannucci et al. (2010) also showed that the FMR does not evolve with redshift up to \(z=2.5\) (Fig. 23). All the data at \(z<2.5\) available to Mannucci et al. (2010) (Savaglio et al. 2005; Shapley et al. 2005; Erb et al. 2006a; Liu et al. 2008; Epinat et al. 2009; Wright et al. 2009; Förster Schreiber et al. 2009; Law et al. 2009; Lehnert et al. 2009) closely follow the surface defined by the local SDSS. In contrast, data at higher redshifts (Maiolino et al. 2008; Mannucci et al. 2009) show much lower metallicities, with a difference of about 0.6dex. This difference is discussed in Sect. 5.2.1.

The actual shape of the local FMR depends on several factors, such as how galaxies are selected and how mass, SFR, and metallicities are measured. As a result, relations with different shapes have been published (e.g., Lara-López et al. 2010; Hunt et al. 2012; Yates et al. 2012), while other authors (Brisbin and Harwit 2012; Nakajima and Ouchi 2014) obtained a shape consistent with that derived by Mannucci et al. (2010) (Fig. 21). A FMR using “direct” metallicities, based on the \(T_{\text {e}}\) method, was derived by Andrews and Martini (2013), who found a significantly stronger dependence on SFR than in Mannucci et al. (2010). To detect the faint auroral lines needed to measure \(T_{\text {e}}\), these authors stack SDSS spectra according to mass and SFR, i.e., assume that metallicity depends only (or mainly) on these two parameters. As discussed in Sect. 3.5, Curti et al. (2017) derived a new \(T_{\text {e}}\)-based calibration using a different stacking scheme, based on similarities of the spectra (same [OII]\(\lambda \)3727/H\(\beta \) and [OIII]\(\lambda \)5007/H\(\beta \)ratios) rather then on the galaxy parameters. Using these calibrations, (Curti 2019) derived the corresponding MZR and FMR . Salim et al. (2014), Telford et al. (2016) and Cresci et al. (2018) critically re-analyzed the FMR using different metallicity indicators, various ways to measure SFR, studying the further dependence on galaxy size, and considering the SFR distance of the galaxies from the MSSF as basic parameter, finding similar results to Mannucci et al. (2010).

A number of authors have studied the dependence of metallicity on mass and SFR looking for linear relations or through the use of principal component analysis (PCA), i.e., using a technique that rotates the axes in a multi-dimensional space to minimize the scatter; being a simple linear transformation, the PCA technique cannot really account for more complex, non-linear correlations, such as the MZR and the FMR. Within the same context, using the SDSS sample Lara-López et al. (2010) proposed the existence of a “fundamental plane” between mass, SFR, and metallicities. In contrast to virtually all other authors, Lara-López et al. (2010) use a different approach in which they derive the mass as a function of SFR and metallicity, deriving a plane with a large scatter, 0.16 dex. The FMR, which has a curved shape at high masses and different dependence on SFR for different masses, can be approximated by a plane by reducing the upper limit in redshift which has the effect of lowering the number of high-mass galaxies present in the sample. The resulting plane is quantitatively very different from the FMR in Mannucci et al. (2010). This Fundamental Plane was later revised and extended toward higher masses by Lara-López et al. (2013), confirming the correlation between the three quantities albeit with a larger scatter than the FMR. Hunt et al. (2012) used a composite sample of metal poor, starburst galaxies at \(0<z<3.4\), including the many GPs, with metallicities measured in various ways, finding a planar correlation of metallicities with mass and SFR at any redshift. A new PCA was computed by Hunt et al. (2016a) using a different and larger sample including a number of high-z objects, testing various metallicity estimators, and comparing with the local SDSS galaxies. The linear relation derived by Hunt et al. (2016a) does not evolve with redshift. As discussed by Cresci et al. (2018) it does not reproduce the local dependence of metallicity on mass and the detailed MZR at high redshifts.

Surveys of local galaxies based on the use of two large, IFU spectrographs, CALIFA (Sánchez et al. 2012) and MaNGA (Bundy et al. 2015), were used to investigate the existence and properties of the FMR in these nearby, well-resolved galaxies. A series of papers based on these data (Sánchez et al. 2013; Hughes et al. 2013; de los Reyes et al. 2015; Barrera-Ballesteros et al. 2017; Sánchez et al. 2017) have questioned the existence of the relation in this sample, and hypothesized that the dependence on SFR observed by many authors on SDSS and other data is due to aperture effects. In these works a metallicity is computed for each HII region detected, a radial gradient is computed, and the value of metallicity at the effective radius is used as representative metallicity of the galaxy. Salim et al. (2014) re-analyzed the same data finding opposite results, i.e., the presence on an anti-correlation between metallicity and sSFR, and Pilyugin et al. (2013) excluded the importance of the aperture effects. As discussed in Cresci et al. (2018) and shown in Fig. 22, a clear, monotonic dependence of the metallicity on SFR is actually present in the CALIFA and MaNGA data at the expected level. More recently also Belfiore et al. (2018, in prep.) have re-analyzed the MaNGA data, finding a clear dependence on SFR, also on resolved scales, as discussed more extensively in Sect.6.

Fig. 22
figure 22

Image reproduced with permission from Sánchez et al. (2017), copyright by the authors

Dependence of metallicity on mass and SFR for the 612 CALIFA galaxies in Sánchez et al. (2017). Galaxies are colour-coded with SFR as indicated in the Colour bar. Lines are the best-fitting MZR in different SFR bins

In the local universe, the FMR was found to hold also for Ly\(\alpha \) local analogs (Lian et al. 2015) and Herschel-selected starburst galaxies (Roseboom et al. 2012), with properties similar to the SDSS galaxies. This is somewhat surprising as the conditions in these galaxies are expected to be quite different from the more common galaxies dominating the objects in the SDSS. In contrast, the central regions of barred galaxies observed in the SDSS do not follow the FMR. The radial motions of the gas produced by the bar creates a temporary central increase of SFR which is not linked to a decrease of metallicity. The metallicity is actually observed to increase, even if the amount of this increase, of the order of 0.02–0.06 dex, is debated and depends on the indicator used and on the element considered (Ellison et al. 2011; Cacho et al. 2014). This metallicity increase is explained with accretion of metal-rich gas toward the central regions of the galaxies (Martel et al. 2013, 2018). As noted in Sect, 5.1.7, interacting galaxies show lower metallicities and higher SFR, in qualitative agreement with the FMR (e.g., Ellison et al. 2013). Nevertheless, they tend to have large metallicity dispersion and some offset from the FMR toward lower abundances (Grønnow et al. 2015; Bustamante et al. 2018).

Summarizing, the gas-phase metallicity not only correlates with mass but also anti-correlates with SFR, i.e. more star forming galaxies showing lower metallicities. The actual shape of this relation depends on how galaxies are selected and mass, SFR and metallicities are measured. This relation seems not to evolve up to \(z=2.5\), as detailed in the next section (Fig. 23).

Fig. 23
figure 23

Image reproduced with permission from Cresci et al. (2012), copyright by the authors

Difference between the metallicity predicted by the local FMR in Mannucci et al. (2010) and the metallicity actually observed in a number of galaxy samples at various look-back times (lower scale) or redshifts (upper scale). This diagram shows that no evolution is observed up to \(z\sim 2.5\). The evolution of the MZR can be explained by the selection of galaxies with progressively higher SFR at higher redshifts

5.3.1 Redshift evolution of the FMR

The issue of the evolution of the FMR, or lack thereof, with redshift has been the subject of much work, both observational and theoretical. The dependence of metallicity on both mass and SFR was measured in many galaxy samples at various redshifts, and the results compared with the predictions based on the local FMR in Mannucci et al. (2010).

Many works found good agreement between the local FMR and metallicities at \(z\sim 1\) (Yabe et al. 2012, 2014; Cresci et al. 2012; Stott et al. 2013; Henry et al. 2013a, b; Maier et al. 2015a, b; Calabrò et al. 2017) and at \(z\sim 2\) (Nakajima et al. 2012; Belli et al. 2013; Nakajima and Ouchi 2014; Maier et al. 2014; Stott et al. 2014; Song et al. 2014; Yabe et al. 2015a; Salim et al. 2015; Kacprzak et al. 2016; Wuyts et al. 2016; Sanders et al. 2018; Hirschauer et al. 2018). Lensed galaxies also played an important role in studying the evolution of the FMR over a larger redshift range and toward lower masses and lower SFR, usually finding good agreement with the FMR, at least for the average metallicity level (Richard et al. 2011; Christensen et al. 2012; Wuyts et al. 2012). Some authors confirm the relation between metallicity and SFR, but with a weak dependence or with some evolution with respect to the local FMR (Niino 2012; Pérez-Montero et al. 2013; Grasshorn Gebhardt et al. 2016; Cullen et al. 2014; Zahid et al. 2014b; Salim et al. 2015; Wu et al. 2016). Other works do not have the sensitivity to confirm or discard evolution (Divoy et al. 2014; de los Reyes et al. 2015).

Several works on high redshift galaxies (Wuyts et al. 2012; Steidel et al. 2014; Yabe et al. 2014, 2015a; Sanders et al. 2015; Wuyts et al. 2014a; Onodera et al. 2015; Wuyts et al. 2016) do not find a clear dependence of metallicity on SFR within their data. It should be considered that the comparison of the metallicity properties of distant galaxies with the FMR of the local (SDSS) galaxies must take into account two separate issues. The first issue is whether the local FMR is able to predict the average metallicity of a galaxy sample with a given average mass and average SFR, i.e., whether high-redshift galaxies on average follow the local FMR or not; the second issue is whether within the high-redshift galaxy sample it is possible to detect the metallicity dependence on the SFR as in the local universe. The latter issue is considerably more difficult to tackle than the former because it requires accurate mass, SFRs and metallicities for each galaxy, with final uncertainties smaller than the expected dependence of metallicity on SFR, and a large dynamic range of SFR. Measurements of the SFR at high-z are often quite uncertain due the effect of dust extinction. Similarly, accurate metallicities are often difficult to obtain at high-z. Moreover, at high-z it is difficult to achieve a large dynamical range in mass and SFR, especially at low-masses where the dependence of metallicity on SFR in large.

About the first issue, Cresci et al. (2018) showed that the papers quoted above that do not find an internal dependence of metallicity with SFR are actually confirming the predictions of the FMR in terms of average metallicity. This is a remarkable result given that the FMR is only based on local galaxies. About the second issue, a few works report the detection of the dependence of metallicity on SFR also within the high redshift data alone. Zahid et al. (2014b), Salim et al. (2015), Kacprzak et al. (2016) and Kashino et al. (2017) find lower metallicities in galaxies with higher SFR at \(z\sim 2\), albeit with uncertainties due to the low S/N, the presence of galaxies with no detections in the faint lines, and the large metallicity scatter. Finally, in the context of the MOSDEF survey Sanders et al. (2018) investigated the evolution of the FMR using a large sample (260 objects) of star-forming galaxies at \(z\sim 2.3\) using several emission-line ratios. In contrast to earlier works by the same group based on smaller data sets (Sanders et al. 2015), this study detected the explicit dependence of metallicity with SFR for a given stellar mass within the data sample itself without the need of a comparison with lower redshift samples, see Fig. 24. The level of metal enrichment of galaxies with a given stellar mass and SFR is similar to that in the local universe, being only \(\sim 0.1\) dex lower, consistent with the no-evolution upper limit derived by Cresci et al. (2018).

Fig. 24
figure 24

Image reproduced with permission from Sanders et al. (2018), copyright by AAS

Observed dependence of metallicity on the sSFR at \(z=2.3\). The plot shows the metallicity difference \(\varDelta \)log(O/H) from the mean MZR as a function of the sSFR difference \(\varDelta \)sSFR from the MSSF. The metallicity shown is computed with the N2O2 strong-line indicator, and the original paper also analyses other indicators, finding similar results. Blue points are individual galaxies, red squares show medians in bins of \(\varDelta \)log(sSFR), stars display the results from stacked spectra, colour-coded by stellar mass. The solid line shows the best fit to the medians, and the dashed line shows the prediction from the cosmological simulations of Davé et al. (2017). A clear dependence of the residuals with sSFR is seen

As mentioned, at \(z>2.5\) various works have found that galaxies deviate from the FMR by being more metal poor than expected (Mannucci et al. 2010; Troncoso et al. 2014; Onodera et al. 2016), at least in the redshift range \(2.8<z<3.6\) where optical metallicity diagnostics are still detectable in the K-band. The observed difference, albeit large (\(\sim \) 0.6 dex) should be taken with care: only a few galaxies have reliable metallicities measured at these redshifts, these galaxies often have very high SFRs and could be in special phases of their evolution, metallicity is usually measured using only one indicator (often R23), and the N2 indicator, which is most often used at lower redshift, is no longer available. Nevertheless is interesting to note that this redshift range \(z\sim 3.5\) where the FMR appears to break down is the same range where the DLA scaling relations seem to change behavior (Møller et al. 2013) and where the cosmic density of HI start to evolve rapidly with respect to the non-evolution at lower redshifts (Prochaska and Wolfe 2009). As mentioned, this is the epoch beyond which galaxies may no longer be in equilibrium, i.e., at \(z>2.5\) probably the infall of (pristine) gas on most galaxies is too fast for them to efficiently transform it into stars, hence most of the accreting gas is being accumulated and resulting into a larger dilution of the metals in the host galaxy.

It would be interesting to extend these studies to even higher redshifts. The use of (rest-frame) optical diagnostics shall await JWST. However, the UV nebular lines of CIII]1909, CIV1549 and OIII]1808, although weak, offer an alternative potential tool to trace the gas metallicity out to very high redshift. Interestingly, the few current detections of these lines at high redshift indicate moderately low metallicities (\(z\sim 0.1\)) even out to \(z\sim 7\) (Stark et al. 2017), while at such high redshift one would expect lower levels of metal enrichment, at least for ‘normal’ galaxies. Yet, one must take into account that these early detections may be biased toward higher metallicities.

It is interesting that far-IR fine structure atomic transitions, such as [CII]158\(\upmu \)m and [OIII]88\(\upmu \)m, redshifted into the millimetre bands of atmospheric transmission, are increasingly being used to trace the metal enrichment in ‘normal’ galaxies even at \(z>7\) (Maiolino et al. 2015; Pentericci et al. 2016; Inoue et al. 2016; Carniani et al. 2017, 2018; Hashimoto et al. 2018c; Tamura et al. 2018) and out to \(z\sim 9\) (Hashimoto et al. 2018c). In some of the galaxies the very high [OIII]88\(\upmu \)m–to–[CII]158\(\upmu \)m ratio is interpreted as indication of low metallicity (\(z\sim 0.1\), Inoue et al. 2016), although this interpretation is subject to degeneracies, as the same high ratio can also be interpreted in terms of ionization parameter (Carniani et al. 2017; Katz et al. 2017). Very interestingly, for some of these galaxies the thermal dust continuum is also detected (Hashimoto et al. 2018b; Tamura et al. 2018) which implies a significant amount of dust, and, therefore, of metals, already in place at such early epochs. Some of these authors infer metallicities as ‘high’ as \({Z} \sim 0.2\) at \(z\sim 8\), which requires a very fast and efficient enrichment process in these low mass (\({M}_{{star}}\) a few times \(10^9~\hbox {M}_{\odot }\)) primeval galaxies.

Summarizing, most of the works confirm or are consistent with the existence of the FMR at high redshift, and find no or a very limited evolution of this relation with redshift up to \(z=2.5\). The situation at even higher redshift is unclear but there are several indications of evolution during the first \(\sim \) 2 Gyrs after the Big Bang. This stability of the FMR is an important observation that models must reproduce, as discussed in the next section.

5.3.2 Origin of the FMR

The origin of the dependence of metallicity not only on mass but also on SFR is debated. Ellison et al. (2008a) proposed a varying star formation efficiency as the most probable explanation for the dependence of metallicity on SFR and galaxy size. Mannucci et al. (2010) proposed that the interplay of infall of metal poor gas and star formation may play a central role in shaping the FMR: on the one hand, infall provides chemically poor gas, lowering metallicity; on the other hand gas accretion delivers additional fuel for star formation, hence enhancing the SFR. As such, the FMR is yet another piece of evidence for the importance and ubiquity of cold gas accretion as a dominant driver of galaxy evolution. As presented in Sect. 5.1.7, the same kind of explanation had been previously proposed to explain the deviations from the MZR of several classes of galaxies, such as interacting systems, ULIRG and Green Peas. The spatially-resolved inverse dependence between SFR and metallicity observed within galaxies, both in the local universe (Sánchez Almeida et al. 2013, 2014, 2018; Belfiore et al. 2018) and at high redshift (Cresci et al. 2010) is also naturally interpreted as the effect of the infall of chemically poor gas, as discussed in Sect. 6.

In a more consistent and quantitative way, the FMR is predicted or explained by the gas-equilibrium models described in Sect. 5.1.3. This is in agreement with the observations that local galaxies that are expected to be far from equilibrium, such as interacting pairs and barred galaxies, show quantitative disagreement wit the FMR. Some models (e.g., Lilly et al. 2013) predict no evolution of the FMR because the equilibrium is based on basic physical processes that can remain stable with cosmic time, while others (e.g., Davé et al. 2011) expect a mild evolution with redshift due to the progressive enrichment of the CGM. The observational results are not precise enough to distinguish between these two scenarios in a robust way (Sanders et al. 2018). An infall with a super-linear dependence of SFR on infalling mass is the explanation proposed by Brisbin and Harwit (2012). In the context of the equilibrium models, Wuyts et al. (2016) modified the model by Lilly et al. (2013) by introducing the evolving gas depletion time obtained by Genzel et al. (2015), obtaining a good fit to the evolution of the MZR starting from the local FMR, see Fig. 25. The modest evolution of \(\sim 0.1\) dex of the FMR found by Sanders et al. (2018) at \(z\sim 2.3\) can be interpreted either as a limited but systematic increase in the mass loading factors of the outflowing gas for a given stellar mass, or with the decrease with redshift of the average metal abundance of the inflowing gas (Davé et al. 2011).

Fig. 25
figure 25

Image reproduced with permission from Wuyts et al. (2016), copyright by AAS

Redshift evolution of the MZR at a stellar mass \(M=10^{10}\) M\(_\odot \). The gray shaded areas are the predictions of the Lilly et al. (2013) model with the variable depletion time by Genzel et al. (2015), calibrated on the local FMR by Mannucci et al. (2010)

Lagos et al. (2016b) and De Rossi et al. (2017, 2018) find a very good match between the observations of the FMR and their hydrodynamical codes up to high redshift, and propose gas fraction as the parameter that is better correlated with metallicity (see next section). Critical parameters to reproduce the observations are the properties of stellar and AGN feedback and the dependence of star formation on metallicity. In the IllustrisTNG simulation Torrey et al. (2018) find the FMR up to high redshift (see Fig. 26) and ascribe its origin to the similar timescales of the evolution of SFR and metallicity. This opens up the interesting possibility of using the observed scatter across the FMR to test the existence of variations of the properties of high redshift galaxies on timescales shorter than cosmic evolutionary times.

Fig. 26
figure 26

Image reproduced with permission from Torrey et al. (2018), copyright by the authors

Dependence of metallicity on mass and sSFR produced at z=1 by the IllustrisTNG simulation. Solid and dashed black lines indicate the median MZR and one sigma scatter, respectively. Colours show how the residuals about the MZR are correlated with SFR. This model produces a low gas retention in the ISM in the local universe, where \(\sim 85\%\) of the metals are outside the ISM

Differences among the various chemical elements can also have a role in defining and explaining the shape of the FMR. Recently, Matthee and Schaye (2018) computed the dependence of the abundance of various elements on mass and SFR in the results of the EAGLE simulation. They found a FMR for the abundance of every element, with a smaller dispersion for the alpha elements, and they also find a 3D relation between mass, SFR and [O/Fe].

Interestingly, Pérez-Montero et al. (2013) and Kashino et al. (2016) showed that the dependence of metallicity on mass seen in SDSS galaxies when using R23 or N2, as in Mannucci et al. (2010), is not present when the N2S2 or N2S2H\(\alpha \) indicators (Dopita et al. 2016) is used. Telford et al. (2016) made a critical analysis of these studies, finding that, in contrast, the FMR is also present in the N2S2 indicator, albeit with a weaker dependence on SFR. If confirmed, this effect would be in agreement with the explanation of the FMR as due to infall of metal-poor gas: N2S2 is based on the ratio between two metal lines, therefore is not sensitive to dilution by pristine gas as both abundances go down simultaneously. In contrast, N2 and R23 are based on ratios between metal, collisionally-excited lines and H recombination line, therefore sensitive to dilution.

5.3.3 The relation between metallicity and gas fraction

According to the explanations of the FMR given by the gas-equilibrium models, the link between metallicity and SFR is actually due to two relations: the increase of SFR with gas content (or, equivalently, with gas fraction, \(f_{\text {gas}}\)), and the decrease of metallicity with \(f_{\text {gas}}\) due to dilution. If this is true, a direct dependence of metallicity on \(f_{\text {gas}}\) is expected, and several models obtain this relation as the most fundamental (Lagos et al. 2016a; Segers et al. 2016; De Rossi et al. 2017).

Peeples et al. (2008) and Peeples et al. (2009) showed that such a correlation is present in the local universe, in the sense that SDSS galaxies with low gas-fractions also have metallicities above the MZR. In the formalism of the FMR, Bothwell et al. (2013), Hughes et al. (2013), Lara-López et al. (2013), Jimmy et al. (2015) and Brown et al. (2018) showed that metallicities anti-correlate with the mass of atomic hydrogen (or, equivalently, with the fraction of mass in atomic hydrogen), and that this correlation is tighter than with SFR. Interestingly, the dependence of metallicity on HI gas fraction persists at high stellar mass, where metallicity does not depend on SFR any more (see Fig. 27; Bothwell et al. 2013). Bothwell et al. (2016a, b) showed that an even tighter correlation can be present with molecular hydrogen, more directly related to the SF activity, although the statistics for the molecular gas in galaxies (for which information on the metallicity is available) are still poor, hence the result needs to be confirmed with larger sample of galaxies having measurements for both the gas metallicity and for the molecular gas content.

Fig. 27
figure 27

Image reproduced with permission from Bothwell et al. (2013), copyright by the authors

Observed dependence of gas-phase metallicity on HI mass, in bins of stellar mass, for the SDSS galaxies

5.4 Metallicity dependence on other physical properties

Besides SFR and gas fraction, other dependencies have been proposed to further reduce the scatter of the MZR or explain the evolution of metallicity with respect to the FMR at \(z>2.5\).

Hoopes et al. (2007) and Ellison et al. (2008a) found a dependence of the gas metallicity on size, according to which, at fixed mass, more compact galaxies are more metal rich. The actual difference between the more compact and more diffuse galaxies depends on mass and is about 0.05–0.2 dex. This relation was later analyzed by Brisbin and Harwit (2012). Wu et al. (2015) and Hashimoto et al. (2018a) found similar dependence of metallicity on quantities related to size, i.e., respectively, galaxy surface brightness and stellar surface density (see also Sect. 6). Yabe et al. (2012) found a similar relation at \(z\sim 1.4\) albeit with a large scatter. The origin of this dependence is discussed by all the papers cited above and generally attributed to the effect of gas infall (Ellison et al. 2008a; Brisbin and Harwit 2012). Wu et al. (2015) showed that the dependence is not mediated by gas fraction, and attribute the effect to the dependence on galaxy mass and surface brightness of infalls, outflows, and their relation to SFR. Sánchez Almeida and Dalla Vecchia (2018) investigated the origin of this effect by using the EAGLE cosmological simulations (see Sect. 4). The relation between mass, metallicity and size is present in the results of the simulation (see Fig. 28), together with the FMR, and concluded that the driving mechanism is the infall of metal-poor gas.

Fig. 28
figure 28

Image reproduced with permission from Sánchez Almeida and Dalla Vecchia (2018), copyright by AAS

Predicted dependence of gas-phase metallicity on stellar mass for galaxies of different size from the EAGLE cosmological simulations

The dependence on both mass and radius was further discussed by D’Eugenio et al. (2018) by studying the relation of metallicity with average gravitational potential \((\varPhi \equiv M_*/R_\text {e})\) and average surface mass density \((\varSigma \equiv M_*/R^2_\text {e})\). They found a more direct relation of metallicity with \(\varPhi \), in agreement with the expected dependence of metallicity on escape velocity, giving a central role to metal losses due to galactic winds. The dependence on observing aperture is discussed in Sect. 6.4.

As discussed in Sect. 3.5.4, the ionization parameter q is an important quantity to understand the nebular spectra of galaxies. This parameter shows some correlation with stellar mass and SFR (e.g., Dopita et al. 2006; Brinchmann et al. 2008; Nagao et al. 2006a). Its influence is often studied placing the galaxies in a plane defined by R23, mainly sensitive to metallicity (but with secondary dependence on q), and O32, sensitive to q but with a secondary dependence on metallicity (e.g., Kewley and Dopita 2002; Nakajima et al. 2013; Shapley et al. 2015, see Sect. 3.5 for the definitions). Based on a sample of local and high-z galaxies, Nakajima and Ouchi (2014) analyzed the relation between the residuals from the MZR and the ionization parameters as measured by O32 and by photoionization models. On the base of a correlation between q and the main integrated parameters of the galaxies, i.e., mass, SFR and metallicity, they extended the FMR into a four-dimensional relation including q, the “Fundamental Ionization Relation” (FIR). This relation is able to account for the lower metallicities at \(z>2.5\).

6 Metallicity gradients in galaxies

6.1 Overall properties of metallicity gradients in galaxy discs

Since the first seminal works by Aller (1942), Searle (1971) and Pagel and Edmunds (1981), the investigation of gradients of metallicity in galactic discs has been subject to continuously growing interest, especially in recent years thanks to the developments of new facilities and extensive surveys that have enabled the measurement of gradients in detail using multiple tracers and also over large samples of galaxies.

In the local Universe the general finding is that, at least within the optical radius, in most spiral galaxies the metallicity decreases exponentially with galactocentric radius. As a consequence, gradient of the metallicity in its logarithmic form [e.g., 12 + log(O/H)] is linear with radius, and can be expressed in units of \(\hbox {dex}~\hbox {kpc}^{-1}\). In the case of the Milky Way, typically radial metallicity gradients range between \(-\,0.06\) and \(-\,0.01\ \hbox {dex}~\hbox {kpc}^{-1}\), depending on the metallicity tracer adopted, as discussed in the following.

Indeed, different metallicity tracers have been used to investigate the metallicity gradients, providing different kind of information, especially because they trace the enrichment at different cosmic times, hence their comparison can provide precious information on the evolution of the gradients with time.

HII regions are used to trace the current metallicity of the gaseous component of galactic discs. The preferred method is obviously through the “direct”-\(T_\text {e}\) method (see Sect. 3.1), which however, due to the weakness of the auroral lines, can be properly mapped only in the MW (Deharveng et al. 2000; Esteban et al. 2005; Rudolph et al. 2006; Balser et al. 2011) and in a few nearby galaxies (e.g., Bresolin 2007; Bresolin et al. 2009a, 2012; Werk et al. 2011; Berg et al. 2012, 2013, 2015a). The strong line method has been extended to probe much larger samples of galaxies, as discussed later on, but with the uncertainties associated with the latter method, as discussed in Sect. 3.5.

Massive stars are also an alternative, accurate method to probe the metallicity of the gas out of which these stars have recently formed. Also in this case, the difficulty of obtaining the high quality and (high resolution) spectra required to measure the abundances has limited this method to the MW and a few other nearby galaxies (e.g., Daflon and Cunha 2004; Bresolin 2007; Davies et al. 2015), see the discussion in Sect. 3.4.

Cepheids have also been extensively used as a tool to investigate the current metallicity gradient, especially in the MW, as in their case, besides being very luminous, the distance is determined with high accuracy (Luck et al. 2011, 2006; Andrievsky et al. 2004). As a result, they provide probably the most accurate metallicity gradient measurements in the Milky Way (Fig. 29). The MW radial gradient of [Fe/H] obtained by Luck et al. (2011) with this method is \(\hbox {d[Fe/H]/d}R_{\mathrm{G}} = -\,0.062 \pm 0.002~\hbox {dex}~\hbox {kpc}^{-1}\). The gradient of [O/H] has nearly the same slope within uncertainties: \(\hbox {d[O/H]/d}R_\mathrm{G} = -\,0.056 \pm 0.003~\hbox {dex}~\hbox {kpc}^{-1}\).

Fig. 29
figure 29

Image reproduced with permission from Luck et al. (2011), copyright by AAS

Metallicity gradient of the MW based on Cepheids measurements (Luck et al. 2011).

Open clusters have instead been used to probe the metallicity of the gas when they formed, typically probing ages of a few Gyr (Friel 1995; Chen et al. 2003; Sestito et al. 2006; Magrini et al. 2009; Lépine et al. 2011). By using the extensive ESO-Gaia survey Magrini et al. (2017) has provided the most extensive mapping of the metallicity of open clusters (and field stars) in the MW, by also differentiating between clusters younger and older than 2 Gyr.

Planetary Nebulae instead probe the enrichment of the gas back to the time when their progenitors formed, on timescales ranging from about 3–4 Gyr to potentially the oldest epochs (\(\sim 13\) Gyr), and they have been extensively used to investigate gradients in nearby galaxies and in the MW (e.g., Maciel and Quireza 1999; Maciel et al. 2003; Henry et al. 2010; Stanghellini and Haywood 2010; Stanghellini et al. 2014), in some works even differentiating among the ages of the planetary nebulae (Stanghellini and Haywood 2018).

In principle the different tracers discussed above enable us to determine the evolution of the chemical gradient in galaxies as a function of lookback time. The general result is that diagnostics tracing metal enrichment on longer time scales tend to give gradients that are flatter than those inferred from HII regions, suggesting that metallicity gradients have become steeper (more negative) with time (Magrini et al. 2016). This is shown in Fig. 30, from Stanghellini et al. (2014), where the evolution of the metallicity gradients is shown as a function of lookback time, for a few galaxies for which this information can be extracted.

Fig. 30
figure 30

Image reproduced with permission from Stanghellini et al. (2014), copyright by ESO

Evolution of the metallicity gradient in the MW and in a few additional galaxies as a function of lookback time, based on different metallicity tracers probing the gas phase at different epochs. This diagram shows that gradients become steeper, more negative, as time flows

This is somewhat in contrast with simple expectations of inside-out galaxy formation (inferred by other tracers), in which the inner regions would be expected to start forming at earlier times, hence having more times to produce metals, than the outer galactic regions (e.g., Davé et al. 2011; Gibson et al. 2013; Prantzos and Boissier 2000; Pilkington et al. 2012). This steepening requires some heavy redistribution of metals in the early galaxy formation, such as powerful feedback effects, prominent radial flows (for instance induced by gravitational instabilities or galaxy interactions) or early stochastic accretion/dilution (e.g. from intensive accretion and minor merging events in the early universe) (Dekel et al. 2013; Dekel and Mandelker 2014; Tissera et al. 2018; Grisoni et al. 2018).

However, one should also be aware that the tracers probing the metallicity gradients across long lookback times, such as PNe or older open clusters, are potentially affected by the effect of radial stellar migration. Indeed, stars may potentially migrate significantly from their original birth site, washing out an originally steep gradient, as a consequence of stellar bars, galaxy interactions or other secular processes. In this scenario the flatter gradients observed at later cosmic times would simply be a consequence of the longer timescale during which stellar migration has been mixing the older stellar populations. Yet, as discussed later, other results seem to independently support the scenario in which metallicity gradients have become steeper with time. Moreover, some detailed modelling of the prominence of stellar migration across the lifetime of galaxies have indicated that this effect is minor and unlikely to substantially affect the slope of the metallicity gradient in galaxies (Spitoni et al. 2015).

6.2 Statistical properties of galactic discs metallicity gradients

The advent of large integral field spectroscopic surveys has made it possible to investigate gradients more systematically for large sample of galaxies, especially exploiting HII regions, whose nebular lines are easier to detect and map in galaxies. By using spatially resolved spectroscopic data from a sample of 306 star-forming discs from the CALIFA survey, Sánchez et al. (2014) found a large spread of metallicity gradients. However, both Sánchez et al. (2014) and Ho et al. (2015) have pointed out that the spread is greatly reduced, and metallicity gradients become comparable, if the radii are normalized to the galaxy effective radius (\({R}_{\mathrm{e}}\)), suggesting that the chemical evolution of galaxies with different masses is governed by the same enrichment processes occurring on local scales.

Ho et al. (2015) and Sánchez-Menguiano et al. (2016) found no evidence that the metallicity gradient (normalized to \({R}_{25}\)) depends on galaxy mass, based on data of star-forming galaxies from the SAMI and CALIFA integral field surveys, However, the size and/or the mass range of these samples may not be sufficient to identify clear trends with mass. Indeed, by using the second MaNGA–SDSS4 data release, comprising integral field spectroscopic data for about 2800 galaxies spanning two orders of magnitudes in stellar mass, Belfiore et al. (2017) later revealed that the metallicity gradients of star-forming galaxies (normalized to \({R}_{\mathrm{e}}\)) actually strongly depend on the stellar mass of the galaxy. Indeed, as illustrated in Fig. 31, the metallicity gradient is nearly flat for low mass galaxies (\({M}_{{star}}\sim 10^9~\hbox {M}_{\odot }\)) and becomes progressively steeper (more negative) for more massive galaxies. If one considers that the sequence in mass somehow reflects an evolutionary sequence, the mass-gradient relationship is at least qualitatively in agreement with the indication that the metallicity gradient was shallower at earlier epochs as inferred by the different tracers of individual galaxies.

Fig. 31
figure 31

Image reproduced with permission from Belfiore et al. (2017), copyright by the authors

Average gaseous metallicity profile in bins of stellar mass for \(\sim \) 2800 galaxies from the MaNGA survey. Metallicities are derived from O3N2 using the calibrations in Pettini and Pagel (2004)

A steepening of the metallicity gradient with stellar mass has later been confirmed by Poetrodjojo et al. (2018) by using data from the SAMI Survey (Bryant et al. 2015). The dependence seems weaker than observed by Belfiore et al. (2017), but the sample of galaxies is also much smaller.

It should be noted that these results on metallicity gradients of large samples of galaxies, based on the extensive MaNGA, CALIFA and SAMI surveys, adopt the strong line methods for measuring the metallicity, with all caveats discussed in Sect. 3.4. However, the absolute calibration scale offset potentially plaguing the strong line method is partly mitigated in the case of metallicity gradients as these imply differential measurements. It is also important to recall, that all these studies confine the measurement of the metallicity gradient to the galactic regions showing evidence for star formation; extensive bulge or inter-arm regions with LIER-like emission or with nebular emission too weak to be probed are excluded from the determination of the metallicity gradients, which may result in either potential bias or may be missing some key information associated with the metal evolution in these specific regions.

The flattening of the metallicity gradient in the central region of the most massive spiral galaxies (Fig. 31, see also Zinchenko et al. 2016) is likely a consequence of the metallicity saturating, and approaching the yield, in the central most metal rich regions. However, it is also possible that, despite the attempt to confine the measurement of the metallicity gradient to star-forming regions, some contamination from the central LIER-like emission in massive galaxies may affect the strong line diagnostics.

Some attempts have been made to explore azimuthal variations of the metallicity (i.e., at fixed galactocentric distance) in galactic discs, both by using the direct \(T_\text {e}\) method (Li et al. 2013; Berg et al. 2015a) and the strong line method (Zinchenko et al. 2016; Sánchez-Menguiano et al. 2017a; Ho et al. 2017, 2018). Generally, the azimuthal variations, if anything, are found to be small (less than 0.1 dex), but significant in a growing number of systems, as shown for example in Fig. 32 (Ho et al. 2018). There are also claims that, in particular, metallicity varies between the arm and inter-arm regions in some galaxies (Ho et al. 2017, 2018), which may reveal fast enrichment by the enhanced star formation along the spiral arms. However, these properties do not seem to be common to most galaxies (Sánchez-Menguiano et al. 2017a). Moreover, one should take into account that the azimuthal metallicity gradients, as well as the arm/inter-arm variations, have so far been assessed primarily through strong line diagnostics and, despite the efforts to take the ionization parameter into account, the small variations observed may still be affected by changes in the physical conditions of the ISM and excitation conditions between arm and inter-arm regions. Studies of azimuthal variations based on direct (\(T_\text {e}\)) measurements are less conclusive in establishing whether there are significant azimuthal metallicity gradients in local galaxies, although the statistics of HII-regions is admittedly poor in this case. Certainly more studies are needed both to expand the samples of galaxies with angular resolution high enough to resolve azimuthal structures and sensitive enough to detect auroral lines in multiple galactic HII regions.

Fig. 32
figure 32

Image reproduced with permission from Ho et al. (2018), copyright by ESO

Left: metallicity distribution of HII regions in the galaxy NGC2997. Right-top: metallicity distribution for the same galaxy in polar projection. Right-bottom: metallicity difference with the average radial gradient. The positions of the spiral arms are indicated

In summary, the study of the radial metallicity gradients is an important piece of information on the current status of the galaxies and on the process dominating their formation. As different populations of the MW (HII regions, young massive stars, open clusters, PNe, etc.) reveal the chemical abundances at different ages and stages of the Galaxies, the evolution of chemical abundance gradient with time can be studies via stellar archeology. A steepening of the gradient with time towards more negative values is detected, and this is explained by the presence of radial redistribution of metals.

6.3 Galactic disk outskirts

While the slope of the radial gradients has generally been probed within the optical radius, high sensitivity observations have made it possible to explore the radial gradients in the outer discs. Observations have shown that at large galactocentric radii (typically at \(R>2R_\text {e}\)), metallicity gradients tend to become very flat and settle to a relatively high value, typically around 0.3–0.5 \({\hbox {Z}}_{\odot }\). This has been confirmed both through the direct method applied to individual galactic discs (Bresolin et al. 2009b, 2012; Goddard et al. 2011; Sánchez-Menguiano et al. 2017b), as illustrated in Fig. 33, and through the stacking of large sample of galaxies (Sánchez et al. 2014). It has been suggested that the flat slope of the metallicity gradient at large radii may be associated with the low efficiency and discontinuous star formation typical of outer discs, similar to what is observed in low-mass galactic discs. However, the really puzzling finding is the relatively high level of enrichment in these outer regions, where the formation of stars has been very low and certainly not enough to bring the content of metals to the observed values. The only realistic explanation is that these outer regions have accreted pre-enriched material, either as a consequence of cooling from the halo, gas stripping from the center due to galaxy interactions, minor merging with enriched satellites, galactic fountains or major outflows from the central (metal rich) regions (Bresolin et al. 2012; Sanchez et al. 2013; Belfiore et al. 2016a).

Fig. 33
figure 33

Image reproduced with permission from Bresolin et al. (2012), copyright by AAS

Schematic representation of the abundance gradients for a sample of local galaxies. Dots are shown at \(R = 0.4 R_{25}\) for each galaxy to represent the characteristic abundances of their inner disks

It is also important to mention that there are some notable exceptions to these trends. Indeed, Moran et al. (2012) have found cases of spiral galaxies whose metallicity gradient drops significantly at \({R}>{R}_{90}\) (the radius enclosing 90% of the r-band light) and they find that this feature is linked to the amount of atomic gas content HI. This may indicate that galaxies with prominent metallicity drop in their outskirts have been recently accreting pristine/low-metallicity gas from the intergalactic medium that, for momentum conservation, is predominantly deposited in the outer regions.

6.4 Spatially resolved scaling relations

Within this context there has been recently an extensive effort in trying to depart from metallicity gradient studies with the classical radial axisymmmetric (or azimuthal) approach, and investigate the variation of metallicities on local scales focusing on the potential correlation with the other local galactic properties, such as stellar surface density and star formation rate surface density. This approach is equivalent to investigating whether scaling relations apply locally.

More specifically, by using CALIFA and MaNGA data, González Delgado et al. (2014) and Barrera-Ballesteros et al. (2016) investigated the spatially resolved dependence of the gas metallicity on local properties of galaxies. They both find a clear correlation between metallicity and stellar mass surface density \(\varSigma _{*}\). Barrera-Ballesteros et al. (2016) investigate also the possible dependence of metallicity on local surface density of star formation rate, \(\varSigma _{\mathrm{SFR}}\), but they claimed that a Z–\(\varSigma _{\mathrm{SFR}}\) relation is not significant, although (at fixed \(\varSigma _{*}\)) their highest \(\varSigma _{\mathrm{SFR}}\) bin does show a clear drop in metallicity. More recently, Belfiore et al. (in prep.) have further investigated in detail the spatially resolved metallicity scaling relations by using the latest MaNGA releases; their analysis confirmed the Z–\(\varSigma _{*}\) scaling relation and revealed also the existence of a clear dependence of metallicity on SFR surface density, i.e., a Z–\(\varSigma _{\mathrm{SFR}}\) anti-correlation, as illustrated in Fig. 34 (left and center).

Fig. 34
figure 34

Image reproduced with permission from Barrera-Ballesteros et al. (2018), copyright by AAS

Left: gas metallicity as a function of stellar mass surface density, averaged across the entire MaNGA sample (black stars) and averaged in bins of star formation rate surface density (coloured lines). Center: gas metallicity as a function of star formation rate surface density, averaged in bins of stellar mass surface density. Courtesy of Belfiore et al. (in prep.). Right: gas metallicity as a function of local gas fraction \(\mu \), as inferred from the Balmer decrement

Barrera-Ballesteros et al. (2018) found also a strong anti-correlation between local metallicity and gas surface density \(\mu \) (inferred from the dust extinction traced by the Balmer decrement) as illustrated in Fig. 34, right. This finding is in good agreement with the Z\(\varSigma _{\mathrm{SFR}}\) anti-correlation, as \(\varSigma _{\mathrm{SFR}}\) and gas surface density \(\mu \) are tightly linked by the Schmidt–Kennicutt relation. It is also in agreement with the results on the gas-FMR by, e.g., Bothwell et al. (2013), see Sect. 5.2.3.

The (anti-)correlation between local metallicity and surface density of SFR (or, equivalently, on gas surface density) may really be driving the global FMR; in the scenario in which local metallicity dilution due to infalling gas both decreases the local metallicity and boosts the local SFR, the cumulative effect can scale up to induce the FMR over the entire galaxy.

Whether the global MZR stems from the local correlation between metallicity and \(\varSigma _*\), or vice-versa, is more difficult to establish. Indeed, total stellar mass and \(\varSigma _*\) are correlated, and therefore it is difficult to establish which of the two is the primary, driving correlation. Barrera-Ballesteros et al. (2016) claim that the local Z–\(\varSigma _*\) is the primary relation, which drives both the global MZR and the radial metallicity gradients in galaxies. However, this certainly cannot be the case at large galactic radii, where the metallicity profile flattens and remains at relatively high levels, while the stellar surface density keeps fading exponentially or even more rapidly. Moreover, recently D’Eugenio et al. (2018) have shown that the gas metallicity is more tightly correlated with the gravitational potential (\(\varPhi \sim {M}_{*}/{R}_{\mathrm{e}}\)) than with the galaxy stellar mass or the stellar surface density, indicating that the gravitational potential is likely the primary mechanism establishing the level of metal content (see also Sect. 5.3). The weak correlation found by Barrera-Ballesteros et al. (2018) between metallicity and local escape velocity suggests that metals lost by winds only play a minor role in shaping the local metallicity. In conclusions, these evidences show that the MZR might be shaped by global rather than local effects.

Summarizing, the scaling relations observed at the global, galaxy-scale level are now also observed at the local level. While the global FMR is probably driven by the local one, the situation for the MZR is uncertain because conflicting evidences are present.

6.5 Interacting galaxies

While most of these studies have focused on regular, isolated galactic discs, a few studies have investigated gradients in interacting/merging systems. It has been found that interacting systems generally have significantly flatter gradients than isolated galaxies, and the effect is stronger in systems that are in a more advanced stage of merging (Fig. 35) (Rupke et al. 2010b; Rich et al. 2012a; Torres-Flores et al. 2014). Moreover, the extended tails resulting from galaxy interaction display remarkably flat gradients out to 70 kpc (Olave-Rojas et al. 2015). Recently Ellison et al. (2018b) has used MaNGA data to study the spatial distribution of the excess of SFR and deficit of metallicity in interacting/merging galaxies, finding that both are more prominent in the inner parts of the galaxies.

Fig. 35
figure 35

Image reproduced with permission from Rich et al. (2012a), copyright by AAS

Gas phase metallicity gradient versus merger stage. Blue points are isolated systems from Rupke et al. (2010b), yellow points are wide pairs, red points are luminous infrared galaxies and black points are from models (Rich et al. 2012a).

The most widely accepted interpretation is that interactions on the one hand make the outer low-metallicity gas lose angular momentum and flow towards the central region of the galaxy causing dilution, while on the other hand metal enriched gas is stripped into extended tails where the metallicity imprint of the original location is rapidly lost and mixed up with gas from other regions (Torrey et al. 2012; Rupke et al. 2010a). This interpretation is consistent with the lower total metallicity observed in these galaxies as discussed in Sect. 5.1.7, and with the observations that recently merged galaxies have larger amounts of atomic gas (Ellison et al. 2018a).

6.6 Stellar metallicity gradients

Several studies have been performed on the metallicity gradients of the stellar population in early type galaxies (e.g., Spolaor et al. 2010; Bedregal et al. 2011; Koleva et al. 2011; Harrison et al. 2011; La Barbera et al. 2012; Goddard et al. 2017). They are often more demanding than gas metallicity gradients, as resolved stellar metallicities require much higher S/N data on the stellar continuum, which is especially difficult to achieve in the outer parts of early type galaxies galaxy due to the rapidly declining surface brightness of the stellar light with galactocentric radius. However, given typically the lack of dust reddening, stellar colours have also been employed to map the metallicity in these passive systems (e.g., Tortora et al. 2010). Fewer studies have been performed on the stellar metallicity gradients of star-forming galaxies (e.g., Morelli et al. 2015; Sánchez-Blázquez et al. 2014; González Delgado et al. 2015; Goddard et al. 2017; Li et al. 2018) both because disentangling age-metallicity degeneracies of the stellar population requires even higher S/N and because the presence of nebular lines makes the analysis of the stellar features more difficult. Moreover, the simple analysis of the stellar metallicity provides light-weighted stellar metallicities, typically dominated by the most recent stellar population. A proper determination of the dominating stellar population, especially in late type galaxies is difficult as it depends on the star formation history and on the stellar indices or wavelength range used to extract the stellar metallicity. Mass-weighted stellar metallicities are potentially more interesting to investigate the metalliFcity gradient of the bulk of the stellar population, but, as discussed in Sect. 2, mass-weighted stellar metallicities are more difficult to infer, especially in the outer regions of galaxies where the S/N is not as high as in the central galactic regions.

Typically, stellar metallicity gradients are shallow. Although the dispersion is large, works based both on CALIFA (González Delgado et al. 2015) and on MaNGA (Goddard et al. 2017) find a significant dependence of the stellar metallicity gradient with stellar mass, with a stronger dependence for late type galaxies than for early type galaxies (Fig. 36). This is a trend similar to that observed for the gas phase metallicity, however the important difference is that for stellar metallicities the gradient is negative but shallower at high masses and becomes positive (“inverted”) at low masses, indicating a change in the formation and accretion histories between low mass and high mass galaxies, and also between late type and early type galaxies. A strong dependence of stellar metallicity gradients with stellar mass for late type galaxies was also confirmed by Lian et al. (2018); they point out that the mass dependence is significantly stronger than for the gaseous metallicity gradients, which is surprising, given that close to equilibrium the stellar metallicity should approach the gas metallicity (Peng and Maiolino 2014a); they interpret the strong difference between stellar and ISM metallicity gradients invoking a variation of either the outflow loading factor, or of the IMF, both in time and radially.

Fig. 36
figure 36

Image reproduced with permission from Goddard et al. (2017), copyright by the authors

Mass-weighted stellar metallicity gradient as a function of stellar mass for early-type (left) and late-type (right) local galaxies from the MaNGA survey

Interestingly, by using MaNGA data, Li et al. (2018) have found that the stellar metallicity gradients depend on stellar velocity dispersion and that they peak (becoming most negative) at intermediate velocity dispersions of about 100 km s\(^{-1}\). This is interpreted as indicating a change in the evolutionary history in galaxies. In particular, metallicity gradients becoming flat at very large velocity dispersions is likely to indicate a growing role of mergers, which redistribute the metallicity in very massive galaxies whose velocity dispersion has been enhanced by the merging history.

Finally, very interestingly, Sánchez-Blázquez et al. (2014) have investigated the relation of stellar metallicity gradients on the presence and strength of stellar bars, finding no evidence for any correlation (in agreement with similar studies tracing the gas metallicity gradients, see Sánchez et al. 2014). This result indicates that stellar migration associated with stellar bars does not play a significant role in shaping the metallicity gradient of galaxies, in agreement with the expectations of some models (Spitoni et al. 2015). However, one should take into account that stellar bars are a recurrent phenomenon in the life of galaxies (about 40% of spiral galaxies in the local universe are barred, Sheth et al. 2008), therefore the lack of a correlation with the current strength of bar in galaxies does not necessarily mean that past barred phases (even in galaxies currently non-barred) have not played a role.

6.7 Metallicity gradients at high redshift

Measuring metallicity gradients at high redshift is obviously much more difficult for various reasons.

Firstly, the steep cosmological dimming of the surface brightness (\(\propto (1+z)^4\)) makes it much more difficult to achieve the S/N required to measure metallicity gradients in the outer parts of galaxies; this effectively translates into measuring high-z metallicity gradients only for the gas phase and only by using strong line methods.

Secondly, it becomes increasingly difficult to spatially resolve galaxies in integral field spectroscopy. In most studies the metallicity gradients are only marginally resolved. Targeting lensed galaxies generally helps (e.g., Yuan et al. 2011; Jones et al. 2013, 2015b; Leethochawalit et al. 2016; Wang et al. 2017), however at the cost of introducing the additional uncertainty associated with the lens modelling. Moreover, even when the lens model is well constrained, gravitational magnification is differential, hence the resulting lensed image and metallicity map are strongly weighted towards the regions close to the lens caustic, therefore potentially resulting into distorted metallicity maps. The use of adaptive optics certainly helps (e.g., Leethochawalit et al. 2016; Perna et al. 2018; Förster Schreiber et al. 2018) although the modest Strehl-ratios cause the sensitivity to drop significantly, especially towards the outer low surface brightness regions. HST grism spectroscopy has also been effectively used to map the metallicities at high-z with HST-like high angular resolution (Fig. 37) (Jones et al. 2015b; Wang et al. 2017), the only problem being the low spectral resolution of the spectra and small wavelength range, which often limit the use of diagnostics and makes the subtraction of the stellar continuum more problematic. In other studies, not exploiting gravitational lensing and from the ground without adaptive optics (e.g., Cresci et al. 2010; Queyrel et al. 2012; Swinbank et al. 2012; Stott et al. 2014; Troncoso et al. 2014; Wuyts et al. 2016), the measurement of metallicity gradients has generally been limited to the larger (hence typically more massive) galaxies therefore with potential bias.

Fig. 37
figure 37

Image reproduced with permission from Wang et al. (2017), copyright by AAS

Metallicity gradient in a lensed galaxy at \(z=1.48\) obtained with HST grism data. Left: metallicity map. Center: uncertainty. Right: radial metallicity distribution (blue points are individual pixels, red points show the average within annuli)

The additional problem is that most nebular diagnostics are shifted to the near infrared bands. The basic, bluest diagnostics enabling the use of the strong line method (e.g., the \({R}_{23}\) parameter) can be traced in the optical band through optical integral field spectrometers only out to \(z\sim 0.8\) (Carton et al. 2018). While integral field spectrometers in the near-IR bands are well developed and available on large telescopes to probe the strong line metallicity diagnostics at \(z>1\), the reduced sensitivity in these bands (primarily because of the higher background, both thermal and from bright OH sky lines) and discontinuous spectral coverage (because of the deep atmospheric absorption bands), makes the measurement of metallicity gradients even more challenging. Often, spectra are obtained in a single spectral band, resulting in limited information. For instance, [NII]/H\(\alpha \) is used to effectively trace metallicity gradients at high-z, as these two lines are conveniently observed in the same band (e.g., Wuyts et al. 2016; Förster Schreiber et al. 2018), but at the cost of introducing all uncertainties associated with this single diagnostic (e.g., dependence on nitrogen enrichment, dependence on ionization parameter, etc), see Sect. 3.5.4.

With all these caveats, extensive studies have been undertaken to constrain the evolution of metallicity gradients at high redshift (e.g., Cresci et al. 2010; Queyrel et al. 2012; Swinbank et al. 2012; Stott et al. 2014; Troncoso et al. 2014; Wuyts et al. 2016; Yuan et al. 2011; Jones et al. 2013, 2015b; Leethochawalit et al. 2016; Wang et al. 2017; Carton et al. 2018; Förster Schreiber et al. 2018). A summary of the observed evolution of the metallicity gradients is shown in Fig. 38 (from new observations and a compilation courtesy of Curti et al., in prep.). There is clearly a large dispersion, part of which is likely due to the observational uncertainties, but it is also likely reflecting a real dispersion of the metallicity gradients during the early phases of galaxy evolution, when the accretion and merging processes were more stochastic and resulting into a more irregular behaviour (e.g., Ceverino et al. 2016). Despite the large scatter, and with the exception of a few rare cases of very steep gradients (Jones et al. 2013, based on low S/N spectra), most studies find that at high redshift the metallicity gradients are on average flatter than observed locally. This is in agreement with the flattening of the gradients in local galaxies when using tracers that probe longer lookback times, i.e., primarily PNe (Stanghellini et al. 2014; Stanghellini and Haywood 2018), as discussed in Sect. 6.1.

Fig. 38
figure 38

Image courtesy of Curti et al. (in prep.)

Overview of the metallicity gradients measured at different redshifts, both in lensed and unlensed galaxies, from different surveys, including local galaxies (whose gradient evolution has been inferred through tracers at different lookback times). Evolutionary tracks of different models are also show (see legend)

It is interesting to note that studies comparing metallicities measured in DLA absorption systems with their optical counterparts (whose metallicity is measured through nebular lines) have also independently inferred gradients that are, on average, quite flat (\(-\,0.022\ {\mathrm{dex}}\ {\mathrm{kpc}}^{-1}\)) even out to galactocentric radii of several tens kpc, in high-z galaxies (Christensen et al. 2014; Rhodin et al. 2018).

However, it is becoming increasingly evident that at such early epochs the concept of “radial” (azimuthally averaged) gradient loses part of its meaning, as expected by some models (Ceverino et al. 2016). Indeed, the metallicity distribution is often very irregular, with large local variations (e.g., Förster Schreiber et al. 2018, and Fig. 37), hence the apparently flat radial gradient simply results from averaging large azimuthal variations into the same radial bins. These large irregular metallicity variations in high-z galaxies probably reflect the chaotic accretion and formation processes during these early phases. Therefore, the comparison with models should not be done simply in terms of radial gradients (which can be partly deceiving of the chemical complexity of these systems), additional information should be considered such as the metallicity scatter.

A population of positive (i.e., “inverted”) metallicity gradients has recently been discovered at high redshift (Cresci et al. 2010; Troncoso et al. 2014; Carton et al. 2018) (Wang et al. 2018). These inverted gradients appear to be particularly common at \(z>3\) (Fig. 39). However, it has been pointed out that rather than simply being inverted, the central metallicity suppression is mostly associated with the central enhancement of star formation rate. Cresci et al. (2010) and Troncoso et al. (2014) suggest that this is a consequence of enhanced inflow of metal-poor gas toward the central galaxies, in their early epoch of formation, which results in both a local dilution of the metallicity and a local enhancement of the SFR due to the larger amount of gas available. Stott et al. (2014) suggest a similar scenario to interpret the correlation that they find at \(z\sim 1\) between metallicity gradient and sSFR: the enhanced inflow of near-pristine gas makes the central region metal poor and also boosts the sSFR. This is in agreement with the infall interpretation of the FMR, see Sect. 5.2.2.

Fig. 39
figure 39

Image reproduced with permission from Cresci et al. (2010), copyright by Macmillan

Example of galaxy with inverted metallicity gradient at \(z=3.06\). Left: map of the [OIII]5007 line emission (approximately proportional to the surface density of SFR). Center: velocity field. Right: metallicity map

Whether the metallicity gradients of high-z galaxies have a mass dependence similar to that observed locally is not totally clear, mostly because of the large scatter. Wuyts et al. (2016) and Stott et al. (2014) do not find a clear mass dependence in their samples at \(z\sim 1\)–2. Carton et al. (2018) find a weak mass dependence in their sample at \(0.1<z<0.8\). Together with the local findings these results suggest that the mass dependence of the metallicity gradients is established at late epochs. Correlations of the gradients with sSFR have been also proposed (Stott et al. 2014; Wuyts et al. 2016), but in all cases only with low significance and in disagreement with other studies (Leethochawalit et al. 2016; Troncoso et al. 2014).

Finally, Carton et al. (2018) find an interesting correlation at \(0.1<z<0.8\) between metallicity gradient and galaxy size, in which small galaxies show a large spread in metallicity gradient (both positive and negative), while large galaxies present a much smaller spread and regular negative gradients.

6.8 Metallicity gradients models

Several models have been proposed to reproduce the metallicity gradients in galaxies, both for the gas and the stars (Mollá et al. 1997; Chiappini et al. 2001; Naab and Ostriker 2006; Mott et al. 2013; Spitoni et al. 2013, 2015; Ho et al. 2015; Kudritzki et al. 2015; Ascasibar et al. 2015; Schönrich and McMillan 2017; Lian et al. 2018). Many models consist in the extension of the analytical “gas regulator” scenario applied to concentric galactic rings, in which there is a radial dependence of the primary parameters, such as gas infall timescale and star formation efficiency (also introducing threshold for star formation). Some models, as discussed above, also introduce variable outflow loading factors and even variable IMF. The main challenge of these models is that they have also to reproduce the inside-out growth of galaxies, i.e., the finding that, based on both the stellar population age gradients and surface density gradients, the central parts of galaxies must have grown earlier and faster than the outer parts. This requires models to have accelerated star formation and enrichment in the central regions, relative to the outer regions, which results into a negative metallicity gradient, as observed in galaxies.

The main problem of this scenario is that, if different galactic annular rings do not exchange metals (as it is the case in most models), then the unavoidable consequence is that the metallicity gradient should flatten with cosmic time, which is opposite to that observed as inferred from the observation in local galaxies diagnostics that trace metallicity gradients at earlier epochs (Fig. 30), from the observed evolution of metallicity gradients at high redshift (Fig. 38) and even based on the simple finding that the metallicity gradient steepens with galaxy stellar mass (Fig. 31) and assuming the stellar mass sequence of galaxies somehow reflects their evolutionary pattern. Possible solutions to comply with the inside-out growth of galaxies and having metallicity gradients that steepen with time is that stellar migration plays a role in mixing stellar metallicities and the production of metal at different radii, from different generations of stars (Spitoni et al. 2015; Schönrich and McMillan 2017), or prominent radial flows of gas that can dilute the central regions at early epochs (Spitoni and Matteucci 2011; Mott et al. 2013; Spitoni et al. 2013), or strong feedback (in the form of outflows) at early epochs has redistributed the metals produced by the central active region in the circumgalactic medium and towards the external region of galactic discs. This is not unreasonable given that recent observations have revealed that as much as 40% of the metals produced by galaxies have been expelled in their halo (see Sect. 9) and that the circumgalactic medium has already been significantly enriched (\({Z}\sim 0.1~{\hbox {Z}}_{{\odot }}\)) by \(z\sim 2\) (Prochaska et al. 2013). Figure 40 shows the effect of introducing radial flow in the analytical models proposed by Mott et al. (2013). Interestingly, in this models the high inflow rate of pristine gas towards the central regions even inverts the gradient centrally, nicely reproducing the inverted gradients observed at high redshift (Fig. 39Cresci et al. 2010; Troncoso et al. 2014).

Fig. 40
figure 40

Image reproduced with permission from Mott et al. (2013), copyright by the authors

Examples of metallicity gradients predicted by the analytical model of Mott et al. (2013), which include radial flows of gas. The model on the left assumes a threshold, inside-out formation, constant star formation efficiency, and radial flows. The model on the right uses a radially variable star formation efficiency. It has to be noted that the metallicity gradient becomes more negative with time (at least within the central 10 kpc, which is the only part observed at high redshift) and also that in the left-side model at early times the metallicity gradient is centrally inverted as a consequence of strong central infall of pristine gas

Zoom-in cosmological simulations have become increasingly popular to investigate the evolution of metallicity gradients in galaxies (Di Matteo et al. 2009; Torrey et al. 2012; Pilkington et al. 2012; Gibson et al. 2013; Tissera et al. 2016a, b, 2017, 2018). However, these high-resolution simulations are computationally expensive. As a consequence, only a limited number of galaxies are generally simulated at this level of detail. Simulations typically obtain reasonably metallicity gradients, although some of them struggle to reproduce the observed mass dependence. A common outcome is that galaxy interactions tend to flatten the metallicity gradients at any epoch (Di Matteo et al. 2009; Ma et al. 2017). There is typically a large scatter and rapid variations in the gradients observed in simulations especially at high redshift; Ma et al. (2017) suggest that the observed gaseous metallicity gradients reflect more the current status of the galaxy rather than tracing its past cosmological evolution. Many simulations tend to have the same problem as analytical models in reproducing the temporal steepening of metallicity in regular, rotating galaxies. Within this context the most successful models are those that introduce enhanced feedback from star formation, which redistributes the metals across the galaxy more effectively in the early, active stages of galaxy formation (Gibson et al. 2013).

6.9 Summary of metallicity gradients in galaxies

In this section we summarize some of the key results regarding the spatially resolved distribution of metals in galaxies.

Radial metallicity gradients tend to steepen in more massive galaxies. This is observed both for the gas phase metallicity and for the stellar metallicity (although with larger scatter).

However, metallicity gradients are not described by a single slope; they generally flatten at high galactocentric radii (suggestive that galactic outskirts have accreted pre-enriched material) and, in massive galaxies, they tend to flatten also in the central region (likely because of metallicity saturation).

There is evidence that the metallicity is not simply a function of galactocentric radius, but it depends on the local physical properties of the galactic disc. More specifically, the gas metallicity is found to correlate with the surface density of stellar mass and to anti-correlate with the surface density of star formation rate and with the surface density of gas. Whether these local scaling relations are responsible for the global metallicity scaling relations has yet to be properly assessed.

The use of indicators that trace the metallicity in galaxies at different lookback times, as well as the direct observation of metallicity gradients at high redshift, indicate that metallicity gradients were flatter in the past (which is also consistent with the local steepening of metallicity gradients as a function of stellar mass, if one consider low mass galaxies progenitors of more massive ones). The fact that radial metallicity gradients steepen with cosmic time is in contrast with the simplest expectations of the chemical enrichment in the scenarios in which galaxies grow inside-out (as confirmed by independent observations), which would expect a flattening of the metallicity gradients with time. Solutions to this issue involve either stellar radial migration, or prominent early radial inflows of low-metallicity gas or feedback (outflow) effects that redistribute the metals produced in the central regions in the circumgalactic medium and towards the outer galactic regions.

It is finally interesting to note the growing evidence for a population of high-z galaxies with inverted (i.e., positive) radial gradients, which may be tracing systems in which prominent inflows of low metallicity gas are taking place or which are associated with galaxy merging/interaction which are effective in driving metal-poor gas from the outskirts towards the central region, hence flattening or even inverting metallicity gradients.

7 Relative chemical abundances

Since different elements are produced by different classes of stars/SNe and released on different timescales (see Fig. 1), the relative abundances between different chemical elements provide precious information on the star formation history and on the IMF.

The ratio between \(\alpha \)-elements, which are primarily produced on short timescales by massive stars through core-collapse SNe, and iron-peak elements, which are primarily produced by SN Ia on longer timescales, is a classical example of tracer of star formation history. Stellar populations characterized by “enhanced” \(\alpha \)/Fe must have formed on short timescale, before that SN Ia had time to enrich the ISM, while stellar populations characterized by “low” or solar-like \(\alpha \)/Fe must have formed over a prolonged phase of star formation. Other chemical elemental ratios, including CN/Fe (Carretero et al. 2004, 2007) have a similar potential of constraining the star formation history.

The detailed abundance pattern of several chemical elements can potentially provide the fingerprints of the specific stellar progenitors that have been responsible for enriching the ISM. Potentially they can even provide the signature of the enrichment by the first population of stars (Frebel and Norris 2015; Caffau et al. 2011).

In this section, after a rapid overview of the metal abundances observed in the Milky Way, we will focus primarily on the chemical abundances observed across the galaxy populations, locally and at high redshift. A detailed analysis of all chemical elements would require a full, dedicated review. We will, therefore, focus on the analysis of the most commonly used chemical elements and, in particular, \(\alpha \)/Fe, N/O and C/O, although we will also mention some other chemical abundance ratios.

7.1 The Milky Way

Extensive spectroscopic surveys have enabled astronomers to map the relative chemical elements of large numbers of stars in the Galaxy (e.g., SEGUE, RAVE, APOGEE, GAIA-ESO, HERMES-GALAH, Yanny et al. 2009; Steinmetz et al. 2006; Gilmore et al. 2012; Freeman 2012; Majewski et al. 2017). High spectral resolution observations have provided unprecedented constraints on the relative abundances of most chemical elements, although on smaller samples of stars. (e.g., Reddy et al. 2006; Bensby et al. 2010; Nissen and Schuster 2010; Johnson et al. 2014; Zoccali et al. 2017, and references therein). The advent of high resolution near-IR spectroscopic surveys has further enabled astronomers to trace the chemical abundances in the inner bulge and for metal-rich cool stars, whose heavy metal lines blending and blanking at optical wavelengths makes it difficult to measure chemical abundances with classical optical spectroscopy (e.g., Rich et al. 2012b; Önehag et al. 2012; Lindgren et al. 2016).

An overview of the [\(\alpha \)/Fe] versus [Fe/H] for the Galactic disc is given in Fig. 41, obtained with APOGEE near-IR, medium-resolution spectroscopic data of \(\sim \) 70 000 red giants (Hayden et al. 2015). \([\alpha \)/Fe]–[Fe/H] diagrams are very useful in identifying stellar populations resulting from different star formation histories. Indeed, as mentioned above, the [\(\alpha \)/Fe] ratio works like a clock of the star formation history, indicating how rapidly star formation has occurred, while the [Fe/H] ratio distributes stellar populations along their temporal evolutionary sequence (of course, depending on the star formation efficiency stellar populations spread more or less quickly along this axis). In Fig. 41 different panels conveniently show the [\(\alpha \)/Fe]–[Fe/H] distribution in bins of galactocentric radius and distance from the Galactic plane, illustrating that the stellar populations in the disc have a bimodal distribution, which has been identified with the thin and thick disc. The thin disc, which dominates the stellar population on the disc plane and at intermediate and large radii, is characterized by a rather flat [\(\alpha \)/Fe]–[Fe/H] distribution indicative of a prolonged star formation history, on timescales of about 7 Gyr, probably also associated with normal star formation efficiency (typical of disc galaxies) (Matteucci and Greggio 1986; Ryde et al. 2016; Micali et al. 2013). The thick disc population, dominating at higher vertical distances from the Galactic plane and in the inner disc, is clearly \(\alpha \)-enhanced (and, on average, more metal poor than the thin disc), indicating that the thick disc was formed faster, on a timescale of about \(\sim 2\) Gyr, and likely with higher star formation efficiency. It is important to highlight that the existence of thick-thin disc bimodality is still debated. In particular, Bovy et al. (2012) claim that there is not a real bimodality but a smooth, continuous distribution if one considers the mass-weighted scale-height distribution of stellar populations.

Fig. 41
figure 41

Image reproduced with permission from Hayden et al. (2015), copyright by AAS

[\(\alpha \)/Fe] versus [Fe/H] in the Galactic disc stars in bins of galactocentric distance and vertical distance from the Galactic plane from a sample of \(\sim \) 70 000 red giants observed by the SDSS-III/APOGEE survey

There has been growing evidence that the bulge stellar population is bimodal in metallicity (e.g., Rojas-Arriagada et al. 2017), hosting both a sub-solar population and a supersolar population. The two populations also have different [\(\alpha \)/Fe] enrichment levels, as illustrated in Fig. 42, top. The low-metallicity component chemical properties are consistent with those of the thick disc. The high-metallicity component appeared similar to the thick-disk based on the optical spectroscopic surveys, but when measured through near-IR data (i.e., consistent with the disc data shown in Fig. 41), more adequate to probe the high metallicity component of the bulge, also the high metallicity component appears clearly \(\alpha \)-enhanced relative to both the thick and thin discs, as illustrated in Fig. 42, bottom (Schultheis et al. 2017). The emerging picture is that the high-metallicity component was formed through the same secular process as the thick disc, and actually associated with the inner Galactic bar, while the low-metallicity, \(\alpha \)-enhanced component was formed quickly (within less than about 0.5 Gyr), with high efficiency (which enabled quick enrichment), as a consequence of the initial gravitational collapse of the galaxy.

Fig. 42
figure 42

Images reproduced with permission, copyright by ESO

[\(\alpha \)/Fe] versus [Fe/H] in the Galactic bulge from the Gaia-ESO (optical) survey (top, Rojas-Arriagada et al. 2017) and from the APOGEE (near-IR) survey (bottom, Schultheis et al. 2017)

The halo is also \(\alpha \)-enhanced but characterized by even more metal-poor stars (Cayrelm et al. 2004; Frebel and Norris 2015) and it has been suggested that its chemical properties are also consistent with the early, fast and efficient collapse of the galaxy (Micali et al. 2013), although more recent data point at two stellar populations, differentiated both chemically and kinematically (Carollo et al. 2007; Fernández-Alvar et al. 2018), corresponding to faster and slower star formation histories, respectively, the former possibly also associated with a top-heavier IMF.

Overall, among analytical models, the abundances and metallicities in the various components of the Milky Way can be explained well in the framework of the so-called two- or three-phase infall model (Fig. 43, Micali et al. 2013; Chiappini and Gratton 1997; Chiappini et al. 2001, 2005) in which the halo and old bulge components have formed very quickly (within less than 0.5 Myr) and very efficiently (above the Schmidt–Kennicutt relation): the thick disc (and lower metallicity component of the bulge) has formed in a second infall event (on intermediate timescales of \(\sim 2\) Gyr) and with even higher efficiency, while the thin disk has formed on much longer timescale (\(\sim 7\) Gyr), following the Schmidt–Kennicutt relation for normal star-forming discs. Both the thick and (especially) the thin disc must also have had a radially variable inflow rate. Alternative models by Schönrich and Binney (2009a, b) reproduce the multiple components as the effect of different star formation conditions in different parts of the disc followed by radial mixing of stars.

Fig. 43
figure 43

Image reproduced with permission from Micali et al. (2013), copyright by the authors

Galaxy’s three-phases infall model (3IM). Left: temporal evolution of the SFR as predicted by the 3IM. The (red) solid portion of the curve refers to the halo phase, the (black) dot-dashed one to the thick-disc phase, the (blue) dashed one to the thin-disc phase. The coloured bands are due to a large number of subsequent, rapid variations. Right: SFR density versus gas density in the halo (red crosses), thick disc (black crosses) and thin disc (blue crosses). Also shown are the fit to the \(\varSigma _{\mathrm{SFR}}-\varSigma _{\mathrm{gas}}\) relation for local spirals and \(z=1.5\) BzK galaxies (grey solid line), the extrapolation of the starburst sequence from the same authors (grey dashed line) and the region of the plot occupied by spiral galaxies data (delimited by the dotted grey lines)

Summarizing, different morphological components of the Milky are characterized by different chemical abundance patterns, especially for what concerns \(\alpha \)/Fe, indicating different formation histories and different formation processes. More specifically, halo, bulge, thick and thin discs are characterized by gradually later ages of star formation, on gradually longer timescales, and likely associated with different star formation efficiencies. Whether these were physically distinct phases (yielding to distinct, different populations) or part of a more smoother, continuous evolution is not totally clear.

7.2 Local galaxies

7.2.1 \(\alpha \)/Fe

As introduced above, [\(\alpha \)/Fe] is sensitive to the number ratio of CC to thermonuclear SNe and, therefore, is a powerful way to study the star formation timescales of galaxies. It is also sensitive to a number of parameters such as the IMF, the assumed stellar yields, the delay time distribution of type Ia SNe, and the differential ejection of metals into the CGM. In simple, close-box, constant-IMF models, [\(\alpha \)/Fe] is expected to evolve from a high value when, at early times and low metallicities, (little) iron production is dominated by CC events, toward a lower value, when iron produced by Ia events becomes dominant. The faster star formation occurs, the higher is the enrichment of \(\alpha \)-elements by CC SNe before that SNIa start polluting the ISM with iron.

Beyond the MW, chemical abundances of individual stars (giants/supergiants) has been determined only for a few galaxies of the Local Group, mostly dwarf satellites of the MW and Andromeda (e.g., Bonifacio et al. 2004; Monaco et al. 2005; Sbordone et al. 2007; Tolstoy et al. 2009; Cohen and Huang 2010; Kirby et al. 2011; Hill and DART Collaboration 2012; Starkenburg et al. 2013; Hendricks et al. 2014). Dwarf spheroidal galaxies are characterized by a distribution on the \(\alpha \)/Fe vs Fe/H diagram that is below the plateau observed in MW halo and thick disc stars, joining the \(\alpha \)/Fe abundance ratio of MW halo stars only at [Fe/H] \(<-1.5\). The presence of a knee in the distribution is still debated (Tolstoy et al. 2009; Kirby et al. 2011; Hendricks et al. 2014). This distribution implies that SNIa have contributed to the enrichment of the stars in these systems at most metallicities, i.e., during most of the formation of these systems, except for the their earliest phases, suggesting a bursty evolution possibly resulting from a sequence of minor accretion or merging events. Dwarf irregulars (e.g., the Small Magellanic Cloud), which are still in the process of actively forming stars, are also characterized by low \(\alpha \)/Fe abundance ratio, indicating that they have been forming stars slowly, stochastically and/or inefficiently (Matteucci and Chiosi 1983; Recchi et al. 2001). As a consequence of their shallow gravitational potential, in dwarf galaxies supernova-driven winds are also expected to play a more important role, with respect to more massive galaxies, in regulating the chemical enrichment history, by removing metals (likely in a differential way, i.e., preferentially \(\alpha \)-elements), by reducing the efficiency of star formation and by contributing to its stochasticity.

Except for the few galaxies in the Local Group, the bulk of the investigation of the chemical abundances of the stellar population in galaxies has been based on spatially integrated spectra, unavoidably implying larger uncertainties and degeneracies. However, despite these caveats, large spectroscopic surveys have enabled us to investigate the relative chemical abundances across a broad range of galaxy masses and environments.

Initial works have investigated chemical abundances using primarily the Lick indices and focusing on early type galaxies to constrain their evolutionary history. Such early works already identified that early-type galaxies have enhanced \(\alpha \) elements compared to the abundance patterns of the stars in the Galactic disc (Worthey 1994; Thomas et al. 1999). It was further found that such \(\alpha \)-enhancement increases steadily as a function of stellar velocity dispersion, which is a tracer of galaxy mass (Trager et al. 2000a; Thomas et al. 2005). Such a trend has been clearly confirmed by the full-fitting (i.e., not limited to the Lick indices) of the spectra of local galaxies (Walcher et al. 2009; Conroy and van Dokkum 2012) as illustrated in Fig. 44, where the abundance of different chemical elements relative to iron is obtained from Sloan spectra stacked in bins of velocity dispersion (Conroy et al. 2014). This trend indicates that more massive galaxies have formed much more rapidly than low-mass galaxies. More specifically, translating the velocity dispersion into galaxy mass, and combining the \(\alpha \)-enhancement information with the age of the stellar population (more massive galaxies are typically older), results into the scenario originally proposed by Matteucci and Tornambe (1987) and Matteucci (1994), further developed by later studies, and summarized in Fig. 45 from Thomas et al. (2010) (see also Thomas et al. 2005), in which more massive galaxies formed at earlier cosmic epochs (a phenomenon often referred to as “cosmic downsizing”), on shorter timescales and (based on models) more efficiently. This scenario, at least in terms of timeline sequence, has been verified through the evolution of the mass function of galaxies at high redshift, illustrating that most massive galaxies were already in place at early cosmic epochs, while lower mass galaxies have evolved more slowly (e.g., Gavazzi and Scodeggio 1996; Cowie et al. 1996; Pérez-González et al. 2008; Muzzin et al. 2013; Santini et al. 2015). Theoretical models and numerical simulations explain this phenomenon in terms of accelerated evolution in the overdense regions of the Universe, where baryons collapse more rapidly in the deepest gravitational potential wells of dark matter. The enhanced star formation efficiency in these dense regions facilitates the rapid formation of stars and rapid gas consumption. Moreover, the strong negative feedback from the resulting supernovae and rapid black-hole accretion (releasing large amount of energy through the luminous quasar phase) result in the rapid quenching of the star formation (Matteucci and Tornambe 1987; Matteucci 1994; Pipino et al. 2011; Segers et al. 2016; De Lucia et al. 2017). Additional theoretical explanations involve the effect of varying IMFs (Fontanot et al. 2017).

Fig. 44
figure 44

Image reproduced with permission from Conroy et al. (2014), copyright by AAS

Elemental abundances relative to iron inferred from the stacked spectra of local early-type galaxies in bins of velocity dispersion (hence dynamical mass), from \(\sigma =90\) km s\(^{-1}\) (blue) to \(\sigma =270\) km s\(^{-1}\) (red)

Early-type galaxies show rather flat radial gradients in terms of [\(\alpha \)/Fe], or slightly positive (e.g., Greene et al. 2013; Roig et al. 2015). If confirmed, the positive gradients of \(\alpha \)/Fe can be explained through models in which the quenching in the outer parts of massive galaxies results primarily through explosion of supernovae, whose cumulative injection of energy is more effective in ejecting gas in the shallower gravitational potential of the galaxy outskirts (Pipino et al. 2008, 2010). However, such outside-in quenching effect would be in contrast with numerical simulations that expect massive ellipticals to grow inside-out through a sequence of minor (dry) mergers from \(z\sim 2\) (Naab et al. 2007, 2009), also supported by observed size growth of elliptical as a function of reshift (van Dokkum et al. 2010), although the latter has also been interpreted as an observational effect in terms of “progenitors bias” (Lilly and Carollo 2016). However, if the positive radial gradients of \(\alpha \)/Fe in massive elliptical galaxies are confirmed with high significance, then this would be problematic to explain in the minor merger scenarios for the size growth of ellipticals.

Reproducing the \(\alpha \)/Fe enhancement in massive galaxies requires processes that enable the rapid production of stars and that then quench star formation on relatively short timescales (\(\sim \) 0.5 to 1 Gyr). In analytical models this is achieved by requiring a high star formation efficiency (which makes the formation of stars and enrichment of \(\alpha \) elements faster) and then a quenching effect that suppresses star formation (either through SNe or AGN feedback, e.g., Matteucci and Tornambe 1987; Matteucci et al. 1998; Romano et al. 2002; Pipino and Matteucci 2004; Pipino et al. 2008). In cosmological simulations the introduction of quasar feedback seems to reduce the lifetime of massive galaxies enough to reproduce the relationship between \(\alpha \)/Fe enhancement and galaxy mass (Segers et al. 2016). However, more recently, De Lucia et al. (2017) have pointed out that AGN feedback alone may not be able to simultaneously reproduce the \(\alpha \)/Fe enhancement and other galactic properties, leaving the problem open. Within this context, we note that it is not necessary to invoke the “ejective” mode of quasar feedback (i.e., removal of gas through massive quasar-driven outflows) in order to achieve a rapid quenching of star formation; indeed, a scenario in which the galaxy is simply“starved” (e.g., because its surrounding halo has been heated) would also result in a rapid cessation of star formation, as in distant powerful starburst galaxies (such as the Submillimetre Galaxies, which are often regarded as the progenitors of local massive ellipticals) the gas depletion times (by the simple effect of the highly efficient star formation) are as short as a few hundred million years.

X-ray spectroscopy offers the possibility of measuring the abundance of several chemical elements, including iron and \(\alpha \)-elements, of hot plasmas (\(>10^6\) K). Therefore, X-ray spectroscopy is extremely important to investigate, for instance, the ISM and CGM heated by SNe and galactic winds, as well as the hot gas in clusters and groups of galaxies, although sensitivity, spectral and angular resolution issues have often limited the exploitation of this technique.

X-ray spectroscopy of the hot phase of the galactic superwind of the prototypical starburst galaxy M82 has revealed a much higher metallicity in the outer parts of the outflow and with a \(\alpha \)/Fe abundance ratio significantly higher than in the host galaxy, confirming that the outflow is associated with hot plasma freshly enriched by recent generation of core-collapse supernovae produced by the starburst event (Ranalli et al. 2008). These results provide the direct observational evidence that starburst superwinds eject metals with velocities of a few/several hundred km s\(^{-1}\), thus directly enriching the CGM and IGM.

Fig. 45
figure 45

Image reproduced with permission from Thomas et al. (2010), copyright by the authors

Representation of the specific star formation rate as a function of look-back time for galaxies of the different masses shown in the labels

X-ray spectroscopy of the hot plasma in galaxy clusters and galaxy groups has generally revealed Solar-like chemical abundances (Fig. 46) and a surprisingly flat radial distribution of the relative abundances (Simionescu et al. 2010, 2015, 2018; Mernier et al. 2017). One should not confuse these abundances with those of the stellar populations in the galaxies belonging to the clusters (which may be \(\alpha \)-enhanced if they formed rapidly); indeed even if star formation may have stopped in the cluster’s galaxies, SNIa keep enriching the ICM (with their typical high Fe/\(\alpha \) pattern) over time. Mernier et al. (2017) estimate that, on average, the fraction of SNIa with respect to the total number of SNe (i.e., SNIa + SNcc) that have contributed to the enrichment of the ICM must be in the range of 20–40%. Interestingly, de Plaa et al. (2017) find that the O/Fe abundance ratio does not depend on ICM temperature, therefore suggesting that the enrichment of the ICM is not related to cluster mass and that most of the enrichment has occurred before the ICM was formed.

Fig. 46
figure 46

Image reproduced with permission from Simionescu et al. (2018), copyright by the authors

Chemical abundances relative to iron inferred for the ICM of the Perseus cluster, obtained with XMM-Newton RGS for O, Ne and MG, and with the Hitomi SXS for Si to Ni

Within the context of X-ray spectroscopy, we mention that excellent high angular resolution maps of the metal enrichment of some individual clusters have been performed, revealing very interesting substructures. For instance, Sanders et al. (2016a) has obtained a detailed map of the metal enrichment of the Centaurus cluster revealing high metallicity blobs on scales of 5–10 kpc, which are likely tracing material uplifted by the AGN hosted in the central galaxy. Summarizing, the \(\alpha \)/Fe abundance ratio is an important clock to investigate the star formation history in galaxies. Massive galaxies are characterized by systematically higher \(\alpha \)/Fe ratio that, together with information on their stellar ages, indicates that more massive galaxies formed on shorter timescales and at earlier cosmic epochs than lower mass galaxies. The rapid star formation in massive systems is generally modelled in terms of a combination of enhanced star formation efficiency and strong feedback effect that rapidly quench star formation. The intracluster medium typically shows solar abundances, with no significant radial variation, reflecting the additional, continuous ejection of iron by SNIa over time also from passive galaxies.

7.2.2 N/O

In late-type galaxies, the analysis of chemical abundances has often focused on the gas phase and on those elements whose abundance can be inferred through the nebular emission lines (although DLA have also been extensively used to explore in detail the circumgalactic medium and outer parts of galactic discs, as discussed later on). Among these, nitrogen is one of the elements that has been subject to many extensive studies. Indeed, its abundance can be inferred from the relative bright doublet of [NII]\(\lambda \lambda \)6548,6584 next to H\(\alpha \). The nitrogen abundance can be inferred ‘directly’ through the \(T_\text {e}\) method, through auroral lines tracing the temperature of the partially ionized zone (e.g., Andrews and Martini 2013; Pilyugin et al. 2010a; Pérez-Montero and Contini 2009; Berg et al. 2011, 2013, 2015a), see Sect. 3.1. As noted earlier, the latter method requires the detection of very faint lines and hence can be applied only to limited samples of nearby galaxies/HII-regions or stacked spectra of galaxies. However, at least in star-forming galaxies, the N/O abundance ratio is nearly proportional to the [NII]\(\lambda \)6584/[OII]\(\lambda \)3727 line flux ratio, as these lines are emitted from nearly the same zone in HII-regions and hence their ratio is little dependent on other factors such as the ionization parameter and shape of the ionizing continuum. Therefore, since these are both relatively bright lines in most star-forming galaxies, the [NII]\(\lambda \)6584/[OII]\(\lambda \)3727 ratio can be used to investigate the N/O abundance in relatively large samples of galaxies, and calibrations of this diagnostic have been provided (e.g., Pérez-Montero and Contini 2009). In the absence of [OII]\(\lambda \)3727 (which requires a relatively large wavelength range to be observed together with [NII]\(\lambda \)6584), the [SII]\(\lambda \)6717,6730 doublet can be used as an alternative proxy of \(\alpha \) elements, so as to infer N/S from [NII]\(\lambda \)6584/[SII]\(\lambda \)6717,6730 (Pérez-Montero and Contini 2009).

Nitrogen is a particularly interesting element to investigate in galaxies, because, in contrast to oxygen and other \(\alpha \)-elements, it is produced primarily by intermediate-mass stars and only with a smaller contribution by massive stars [possibly enhanced in the presence of stellar rotation (Vangioni et al. 2018; Vincenzo and Kobayashi 2018, and references therein)]. Therefore, the N/O ratio provides precious information on the evolutionary stage of the galaxy. Moreover, nitrogen has also a “secondary” component, whose production increases with metallicity (Edmunds and Pagel 1978); indeed, being a product of the CNO cycle, its abundance increases at expenses of the C and O abundances. As a consequence, at high metallicities the nitrogen abundance is expected to evolve quadratically with the metallicity, \(\hbox {N/H} \propto (\hbox {O/H})^2\) or, equivalently, \(\hbox {N/O} \propto \hbox {O/H}\).

The nitrogen-to-oxygen abundance ratio in nearby galaxies has been investigated by multiple studies (e.g., Edmunds and Pagel 1978; Vila Costas and Edmunds 1992, 1993; Thuan et al. 1995; van Zee et al. 1998; Pérez-Montero and Contini 2009; Pilyugin et al. 2010a, 2012; Pérez-Montero et al. 2013; Andrews and Martini 2013; Berg et al. 2011, 2013, 2015a; Belfiore et al. 2015, 2017). Figure 47, top shows the nitrogen and oxygen abundances inferred from the direct-\(T_\text {e}\) method from SDSS galaxy spectra stacked in bins of SFR and mass; hollow points are in bins of stellar mass, while coloured points are further split in bins of SFR (Andrews and Martini 2013). At low metallicities the N/O abundance is relatively constant (if one ignores galaxies with high SFR, which will be discussed later); this is the region where nitrogen is thought to mostly have a ‘primary’ contribution from massive stars. At \(12 + {\text{{log}(O/H)}} > 8.3\)–8.4, N/O increases steeply with metallicity; this is the region where secondary nitrogen production is thought to take over and where intermediate mass stars start to contribute significantly.

Fig. 47
figure 47

Image reproduced with permission from Andrews and Martini (2013), copyright by AAS

N\(^+\)/O\(^+\) ratio as a function of direct-method oxygen abundance (a) and M\(_*\) (b). Open circles are obtained from SDSS stacks binned in stellar mass, coloured symbols result from stacks in mass and SFR

A similar trend has been clearly observed also for the nitrogen abundance of DLA, which populate mostly the low metallicity plateau (e.g., Pettini et al. 2002a; Centurión et al. 2003; Zafar et al. 2014). However, for DLA there is evidence for a bimodal distribution of this plateau, with most systems clustering at [N/\(\alpha \)]\(\sim \,-\,0.85\) and a smaller fraction of them (\(\sim \) 25%) clustering around [N/\(\alpha \)]\(\sim \,-\) 1.4. Such bimodal distribution may be related to the N/O spread observed in emission line galaxies (especially the strongly star-forming ones) at low metallicities (Fig. 47).

For local galaxies most studies have focused mostly on the use of nebular emission lines in star-forming galaxies. We shall warn that a non-negligible number of studies investigate the distribution of galaxies on the N/O versus O/H diagram using a nitrogen-based strong line diagnostic as a tracer of O/H, such as N2 or O3N2. This is a fundamental mistake that should be avoided by any means. Indeed, the use of the same information (the [NII] line flux) on both axes unavoidably introduces artificial correlations between N/O and O/H. Moreover, strong line diagnostics based on theoretical photoionization models (and several Bayesian methods) assume a-priori a relationship between N/O and O/H; therefore, these strong line diagnostics cannot be used to investigate the N/O vs O/H trends.

It has been pointed out by several authors that the N/O abundance ratio is a strong function of the galaxy mass (Pérez-Montero and Contini 2009; Andrews and Martini 2013; Masters et al. 2016). This is, for instance, shown in the stacked analysis by Andrews and Martini (2013) in Fig. 47. Such a trend with stellar mass is regarded as a consequence of the fact that massive galaxies are more evolved; hence the nitrogen enrichment contribution by intermediate mass stars has been more prominent. However, it is also probably a secondary product of the mass–metallicity relation; indeed the higher metallicity of massive galaxies likely boosts the production of secondary nitrogen.

Chemical evolutionary models have been proposed to interpret the evolution of the N/O abundance (e.g., Matteucci 1986; Garnett 1990; Coziol et al. 1999; Henry et al. 2000; Chiappini et al. 2005; Köppen and Hensler 2005; Torres-Papaqui et al. 2012). The most recent effort in this area is from Vincenzo et al. (2016) where the observational data are compared with the different model predictions by varying different parameters such as gas inflow properties, efficiency of star formation, outflow loading factor and also including the scenario of differential outflow rates, in which oxygen is expelled more preferentially by the SN-driven winds than nitrogen. Some of these models are shown in Fig. 48, overplotted on the density distribution observed in several thousand galaxies from the SDSS survey (colour shaded area). The broad distribution of galaxies in the N/O versus O/H diagram implies that different galaxies have evolved through different paths. However, models show that, on average, the global population of star-forming galaxies require an initial phase in which star formation in fuelled by gas infall over a timescale of 1 Gyr and that outflows start being effective in quenching further enrichment when the galaxy has reached a metallicity close to solar. These data seem to constrain the average past star formation efficiency to a value of about \(\nu \sim 1.5\)\(2\ {\mathrm{Gyr}}^{-1}\).

Fig. 48
figure 48

Image reproduced with permission from Vincenzo et al. (2016), copyright by the authors

Models of nitrogen enrichment in galaxies overplotted onto the SDSS data (colour shaded area). a Effect of varying the infall mass \(M_{\mathrm{inf}}\); b effect of varying the infall time-scale \(\tau _{\mathrm{inf}}\); c effect of varying the star formation efficiency \(\nu \); d effect of varying the outflow loading factor \(\omega \). In panel d, the dashed lines correspond to a non-differential outflow (where both N and O are expelled with the same efficiency) while the solid lines refer to the reference assumption of a differential outflow where N is not expelled (with \(\omega _{\mathrm{N}} = 0\))

It is interesting to note that the characteristic shape of the N/O versus O/H diagram for star-forming galaxies can be very useful to identify secondary evolutionary effects at play when galaxies deviate from this sequence. In particular, observations have shown that a significant fraction of galaxies and star-forming regions tend to scatter toward the high N/O region at fixed O/H (e.g., Köppen and Hensler 2005; Belfiore et al. 2015). One of such examples is illustrated in Fig. 49, which shows the N/O versus O/H diagram for the spatially resolved star-forming regions of a galaxy in the SDSS4-MaNGA sample (Belfiore et al. 2015). Such deviations can be explained through different scenarios: (1) a burst of star formation with increased star formation efficiency that, as shown in Fig. 48c, at later time boosts the N/O abundance relative to the sequence; (2) the infall of pristine/metal-poor gas at late times; such an event dilutes the overall metallicity leaving unaffected the N/O abundance, hence moving the galaxy/region horizontally on the diagram (Köppen and Hensler 2005; 3) a fountain scenario in which metal-rich gas with high N/O abundance is ejected from the central region and deposited on the outer galactic regions, which are more metal-poor and have lower N/O, resulting in a mixing sequence. Each of these scenarios is characterized by a different pattern on the N/O versus O/H diagram. For instance, in the case of the galaxy shown in Fig. 49 a fountain/mixing scenario seems to describe well the deviations from the N/O main sequence in this system. Another important mechanism that has been proposed to explain local enhanced nitrogen abundance in some star-forming regions (especially thanks to IFU techniques) is the enrichment of nitrogen by Wolf–Rayet outflows (Walsh and Roy 1989; López-Sánchez et al. 2007; James et al. 2009; Pérez-Montero et al. 2011; Monreal-Ibero et al. 2012).

Fig. 49
figure 49

Image reproduced with permission from Belfiore et al. (2015), copyright by the authors

Spatially resolved N/O vs O/H diagram for one galaxy in the SDSS4-MaNGA sample. Points are colour-coded according to their \(D_{\mathrm{n}}(4000)\) parameter (strength of the 4000 Å break), which is a tracer of the age of the stellar population (redder points are older). The contours show the distribution in the SDSS sample (with the same calibrators adopted for the galaxy). The blue and black solid lines with diamonds show the fountain mixing effect in which metal-rich gas expelled from the central region is mixed with metal-poor gas from the outer regions, for two different values of the lower metallicity. The deviation from the main N/O–O/H main sequence observed in this galaxy can be reproduced fairly well through such fountain mixing simple model. The correlation of the deviation with the \(D_{\mathrm{n}}(4000)\) parameter can be explained with the fact that older regions are more gas poor, and hence the dilution effect is more effective

The gradient of N/O abundance has also been investigated recently (Berg et al. 2013, 2015a; Belfiore et al. 2017; Esteban and García-Rojas 2018; James et al. 2009, 2013; Westmoquette et al. 2013; Kumari et al. 2018). Figure 50 shows the N/O abundance ratio radial gradient for star-forming galaxies from the SDSS4-MaNGA survey in bins of stellar mass (Belfiore et al. 2017). The systemically increasing nitrogen abundance with galaxy mass, the tendency for the gradient to steepen with galaxy mass and to flatten in the outer parts are all trends similar to those seen in the metallicity gradients (Fig. 31). However, one important difference is that the N/O gradient does not flatten in the central region as instead observed for the O/H abundance. Within the context of the inside-out growth of late-type galaxies, this finding indicates that the central regions of massive galaxies have locally evolved to an equilibrium metallicity (saturated around the yield), while the nitrogen abundance continues to increase as a consequence of both the delayed secondary nucleosynthetic production and the contribution from intermediate mass stars.

Fig. 50
figure 50

Image reproduced with permission from Belfiore et al. (2017), copyright by the authors

N/O abundance ratio radial gradient for star-forming galaxies from the SDSS4-MaNGA survey in bins of stellar mass

Summarizing, the N/O vs O/H diagram of galaxies shows a dual behaviour, with a plateau at low metallicities and a steeply increasing trend at high metallicities, which can be interpreted in terms of primary production of nitrogen at low metallicities, while at high metallicities can be interpreted as the delayed production by nitrogen by intermediate-mass stars together with the “secondary” production channel; hence the N/O vs O/H diagram can be interpreted as an evolutionary sequence. The nitrogen abundance also correlates significantly with galaxy stellar mass, which is interpreted in terms of more massive galaxies being more evolved; hence intermediate mass stars have had more time to enrich the ISM with nitrogen. Individual galaxies or galactic regions may show significant deviations from this trend (in particular by showing enhanced nitrogen abundance), which can be explained by models in terms of different effects (such as metallicity dilution by accreting near-pristine gas, variation of star formation efficiency, differential outflows effects, enhanced enrichment by Wolf–Rayet stars).

7.2.3 C/O

The abundance of carbon provides additional important information on the evolutionary stage of galaxies as carbon is primarily released by intermediate mass stars (although Wolf–Rayet stars, whose progenitors are thought to be high mass stars, are regarded as additional important contributors of carbon, Dray et al. 2003; Dray and Tout 2003). Therefore, significant carbon enrichment is generally delayed with respect to \(\alpha \) elements.

In the local universe the carbon abundance has been investigated in galactic stars through medium/high-resolution surveys (Gustafsson et al. 1999; Bensby and Feltzing 2006; Akerman et al. 2004; Spite et al. 2005; Fabbian et al. 2009; Nieva and Przybilla 2012; Nissen et al. 2014; Tautvaišiene et al. 2016), although one should be aware that when observed in evolved giant and supergiant stars complex internal mixing and dredge up processes make it difficult to properly interpret the observed abundances and to use them to trace galaxy evolution. Carbon abundance in stellar spectroscopy has been extended also to some nearby galaxies, and recently Conroy et al. (2014) have extended the analysis to thousands of galaxies from the SDSS, through fitting of stacked spectra (see Fig. 44).

The gas phase metallicity has been more difficult to determine as carbon does not have strong transitions at optical wavelengths. However, deep observations of bright HII regions, both in our galaxy and external galaxies, have enabled the detection of CII recombination lines in the optical, enabling the measurement of the carbon abundance for some of them (Esteban et al. 2002, 2009, 2014; Peimbert et al. 2005; López-Sánchez et al. 2007; García-Rojas and Esteban 2007). HST has enabled sensitive spectroscopy in the UV, where collisionally excited lines of carbon are present, especially CIII]\(\lambda \lambda 1907,1909\) (but also CIV1459, for galaxies with harder ionizing spectrum, such as AGNs), which have been used to estimate the carbon abundance in HII regions, although often requiring photoionization modelling (Garnett et al. 1997, 1999; Kobulnicky et al. 1997; Kobulnicky and Skillman 1998; Izotov and Thuan 1999; Berg et al. 2016a; Pérez-Montero and Amorín 2017; Peña-Guerrero et al. 2017). As mentioned in Sect. 3.5.2 if information on the gas temperature is available, then the carbon abundance can be inferred with higher accuracy (Garnett et al. 1995, 1999). The UV spectral range also contains absorption features from UV resonant lines which can further be used, through HST data, to constrain the carbon abundance of the ISM (primarily using the CII1334 ISM absorption, James et al. 2014a) and of the young stellar population (Leitherer et al. 2011; Leitherer 2011).

Figure 51, from Berg et al. (2016a), summarizes some of the main findings on the C/O versus O/H diagram for different galactic systems and components, specifically halo and disk MW stars, HII regions measured either through recombination lines or through UV collisionally excited lines (as well as high-z DLA, which will be discussed in the next section). The plot shows large dispersion. However, at least at metallicities higher than \(12+\log ({\mathrm{O/H}})=7\), it resembles the trend observed in the N/O versus O/H diagram, with a flat relation at low metallicity and a steeply increasing abundance at high metallicities. The latter trend has led some authors to suggest that carbon may also have a secondary production, i.e., yields that are strongly metallicity-dependent; alternatively, or in addition, the delayed release of C may also mimic a secondary production effect (Garnett et al. 1999; Henry et al. 2000; Carigi 2000; Chiappini et al. 2003). The similar behaviour of carbon and nitrogen is further supported by the observed C/N trend (Fig. 52), which is constant with metallicity (Berg et al. 2016a) and which has strengthened the idea that carbon follows an enrichment pattern similar to nitrogen.

Fig. 51
figure 51

Image reproduced with permission from Berg et al. (2016a), copyright by AAS

C/O abundance ratio as a function of O/H for Galactic and extragalactic systems, compared with various models. Blue dots are HII regions in dwarf galaxies observed by Berg et al. (2016a) with HST/COS. Purple dots are other objects with direct oxygen abundances and C/O abundances determined from UV CELs. Green filled squares are star-forming galaxies with abundances based on RLs. Triangles are MW halo stars, while 4-pointed stars are disk stars. Finally, orange diamonds are DLAs, and lines are the results of three enrichment models, see Berg et al. (2016a) for details

Fig. 52
figure 52

Image reproduced with permission from Berg et al. (2016a), copyright by AAS

C/N abundance ratio as a function of O/H for the same systems as in Fig. 51. The dashed line is the weighted average of the significant detections based on CEL

The stacking analysis of Conroy et al. (2014) has revealed that carbon is enhanced relative to iron in more massive galaxies (i.e., with larger velocity dispersion, Fig. 44), an effect similar to the \(\alpha \) enhancement in massive galaxies, although slightly less extreme, indicating that carbon is capturing star formation on intermediate temporal scales.

Summarizing, the C/O vs O/H diagram of galaxies shows a similar dual behaviour as for nitrogen, with a plateau at low metallicities and a increasing trend at high metallicities. These similarities with the N/O diagram suggests that carbon shares a common origin with nitrogen (as confirmed by the constant C/N ratio). In particular, at high metallicities carbon is enriched with delay by intermediate mass stars and may also have a “secondary” component. Therefore, the C/O vs O/H diagram describes a temporal sequence. In the next section we will discuss the peculiarity of the C/O abundance ratio at very low metallicities.

7.3 High redshift and the very low metallicity regime

At high redshift the measurement of the relative chemical abundances are obviously much more difficult to obtain as the limited S/N of distant galaxies often prevents detecting the required multiple spectral diagnostic features. This is especially true for what concerns measuring the [\(\alpha \)/Fe] ratio in the spectra of stellar populations of high-z galaxies, as this requires excellent S/N on the continuum. Indeed, currently, the [\(\alpha \)/Fe] ratio has been measured only in a few cases by exploiting individual spectra (Lonoce et al. 2015; Kriek et al. 2016) or by using stacking of star-forming or quiescent galaxies (Onodera et al. 2015; Steidel et al. 2016). The general result is that these distant galaxies are all \(\alpha \)-enhanced, in some cases even relative to the lower redshift passive galaxies of the same mass (e.g., Fig. 53), indicating that these galaxies have formed quickly, with little iron pollution from SNIa, generally inferring star formation timescales of only 0.5–1 Gyr.

Fig. 53
figure 53

Image reproduced with permission from Kriek et al. (2016), copyright by Macmillan

[Mg/Fe] abundance for a massive galaxy at \(z=2.1\) as a function of mass (left) and as a function of [Fe/] (right), compared with lower-redshift galaxies (see legend). The black dashed line and the red arrows are the results of two evolutionary models.

The N/O abundance ratio in high-z star-forming galaxies has been investigated little so far, due to the need to sample a broad range of nebular diagnostics (hence multi-band near-IR/optical observations are needed, which may be subject to differential slit losses if not performed with IFU), both to measure N/O (which requires measuring [NII] and [OII]) and to measure the metallicity O/H with diagnostics that do not use nitrogen. One of such attempts is shown in Fig. 54 based on a near-IR spectroscopic survey of galaxies at \(z\sim 2.3\) (Steidel et al. 2016; Strom et al. 2017a) Generally high-z star-forming galaxies (symbols) follow the same relation as local galaxies (contours), although with somewhat larger scatter, some tendency of being more nitrogen-rich at a given O/H. This seems confirmed (especially in terms of larger scatter) by a smaller sample, but based on \(T_\text {e}\) measurements, at lower metallicities, in the work by Kojima et al. (2017). The larger scatter, if confirmed with higher statistics in future surveys, may be a consequence of enhanced star formation efficiency at such early epochs (Vincenzo and Kobayashi 2018), or may reflect frequent and prominent inflows of near-pristine gas in high-z galaxies (which dilute the metallicity, hence O/H, but affect little N/O), as expected by many theoretical models (see the discussion in Sect. 7.2.2).

Fig. 54
figure 54

Image reproduced with permission from Strom et al. (2017a), copyright by the authors

N/O versus O/H for a sample of star-forming galaxies at \(z\sim 2.3\) (symbols) compared with local star-forming galaxies (contours)

The C/O ratio was studied by Steidel et al. (2016) in the stacked spectrum of \(\sim 20\) galaxies at \(z\sim 2.4\) and by Amorín et al. (2017) in individual galaxies at \(z\sim 3\). They found values in agreement with the values in local stars and HII regions, although with large scatter.

As discussed in Sect. 3.6, high column density absorption systems (DLA) generally provide the most accurate determination of the relative (and often absolute) chemical abundances (e.g., Prochaska and Wolfe 1999; Berg et al. 2015b, 2016b; Fumagalli 2014; Pettini et al. 2008; Henry and Prochaska 2007; Dessauges-Zavadsky et al. 2006; Wolfe et al. 2005; Prochaska et al. 2003; Rafelski et al. 2012; Neeleman et al. 2013; Jorgenson et al. 2013), although they probe a broad range of environments, whose connection with galaxies is generally not fully clear (likely ranging from outskirts of galactic discs to clumps in the intergalactic medium). In contrast to early claims that DLA abundances may resemble the chemical pattern of stars in the MW halo, hence that DLA may probe the formation of galactic haloes, more recent studies have shown that the chemical and kinematic properties of DLA are more similar of those seen in dwarf galaxies of the Local Group (Fig. 55; Cooke et al. 2015; Berg et al. 2015b; De Cia et al. 2016), hence DLA might be tracing the early formation of dwarf satellite galaxies (also based on their similar velocity dispersion).

Fig. 55
figure 55

Image reproduced with permission from Cooke et al. (2015), copyright by AAS

[\(\alpha \)/Fe] versus [Fe/H] for high redshift DLA (black symbols) compared with local dwarf spheroids

In many DLA it is possible to trace the detailed chemical enrichment pattern of several elements (Fig. 56), in most cases confirming that these result from the enrichment of core-collapse supernovae from massive stars (Prochaska et al. 2003; Dessauges-Zavadsky et al. 2006).

Fig. 56
figure 56

Image reproduced with permission from Dessauges-Zavadsky et al. (2006), copyright by ESO

Detailed chemical abundances inferred for a sample of 11 DLAs at \(z\sim 2\)

It is interesting the finding that, for both DLAs at high-z and for halo stars, the carbon abundance relative to \(\alpha \) elements increases systematically at very low metallicities (\(12+\log ({\mathrm{O/H}})<7\), e.g., Figs. 51 and 57 Cooke et al. 2012; Berg et al. 2016b; Lehner et al. 2016) which has been interpreted as possible signature of enrichment by PopIII stars (Carigi and Peimbert 2011), although such high values of C/O at low metallicities can also been explained by other models without invoking the contribution of PopIII-like yields, but simply different carbon yields in the metal poor (PopII) regime (Mollá et al. 2015). A similar enrichment has been observed in a population of metal-poor halo stars collectively known as carbon-enhanced metal-poor [CEMP] stars (Beers and Christlieb 2005).

Lyman Limit Systems (LLS), which are characterized by lower absorbing column density relative to DLA, appear to deviate from the trend observed in MW stars and DLA, by being carbon enhanced at metallicities \(-2< [\alpha /{\mathrm{H}}] < -1\). Lehner et al. (2016) suggest that this indicates that LLS trace gas clouds enriched by preferential ejection of carbon from low metallicity galaxies.

Even more interesting is the discovery of extremely low metallicity DLA and LLS systems (\([\alpha /{\mathrm{H}}]<-3\)) that are also carbon poor (Fig. 57) (Lehner et al. 2016; Crighton et al. 2016; Cooke et al. 2017). These are considered the best candidates as tracers of gas polluted by PopIII stars. Overall, their chemical enrichment pattern is well-reproduced by pollution from supernovae originating from PopIII stars, with progenitor masses of about \(20~\hbox {M}_{\odot }\). The abundance pattern of these systems is partly reminiscent of some of the most metal poor halo stars. The most metal poor of these stars has a metallicity of \(4.5 \times 10^{-5}~\hbox {M}_{\odot }\) (Caffau et al. 2011) and its chemical abundance pattern is well reproduced by a metal-free progenitor with mass of about 20–30 \(\hbox {M}_{\odot }\) (Fig. 58) (Schneider et al. 2012b). Interestingly, both in the case of DLA/LLS and halo stars the abundance pattern excludes the scenario of pollution by very massive PopIII progenitors, which would result in hypernovae. It is also very interesting that the extremely metal poor halo stars have sub-solar stellar masses, at metallicities well below the critical value that, according to standard model, would allow cooling and fragmentation of the gas that would enable the formation of low-mass stars (shaded region in Fig. 57; Frebel et al. 2007). While initially very puzzling, Schneider et al. (2012a, b) pointed out that small amount of dust formed in the ejecta of PopIII SNe would be enough to enable the cooling and fragmentation of the gas that would result into the formation of the first generation of extremely metal-poor low-mass stars.

Fig. 57
figure 57

Image reproduced with permission from Lehner et al. (2016), copyright by AAS

[C/\(\alpha \)] vs [\(\alpha \)/H] for high-redshift DLAs and Lyman Limit Systems, compared with Galactic stars. The hatched orange region identifies the area within which gas may have been primarily polluted by PopIII stars (Frebel et al. 2007)

Fig. 58
figure 58

Image reproduced with permission from Schneider et al. (2012b), copyright by the authors

Abundance pattern observed in the most metal poor halo star (black points ,Caffau et al. 2011) compared with yields expected from PopIII stars (red points).

We conclude this section by highlighting that millimetre/submillimetre observations of molecular transitions at high-redshift are now sensitive enough to provide valuable constraints not only on the composition of molecular species in the ISM of primeval galaxies, but also on the relative abundance of different atomic isotopes, which can provide precious information on the properties of the stellar populations responsible for the early chemical enrichment. Zhang et al. (2018b) have measured the \(^{13}\hbox {C}/^{18}\hbox {O}\) abundance ratio in a sample of distant lensed, starburst galaxies (\(z\sim 2\)–3), by measuring multiple transitions of the \(^{13}\hbox {CO}\) and \(\hbox {C}^{18}\hbox {O}\) isotopologues of carbon monoxide. Since \(^{13}\hbox {C}\) is mostly produced by low/intermediate mass stars (\({M}_{*}<8\,\hbox {M}_{\odot }\)) while \(^{18}\hbox {O}\) is mostly produced by massive stars (\({M}_*>8\,\hbox {M}_{\odot }\)), the \(^{13}\hbox {C}/^{18}\hbox {O}\) ratio is sensitive to the shape of the IMF. From their measurements. Zhang et al. (2018b) find that these early systems are likely characterized by a top-heavy IMF.

Summarizing, galaxies with \(\alpha \)-enhanced stellar populations are already seen at high redshift (\({z}\sim 1-2\)), indicating that these systems have formed rapidly, likely within 0.5–1 Gyr. The N/O abundance ratio of distant galaxies generally follow the local relation, but with larger dispersion and with a larger fraction of nitrogen-enhanced (relative to oxygen) galaxies, which may indicate that these galaxies have experienced enhanced star formation efficiency or absolute metallicity dilution by infalling near-pristine gas. The \(\alpha \)/Fe properties of DLAs (together with their velocity dispersion) suggest that they may trace the early formation of dwarf galaxies and generally be primarily enriched by massive stars (also based more broadly by a wider set of chemical abundances). The properties of very metal poor DLA with sub-solar carbon abundances suggest that these may trace the early enrichment by the first generation of stars (PopIII).

8 Metallicity and chemical abundances in AGN

Active galactic nuclei (AGN), powered by supermassive accreting black holes, in their various manifestations, have been extensively investigated to probe the metallicity in galactic nuclei, in their host galaxies and even in the CGM. Given that they can reach very high luminosities (quasar phase) they have been effectively used to probe the metallicity of their circumnuclear region and of their host galaxies out to very high redshift.

The nebular emission lines in AGN are primarily divided in two classes, “broad” lines and “narrow” lines.

Broad emission lines, with line widths up of a few to several thousands km s\(^{-1}\), are emitted by a nuclear region typically smaller than a fraction of a parsec (the so-called Broad Line Region, BLR). The density of the clouds in the BLR is so high (\(\sim 10^{11}~{\mathrm{cm}}^{-3}\)) that “forbidden” transitions (e.g., [OIII]5007, [OII]3727, [SII]6730, etc.) are not detected; indeed the critical densities of all these collisionally excited transitions are well below the typical density of the BLR, implying that in this regime their emissivity increases only linearly with gas density, in contrast to the permitted lines (such as hydrogen Balmer recombination lines) whose emissivity increases quadratically with density even in the extreme conditions of the BLR.

The narrow lines have widths more comparable to those typically observed in the host galaxies (a few 100 km s\(^{-1}\)), although typically broader because often associated with outflows, and extend on scales ranging from a few 100 pc to several kpc (the so-called Narrow Line Region, NLR).

Metallicity determinations of the BLR and NLR have mostly relied on photoionization models (e.g., Hamann and Ferland 1999; Nagao et al. 2006c), although attempts have been made to use the direct-\({T}_{\mathrm{e}}\) method (Dors et al. 2015) but which have revealed the inadequacy of this method for AGNs ( the origin of such “temperature problem” in AGNs is not yet clear).

Some of the broad line ratios from metal transitions in the UV have been proposed as sensitive metallicity tracer of the BLR. The nebular emission ratio \((\hbox {SiIV}\lambda 1397+\hbox {OIV}\lambda 1402)/\hbox {CIV}\lambda 1549)\) has been proposed to be the most stable against distribution of gas densities and ionization parameter in the BLR clouds, and also in terms of hardness of the ionizing continuum (Nagao et al. 2006a). The ratios \(\hbox {NV}\lambda 1240/\hbox {CIV}\lambda 1549\) and \(\hbox {NV}\lambda 1240/\hbox {HeII}\lambda 1640\) have also been proposed (e.g., Hamann and Ferland 1999; Dietrich et al. 2003b; Nagao et al. 2006c; Matsuoka et al. 2011b; Wang et al. 2012), but they are more sensitive to ionization parameter, shape of the ionizing continuum and are primarily sensitive to the nitrogen abundance rather than metallicity.

The general finding is that the metallicity of the BLR in quasars is very high, nearly always supersolar and up to several times solar (Hamann and Ferland 1999; Dietrich et al. 2003b; Nagao et al. 2006a; Jiang et al. 2007; Juarez et al. 2009; Simon and Hamann 2010; Matsuoka et al. 2011b; Wang et al. 2012; Shin et al. 2013; Xu et al. 2018). Such high values of the metallicity have posed questions on whether the photoionization modelling of the extreme environment characterizing the BLR is appropriate. However, very high nuclear metallicities (a few/several times solar) are also confirmed by the iron emission and absorption features observed in the X-ray emission coming from the nuclear region (e.g., Jiang et al. 2018). Yet, such high metallicities in the nuclear region of AGN are not really unexpected. Indeed the very high densities and large amount of gas in the central region of AGN likely foster rapid star formation and quick enrichment on the ISM. Moreover, it is important to bear in mind that the mass of gas in the BLR is very small, a few times \(10^4~\hbox {M}_{\odot }\). As pointed out by Juarez et al. (2009), such small mass can be quickly enriched with less than a SN explosion every \(10^4\) yrs.

The somewhat puzzling result is that the metallicity of the BLR does not seem to evolve with redshift (Fig. 59, top) and such lack of evolution appears to persist out to the most distant quasars known, at \(z=7.5\) (Fig. 59, bottom) (Nagao et al. 2006c; Juarez et al. 2009; Mortlock et al. 2011; Bañados et al. 2018; Xu et al. 2018).

Fig. 59
figure 59

Image reproduced with permission from Bañados et al. (2018), copyright by Macmillan

Top: average metallicity of the BLR in quasars (from stacked spectra) as a function of redshift. Image reproduced with permission from Nagao et al. (2006c), copyright by ESO. Bottom: spectrum of the most distant quasar currently known (\(z=7.5\)) compared with the average spectrum of intermediate redshift quasars from SDSS, illustrating that the two are nearly identical, suggesting similar chemical enrichments of the BLR

The lack of redshift evolution of the metallicity which, in contrast to what is observed for galaxies, seems to remain high at all redshifts, is likely a consequence of the mass–metallicity relation combined with selection effects. Indeed, in order to be selected in large-scale surveys, the quasar has to be luminous enough to pass the sensitivity threshold of the survey; this generally implies that (even if accreting at the Eddington limit) the black hole must have already become fairly massive. If some form of black hole/galaxy relation is already in place at high redshift, this implies that, at any redshift, the host galaxy must already be massive, hence typically display high metallicity when the quasar enters into the survey (at any epoch). This combination of effects, explaining the lack of redshift evolution of the BLR metallicities, was discussed more quantitatively in Juarez et al. (2009).

The fact that the metallicity in AGN is linked to the mass–metallicity relation of the host galaxy was first hinted by the fact that the metallicity of the BLR scales with AGN luminosity (Hamann and Ferland 1999; Nagao et al. 2006c; Xu et al. 2018). Indeed, if the AGN luminosity is a function of the black hole mass (assuming an average \(L/L_\mathrm{Edd}\) ratio) and the black hole mass is linked to the host galaxy mass through BH-spheroid relation, then one would expect the AGN nuclear metallicity to scale with the AGN luminosity as a consequence of the mass–metallicity relation of the host galaxy (Juarez et al. 2009). A more clear evidence of this, which by-passes the use of AGN luminosity, is the more direct relationship between BLR metallicity and black hole mass obtained by Matsuoka et al. (2011b) and Xu et al. (2018), as illustrated in Fig. 60 (see also Ludwig et al. 2012, for a potential extension of the relation to low masses.)

Fig. 60
figure 60

Image reproduced with permission from Xu et al. (2018), copyright by the authors

Metallicity of the BLR as a function of black hole mass for quasars divided into bins of redshift (see legend)

Matsuoka et al. (2011b) report the lack of correlation between metallicity and Eddington ratio \(L/L_\mathrm{Edd}\). A correlation with the Eddington ratio is seen only when using the NV/CIV and NV/HeII ratios, which are primarily tracing the nitrogen enrichment. Since nitrogen enrichment is delayed with respect to \(\alpha \) elements, the latter correlation has been interpreted as indication that black hole accretion is triggered with a delay of a few 100 Myr with respect to the onset of star formation. This is delay that has been suggested also in local AGNs (Davies et al. 2007) and it has been interpreted as a consequence of the initial strong turbulence induced by SNe, which may prevent effective accretion onto the BH, while at later epochs the more gentle stellar winds may be effective in removing angular momentum from the gas (hence enabling it to move towards the centre) without introducing excessive turbulence or gas removal through SN-driven winds.

Many authors have used the flux ratio of the MgII 2798Å doublet relative to the UV FeII “bump” (due to a blending of multiplets at 2200–3090 Å), with the goal of constraining the redshift evolution of the \(\alpha \)/Fe ratio in the Broad Liner Region (e.g., Dietrich et al. 2003a; Maiolino et al. 2003; Freudling et al. 2003; Iwamuro et al. 2004; Jiang et al. 2007; De Rosa et al. 2011, 2014; Calderone et al. 2017; Mazzucchelli et al. 2017). These various studies find no evidence for a redshift evolution of the MgII/FeII flux ratio out to most distant quasars at \(z\sim 6.5\). If one assumes that the flux ratio is, to a first order, a proxy of the \(\alpha \)/Fe abundance ratio, then the lack of evolution would imply that the relative contribution of SNIa and core-collapse does not change with redshift. However, ones has to take into account that both emission features, and especially the FeII “bump”, are primary coolants of the BLR and, therefore, their flux does not really scale linearly with abundance, but it rather tends to adjust in order to keep the thermal equilibrium of the BLR clouds (Verner and Peterson 2004; Verner et al. 2004), therefore the lack of evolution of the MgII/FeII ratio may simply reflect the “thermostatic” role of the associated transitions.

Similar studies have been performed to investigate the metallicity of the NLR, i.e., on much larger galactic scales in AGN hosts. In this case studies have generally focused on type 2 AGNs, in which the BLR (whose strong, broad lines would otherwise prevent a proper disentangling of the flux of the narrow lines) is obscured along the line of sight. Studies have both exploited optical narrow nebular lines ratios (especially in local galaxies) and UV nebular lines (especially in distant galaxies, whose UV lines are redshifted into the optical bands) and by using photoionization models to infer the metallicity from a combination of line ratios. Initial claims about highly supersolar metallicities in the NLR (Groves et al. 2006) have been revised downwards; however, it is still true that most NLR appear to be metal rich, with metallicities around solar or super-solar (Nagao et al. 2006b; Matsuoka et al. 2009, 2011a; Stern and Laor 2013; Coil et al. 2015; Dors et al. 2015, 2017; Castro et al. 2017). Matsuoka et al. (2009) suggest that in the redshift range 1–4 there is little cosmic evolution of the NLR metallicity, although admittedly the statistics are much poorer than for the BLR metallicities; however, Coil et al. (2015) have pointed out that, in their sample of type 2 AGNs at \(z\sim 2.3\), the NLR metallicity (inferred from rest-frame optical diagnostics) is lower than the metallicity in the NLR of local AGN. They also point out that the metallicity of the NLR in their sample of type 2 AGN is higher than in a matched sample of star-forming galaxies.

The solar/super-solar metallicities in the NLR of AGNs, already at high redshift, can be partly explained in terms of dust destruction in the NLR, which releases metals into the ISM (Nagao et al. 2006b; Matsuoka et al. 2009; Dors et al. 2014), but probably is also partly due to the fact that the NLR is often also associated with galactic (AGN-driven) outflows originating from the central, metal-rich region of the galaxy. Indeed, several studies have also directly investigated the outflowing gas in quasars and AGNs, especially through absorption lines (and especially in Broad Absorption Line Quasars, where prominent blueshifted absorption troughs probe powerful winds), revealing high metallicity gas being expelled on kpc scales (Hamann and Ferland 1999; Simon and Hamann 2010; Ganguly et al. 2003, 2006; D’Odorico et al. 2004; Gabel et al. 2006; Arav et al. 2007; Borguet et al. 2012; Shin et al. 2017).

However, it has been suggested that the NLR also follows the (host galaxy) mass–metallicity relation (although offset towards higher values), either directly (Matsuoka et al. 2018) or indirectly through the AGN luminosity or black hole mass as indirect tracers of the host galaxy mass (Matsuoka et al. 2009; Ludwig et al. 2012).

We finally mention that, while the NLR tends to be generally metal rich, the investigation of the outer, most extended region of the NLR (sometimes referred to as Extended Narrow Line Region, ELR) does reveal low metallicity regions (e.g., Fu and Stockton 2009; Husemann et al. 2011) indicating that the outer parts of the NLR probe gas in the outer galaxy that are still poorly enriched.

Of growing interest is becoming the technique of probing the metallicity of the halo of quasars through the analysis of associated absorption systems detected in the spectrum of a nearby (in projection) background quasar (Prochaska et al. 2013, 2014). This technique has revealed that the circumgalactic medium of quasars at \(z\sim 2\) hosts significant quantities of cold gas (\(10^{10}~\hbox {M}_{\odot }\)) significantly metal enriched (\({Z}>0.1~{\hbox {Z}}_{\odot }\) out to the virial radius (\(r_{\mathrm{vir}}\sim 160\) kpc), implying that by \(z\sim 2\) feedback has already been quite effective in enriching the CGM of massive galaxies and also implying that the assumption of pristine gas accretion in many models may be inappropriate.

In summary, the BLR show very high metallicities at any redshift, with no indication of evolution with time. Selection effects are probably contributing to hide signs of redshift evolution; nevertheless, this points out toward very early enrichments of the central regions of massive galaxies. The metallicity of the BLR correlates strongly with the black hole mass, which is likely a result of the mass–metallicity relation of the host galaxy combined with the black hole–host galaxy mass scaling relation. The NLR also show high metallicities, though much lower than those observed in the BLR. The NLR metallicity also shows no evidence for redshift evolution. The interpretation is complex because contributions to the metallicity of the NLR are expected to come from the AGN-driven outflows, from the effects of dust destruction, and also to be linked to the mass–metallicity relation of the host galaxy. The investigation of absorption systems in the proximity of quasars (by exploiting pairs of quasars that are close in projection) has revealed large amounts of metal-enriched cold gas in their halos, suggesting that quasar activity has polluted significantly their circumgalactic medium through outflows, already at early cosmic times.

9 Metal budget

Since metals are produced by the star formation activity across the cosmic epochs, comparing the total amount of metals seen in the various phases to the cosmic evolution of stellar mass and SFR is useful to investigating consistency of the interpretation of these various independent observational results, and in particular the validity of the underlying assumptions and models about the production of chemical elements and their transfer among the various galactic and intergalactic phases. In other words, the metal budget is fundamental information whose evolution should be in agreement with the other independent observations of galaxy evolution and should be matched by models. For this reason it has been subject of considerable work (Pei and Fall 1995; Edmunds and Phillipps 1997; Pettini 2004, 2006; Ferrara et al. 2005; Bouché et al. 2007; Gallazzi et al. 2008; Zahid et al. 2012; Peeples et al. 2014; Madau and Dickinson 2014).

Determining the total mass in metals is a challenging goal, as heavy elements are produced in galaxies and then dispersed in the universe in different forms, often in states that, as we have discussed in this review, are difficult to observe. Moreover, in general observational studies measure the metallicity of the various phases, i.e., the abundance of metals relative to the total content of baryons. Therefore, even in those cases where the metallicity is well constrained, inferring the absolute total content of metals implies having a good knowledge of the total mass associated with the same phase (stars, ISM, CGM, ICM, WHIM, IGM), which are all distributed on different scales and forms, hence adding to the problem a additional level of complexity and uncertainty.

Attempts to infer the total mass in metals in the local universe have implied combining the contribution to the metal budget from all these components. More specifically:

  • Stars and ISM. Galaxies contain large amounts of metals locked into stars and star remnants and dispersed into the ISM. The mass of metals contained in galaxies can be obtained by integrating the mass function of galaxies convolved by the metallicity and gas fraction as a function of mass. As the chemical abundance ratios and the overall metallicities depend on galaxy types, it is necessary to make the computation dividing galaxies in bins of morphology or SF histories. Roughly, mass-weighted solar metallicities are obtained for this component when assuming a Salpeter IMF (Calura and Matteucci 2004; Gallazzi et al. 2008).

  • CGM and IGM. The CGM is enriched by galactic winds and plays a critical role in the current “equilibrium models” (see Sect. 4) because it constitutes a reservoir of metals extracted from the galaxy by winds and that could rain back onto the galaxy. Similarly, the IGM is thought to be enriched by galactic winds, especially those escaping from low-mass galaxies. As discussed in the previous sections, the chemical abundance of the CGM and IGM is obtained from the absorption features of several ionization species of various elements and shows significant evolution with redshift, rising to about 0.1 solar in the local universe and containing about 10% of the metals produced (Meiksin 2009; D’Odorico et al. 2010, 2016; Simcoe et al. 2011; Shull et al. 2014). While the amount of metals in galaxies declines steadily with redshift, the amount of metals in the CGM/ISM as traced by DLAs seems to remain roughly constant out to \(z\sim 4\) (Prochaska et al. 2013; Rafelski et al. 2014), implying that the fraction of metals in the CGM/IGM is much larger at high redshift than locally. It had been claimed that the amount of metals in DLAs shows a decline at \(z>4\) (Rafelski et al. 2014); however, this has not been confirmed by more recent studies (De Cia et al. 2018; Poudel et al. 2018).

  • The intracluster medium (ICM). As already mentioned, X-ray observations have revealed that the intracluster medium is highly enriched with all metals (generally with solar-like relative abundances, Sect. 7.2.1). Metals in clusters are mainly produced in the evolved populations of early-type galaxies (Matteucci 2012), which are enriched in Fe by ongoing production of type Ia SNe, also long after the end of star formation (Maoz and Mannucci 2012). The metal transfer toward the ICM is due either to AGNs, to SN explosion or to ram stripping, and the system as a whole evolves nearly as a closed box. As already mentioned, extensive works have been undertaken to estimate the content of metals in the ICM, especially with the advent of high-resolution spectroscopy (Mushotzky and Loewenstein 1997; Balestra et al. 2007; Blanc and Greggio 2008; de Plaa 2013; Molendi et al. 2016; Mernier et al. 2016, 2018; Hitomi Collaboration 2017; Simionescu et al. 2018). The metallicity of the ICM is generally very high (supersolar), and therefore, despite the small contribution to the total baryonic mass (\(\sim 4\%\)), it is expected to contribute significantly to the total metal budget (see discussion later in this section).

  • The warm–hot intergalactic medium (WHIM). The WHIM is the warmer/hotter phase of the IGM, thought to result from the gravitational shock-heating of the intergalactic medium in the local universe to temperatures of \(\sim 10^5\)\(10^7\) K, where most of the baryons in the local universe are though to reside (Nicastro et al. 2017). It is thought to be the local (hotter) counterpart of the (cooler) IGM at high-z observed through the Ly\(\alpha \) Forest. Its content of metals has been inferred through UV and X-ray absorption spectroscopy (e.g., Tripp and Savage 2000; Fang and Bryan 2001; Prochaska et al. 2004; Cooksey et al. 2008; Fang et al. 2010; Zappacosta et al. 2010, 2012). Typically the inferred metallicity is of the order of 0.1 \({\hbox {Z}}_{\odot }\) and its contribution to the local, total metal budget is less than 5%.

Amid these issues and uncertainties, the published values for the total amount of metals show a significant scatter, with values around \(\varOmega _Z\sim 10^{-4},\) (Pei and Fall 1995; Madau and Shull 1996; Zepf and Silk 1996; Mushotzky and Loewenstein 1997; Madau et al. 1998; Pagel 2002; Dunne et al. 2003; Calura and Matteucci 2004; Gallazzi et al. 2008) where \(\varOmega _Z\) is the density of metals normalized to the critical density of the Universe for \(h=0.7\), \(\rho _c=1.36\times 10^{11}\,{\hbox {M}}_\odot \) Mpc\(^{-3}\).

The relative distribution of metals among these different components is clearly still subject to significant uncertainties.

Expectations from the integrated cosmic production of stars can be achieved by integrating the star formation rate density as a function of redshift and convolving it with the yields per stellar generation (Mollá et al. 2015; Vincenzo et al. 2016). This is obviously subject to additional uncertainties, not only associated with the measurement of the evolution of the star formation history (Madau and Dickinson 2014), but also with our knowledge on the return fraction of metals to the gas phase and with our yet limited knowledge of the IMF and its potential variations. Bearing in mind all these uncertainties, the expected average metallicity of the local universe is inferred to be \(Z\sim 0.09\) solar for a Salpeter IMF and to decrease by one order of magnitude by \(z=2.5\) (Madau et al. 1998; Pettini 2006; Madau and Dickinson 2014).

There is reasonable agreement between the expected and measured metal budget in the local universe. However, given large uncertainties in both the measured and expected content of metals the agreement is not too surprising and really not very constraining of any of the underlying assumptions.

It is has been perhaps more instructive to investigate the metal budget in individual systems, as this exercise may provide information on processes resulting in loss of metals, or even provide constraints on the yield of metals.

For instance, Renzini and Andreon (2014) compare the amount of iron in the ICM (as inferred from X-ray observations) with the amount of iron expected to be produced by the stellar populations of galaxies within the cluster, based on empirical yields of iron. While they find a good agreement for intermediate mass clusters (\({M}_{500}\approx 10^{14}~{\hbox {M}}_{\odot }\)), in more massive clusters they reveal a clear tension, in the sense that the ICM contains much more iron mass (up to a factor of \(\sim 6\)) than that produced by stars in galaxies, revealing higher rates of type Ia SNe in clusters (e.g. Mannucci et al. 2008; Friedmann and Maoz 2018, and references therein), or issues either in the metallicity measurements or with our knowledge of the yields.

Calura and Matteucci (2004), Bouché et al. (2007), Peeples et al. (2014) and Tumlinson et al. (2017) present an extensive analysis of the metal budget in galaxies and including their CGM, by exploiting the extensive results from the COS-Halos project (Tumlinson et al. 2011; Prochaska et al. 2017). They show that, nearly independently of mass, only about 20–25% of metals produced in stars remain in galaxies (in stars or in the ISM). They infer that, for \({L}_*\) galaxies, as much as 40% of the metals produced by stars are deposited in the halo (CGM, within a radius of \(\sim \) 150 kpc), while the remaining must be lost into the intergalactic medium.

Finally, it is interesting to note that the quality of the data is becoming good enough to enable a spatially resolved metal budget in galaxies. For instance, Belfiore et al. (2016a) use spatially resolved metallicities for stars and gas combined with spatially resolved maps of the gas content and surface brightness to illustrate that within the central 7 kpc (\(\sim 3~\hbox {R}_{\mathrm{e}}\)) of the well studied galaxy NGC628 about 50% of the metals have been lost (somewhat in tension with the result obtained by Peeples et al. (2014), unless many more metals are lost at larger radii). Interestingly, Belfiore et al. (2016a) also find that the fraction of metals lost increases to about 70% in the central kpc of the galaxy (a similar result was found by, Greggio and Renzini 2011), suggesting that such metals were ejected either by the SNe associated with the early central burst of star formation, associated with the formation of the bulge, or by AGN/quasar driven winds, during the past evolution of the central region of the galaxy.

Very recently, Telford et al. (2018) have performed a very similar, extensive and detailed analysis of the spatially resolved metal budget in M31, finding very similar results, i.e. a higher loss of metals from the central region of M31. Very interestingly, they also find that during the past 1.5 Gyr some of the metal lost from the central region have been redistributed in the galactic disc outskirts.

10 Conclusions

In this paper, we have tried to review the measuring methods, the observational results, and the implications for models of galaxy metallicity evolution. It is the result of many years of effort by many researchers, sometime using dedicated instruments, only part of this effort is reproduced here.

10.1 Summary

The study of chemical abundances in galaxies is a complex and extended field with many open problems and conflicting results. Nevertheless, a few clear points are emerging about methods, observations and models:

  • Stellar metallicities are now routinely measured using UV and optical spectra. Spectrophotometric models of increasing precision, complexity and spectral resolution, use the full information contained in the spectra to derive the metallicity together with other parameters of the stellar population. Simplified methods exist that use particular features to derive metallicities, and these methods are more apt to study large samples of galaxies with lower resolution spectra.

  • The absolute scale of the gas-phase metallicity is still uncertain because of the difference among the three main \(T_{\text {e}}\) method and photoionization models). Discrepancies have been reduced, but they still persist. The “direct” method based on measuring \(T_{\text {e}}\) is currently the most reliable and seems to be in agreement with the metallicity of young stars.

  • These methods have been used to calibrate a large range of strong line ratios diagnostics, which can be applied to faint, distant galaxies. The difference among these secondary calibrations are dominated by the different method used for the primary calibration, photoionization model vs. direct ‘\(T_{\text {e}}\)’ method.

  • Most heavy elements are generally largely depleted onto dust grains; therefore, the evolutions of gas-phase metallicity and of dust are linked together. Dust depletion is a very important, often neglected source of uncertainty in the study of gas-phase abundances.

  • Both stellar and gas-phase metallicities follow a well-defined mass–metallicity relation (MZR) in which the metallicity of galaxies increases with stellar mass. The observed MZR evolves with redshift, with metallicity decreasing at any mass, although more rapidly at lower masses.

  • The gas metallicity of galaxies has also other secondary dependencies. The most important of which is the anti-correlation between gas metallicity and SFR (or gas content, which is related to the SFR), which is called (together with the mass dependence) the Fundamental Metallicity Relation (FMR). This relation has no or a very limited evolution with redshift up to \(z\sim 2.5\). Most authors describe the and a possible strong evolution at \(z\sim 3.5\). FMR as an effect of gas infall, providing further evidence for the ubiquity and importance of cold gas accretion in shaping galaxy evolution, and explain the absence of evolution as the effect of the same dominant physical processes at z < 2.5.

  • Environment also has a secondary effect on the metallicity of galaxies, but only for satellites, whereby satellite galaxies in denser environments (e.g., group and clusters) tend to be more metal rich than galaxies in low-density environments. This is probably a consequence of multiple different effects (such as “strangulation” and accretion of metal-enriched gas).

  • Understanding the redshift evolution of the MZR and of the FMR is made challenging and uncertain by the evident evolution of the ISM properties seen in the excitation BPT diagnostic diagrams. The dominant cause of this evolution is not clear, and most likely a combination of different effects (higher pressure, harder ionizion continua, higher ionization parameter, and variation of the N/O abundance ratio relative to local galaxies). How this evolution affects the determination of metallicities is not yet clear.

  • The metallicity evolution of DLAs provides independent information about the evolution of galaxies and of their CGM. If their velocity dispersion is taken as a proxy of their mass, then the metallicity of DLAs follows a mass–metallicity relation as well. Although it not straightforward to link absorption-selected systems to emission-selected galaxies, both classes of objects identify the redshift range \(2<z<3\) (coincident with the peak of star formation density) as a turning point in galaxy evolution, as this is the redshift range where the evolution of most scaling relations change significantly.

  • The metallicity distribution inside galaxies contains a wealth of information about the spatially resolved processes of galaxy formation. The radial metallicity gradients of local galaxies steepen as a function of galaxy stellar mass (at least within the central \(\sim {2R}_{\mathrm{eff}}\)). However, the outskirts of galaxies show very flat metallicity gradients, which imply accretion of pre-enriched gas from the halo.

  • The scaling relations between galaxy metallicity, stellar mass, star formation rate and gas content are also found locally, on spatially resolved scales, in the form of correlations between metallicity and surface density of stellar mass, surface density of SFR and surface density of gas. However, it is not yet clear whether the local scaling relations are totally driving the global, galaxy-wide ones.

  • Based both on diagnostics that trace the metallicity at different lookback times in local galaxies and the direct observation of metallicity gradients at high redshift, there is a clear indication that the radial metallicity gradients of galaxies have become steeper with cosmic time. This result is difficult to explain in the context of inside-out formation of galaxies and it requires radial migration of stars, radial inflows of low-metallicity gas or radial redistribution of metals at the early stages of galaxy formation.

  • There is growing evidence for the existence of inverted (i.e., positive) metallicity gradients in distant galaxies. These may trace accretion of near-pristine gas or radial inflow of metal-poor gas induced by galaxy mergers or interactions.

  • Different elements provide different information on the evolutionary stage and star formation history of galaxies. Oxygen, Nitrogen, Carbon and the \(\alpha \)-elements are particularly useful and are studied in detail to understand galaxy formation and evolution because they sample different time scales.

  • The \(\alpha \)/Fe abundance ratio shows that the different components of the MW (halo, bulge, thick and thin disc) have formed not only at different epochs, but also on different timescales and with different star formation efficiencies.

  • Similarly, the \(\alpha \)/Fe abundance ratio in local galaxies indicates that more massive galaxies formed faster (probably with higher star formation efficiency) and at earlier times than lower mass galaxies, a phenomenon which is known as galaxy “downsizing”.

  • The N/O and C/O abundance ratios, together with O/H, are used to obtain further information on galaxy evolution as nitrogen and carbon are characterized by longer production timescales than \(\alpha \)-elements, and, therefore, provide important information on the evolutionary stage of galaxies and on other galaxy evolutionary processes (such as gas accretion and efficiency of star formation).

  • X-ray observations have shown that the intracluster medium is highly enriched and with solar-like chemical abundances. The estimated global content of metals in massive clusters exceed significantly the amount of metals that are expected to be produced by the clusters’ galaxies. This is still an unsolved, open problem.

  • AGN are currently one of the few ways to study metallicity in galaxies up to very high redshifts (\(z>7\)). The metallicity information can generally be extracted only for gas in the so-called Broad Liner Region (BLR, on sub-parsec nuclear scales) and for the Narrow Line Region (NLR, on galactic scales).

  • The BLR metallicity scales with black hole mass at all redshifts, which is likely tracing a combination of mass–metallicity relation in the host galaxy and \({M}_{\mathrm{BH}}\)\({M}_{\mathrm{gal}}\) relation, both already in place at early cosmic epochs..

  • The metallicity of the BLR is generally very high (often a few/several times super-solar) and does not evolve with redshift. The latter is probably a consequence of observational selection effects (quasars are detected only when their black holes, and therefore their galaxies, are massive enough, and, therefore, already highly enriched).

  • The metallicity of the NLR is lower than in the BLR, but still higher than in normal galaxies, probably as a consequence of dust destruction and enrichment by quasar-driven outflows. Also the NLR shows little evolution with redshift.

  • Different types of models have been developed to reproduce and explain all these and other observations. The main contribution of the metallicity and chemical abundance analysis is to put strong constraints on the role and the properties of gas infall, galactic winds, stellar and AGN feedback, stellar yields, amount of re-accretion of gas and IMF. Both analytic and numerical models can now account for most of the observational results, but they often disagree on the dominant effects.

  • The distribution of metal mass in the various components of the universe is still not well constrained; nevertheless, it appears that most of the metals have left the parent galaxies.

10.2 Open issues

Although impressive progress has been achieved in recent years, it is clear that several outstanding problems and uncertainties have yet to be faced by observations and theory.

The difficulty of accessing direct gas-phase metallicity tracers of the gas phase for the vast majority of galaxies, especially at high redshift, hence having to rely on strong-line diagnostics or photoionization models remains one of the main issues. Indeed, the latter methods are still prone to degeneracies with other galaxy parameters such as ionization parameter, shape of the ionizing continuum, relative chemical abundances, geometry, pressure and density of the ionized clouds, and the broad distribution of all these parameters inside galaxies; these issues make the comparison between different classes of galaxies and with models still difficult.

Within this context, it is not yet clear how much the evolution of the average excitation conditions of star-forming regions at high redshift (i.e., the evolution of the BPT diagrams) affects the metallicity determination of the gas phase through models or strong line diagnostics calibrated locally.

Similarly, it is still not clear if and how the FMR evolves at high redshift and, if so, why. As discussed, galaxies out to \(z\sim 2.5\) follow the same FMR observed locally once consistent calibrations and a consistent formalism are adopted. However, the dispersion is larger at high redshift than locally and it is not clear why, i.e., whether this is a consequence of observational uncertainties or it truly reflects an evolution of the galaxy evolutionary mechanisms. Even more intriguing is the evolution of the FMR beyond \(z\sim 2.5\), which is accompanied by a similar regime change in other properties of galactic and intergalactic properties around the same redshift (such as the evolution of the content of neutral atomic gas in galaxies, and the evolution of the DLA’s mass–metallicity relation); this change of regime at \(z\sim 2.5\) remains to be fully explained by models.

The redshift evolution of other scaling relations at high redshift, such as with gas content (either molecular or atomic) or with environment, cannot yet be investigated, due to the lack of data on these other properties, jointly to metallicity determinations.

Mainly because of sensitivity limitations, the measurement of the stellar metallicities at high redshift is still limited to very small samples, preventing us to investigate the evolution of scaling relations and the comparison with models at high redshift.

The investigation of metallicity gradients at high redshift is still largely hampered by the lack of spatial resolution. Gravitational lensing, the use of adaptive optics and HST have delivered nicely resolved metallicity maps for small samples, but we are still far from achieving the statistics obtained locally, which has made it possible to investigate metallicity gradients as a function of galaxy properties. Moreover, it is becoming clear that at high redshift metallicity variations within galaxies do not follow the same, simple radial behaviour as in local galaxies, but more complex non-radial variations, implying that a re-thinking of the metallicity gradients characterization is needed.

Models and numerical simulations have made an excellent progress during recent years and can nicely reproduce many of the observed metallicity properties in galaxies, locally and at high redshift. However, all models still suffer from a number of degeneracies and a-priori assumptions that are difficult to control or verify. The IMF (both in terms of shape and cut-offs) is one of the critical input parameters of models, which drastically affects their results. Models obviously depend even more critically on their particular choice of yields and enrichment delay times of stars with different masses. Result from models also depends on the assumed dependence of the outflow loading factor and star formation efficiency as a function of galaxy mass, star formation rate and AGN activity.

Numerical simulations are still limited by lack of high enough resolution to correctly model the subgrid physics associated with star formation and feedback by SNe and AGN. It is observed that the shape of some of the metallicity scaling relations depends on the adopted resolution. Higher resolution simulations can better capture the baryonic physics, but the unavoidably smaller volumes sampled by these simulations result both in a potential bias towards lower density environments and a shortage of massive systems, which may affect the comparison with observations, especially in terms of dependence on environment.

Another outstanding problem consists in the way models and observations are compared. As the actual shape of MZR and FMR depends on how the galaxies are selected, it is necessary to select the model galaxies for the comparison in the same way. This is only rarely done, and usually biased observed samples are compared with volume-limited model samples.

However, major progress on these various fronts is expected in the near future, thanks to the development of new observing facilities and new generation of models, as discussed in the next section.

10.3 Future prospects

In order to properly constrain models, and to advance our understanding of the mechanisms driving galaxy evolution, additional and more accurate observational data are needed.

The next generation of observing facilities and surveys will certainly enable a major step forward in this area of research in the next decade.

The James Webb Space Telescope (JWST) holds some of the major expectations. Its unprecedented sensitivity in the near/mid-IR bands, coupled with its high angular resolution, and multiple spectroscopic modes (including a multi-object spectroscopic mode, integral field units and slitless spectroscopy) will enable astronomers to probe nebular emission lines that are metallicity diagnostics at high redshift, out to the re-ionization epoch and beyond, for several thousands of galaxies, including very low mass systems. Most importantly, JWST will enable astronomers to directly detect auroral lines in hundreds of individual galaxies, hence directly measure the metallicity and recalibrate the strong line diagnostics at different epochs. JWST will also deliver high-fidelity maps of the distribution of metals in hundreds of galaxies. The expectation of JWST is to trace the metallicity evolution of metals back to the first generation of stars, i.e., the so-called PopIII stars, formed out of primordial pristine gas.

On a similar time-frame the Extremely Large Telescopes (such as the Giant Magellan Telescope, GMT, the thirty-meter Telescope, TMT, and the European Extremely Large Telescope, E-ELT), with their huge collecting areas will deliver very high signal-to-noise spectra at intermediate spectral resolution of the stellar continuum in large samples of distant galaxies, therefore, enabling a major leap forward in the characterization of the stellar metallicities, relative chemical abundances (especially \(\alpha \)/Fe) and the associated scaling relations at high redshift. The determination of metallicity gradients will also greatly benefit from the leap in angular resolution delivered by these telescopes, together with their adaptive optics systems. High-resolution spectroscopy, the technique that is most severely affected by photon starving, will probably be the area that will most benefit from the huge collecting area of these telescopes. The number of new systems that will be observable at high spectral resolution will increase by orders of magnitude (thanks to the steep luminosity function of quasars), hence enabling an unprecedented mapping of chemical elements across the universe through absorption systems, with the ultimate goal of finding the chemical signatures of the first generation of stars.

On shorter timescales, the advent of the next generation of high-multiplexing, optical, multi-object spectrographs on 4 m class telescopes, such as WEAVE on WHT and 4MOST on ESO/VISTA will allow to expand the number of observed galaxies to several millions, obtain spectra at higher resolution, and increase the redshift range sampled. Even more interesting, the new near-IR multi-object spectrographs on 8 m class telescopes, such as PFS at Subaru and MOONS at the VLT, will deliver Sloan-like surveys at high redshift by providing near-IR spectra for millions of galaxies out to \(z\sim 2\) and gas metallicities for hundreds of thousands of them. This will enable astronomers to explore the redshift evolution of the metallicity scaling relations with unprecedented statistics. The expected leap in statistics, volume and completeness for distant galaxies will make it possible to investigate, for the first time, the environmental effects on the metallicity scaling relations at high redshift. Stacking of hundreds/thousands of spectra is also expected to enable the investigation of the stellar metallicities and also to detect the auroral lines to recalibrate the strong line diagnostics at high redshift. The current generation of 8m-class telescopes with the next generation of adaptive-optics assisted spectrographs, such as VLT/ERIS, will be used to obtain spatially resolved metallicity for a large number of galaxies and with an higher level of accuracy.

For the warm/hot gas phase the short-lived Hitomi mission has given just a glimpse of the wealth of information that can be obtained through high-resolution high-sensitivity X-ray spectroscopy. The X-ray Imaging and Spectroscopy Mission (XRISM), to be launched in 2021 by JAXA and NASA with a European Space Agency (ESA) participation, will provide observing capabilities similar to the Hitomi satellite and will, therefore, enable us to obtain detailed and accurate measurements of the hot plasma in galaxy clusters, in galactic halos and in galactic winds. This will finally enable astronomers to both derive a much more accurate budget of metals and to directly witness the metal enrichment of the CGM in different classes of galaxies. XARM will pave the way to Athena, the large X-ray observatory to be launched around 2028, which will trace the metallicity and chemical abundance of the hot gas even in distant systems, thanks to its unprecedented sensitivity.

The Atacama Large Millimeter Array (ALMA), which has recently entered in full operation is already delivering exceptional results. In the coming years it is expected to provide a detailed census of the molecular gas in galaxies across the cosmic epochs. The capability of tracing the transitions of multiple molecular species, involving several different elements and associated with different isotopes, will provide unique constraints on the star formation history and on the IMF, both locally and at high redshift.

The Square Kilometre Array (SKA), among many expected ground-breaking results, will finally provide a census of the content and distribution of atomic neutral gas in galaxies at high redshift, which is the key ingredient, still largely unconstrained, to understand galaxy evolution at high redshift and the role played by the HI gas reservoir in distant galaxies.

On longer timescales SPICA space mission, if selected by the European Space Agency (ESA) , will offer a sensitivity improvement by orders of magnitude in the mid- and far-infrared spectroscopic ranges. By measuring fine-structure transitions of several chemical elements, SPICA will enable astronomers to trace the metal enrichment in several hundreds, or even thousands high-z galaxies, without being affected by dust extinction, hence probing, in an unbiased way, also the heavily obscured population of galaxies.