SDSS-V LVM: Mapping of plasma temperatures in the Lagoon Nebula reveals the origin of long-standing chemical abundance discrepancies

Introduction

The chemical composition of galaxies reflect the cumulative history of star formation, stellar nucleosynthesis, astration, supernova explosions, gas and dust chemical redistribution and depletion processes, and the inflow, mixing, and dispersal of gas in and out of galaxies across cosmic time. These processes underpin a wide range of phenomena in astrophysics, making chemical abundances a key parameter in stellar and galaxy evolution models ([45, 36, 44, 27, 8, 28]).

In star-forming galaxies, H II regions mark sites of recent star formation, where radiation from young massive stars ionizes the surrounding neutral gas. The resulting plasma emits a spectrum rich in emission lines from multiple elements, whose relative intensities provide primary diagnostics of the physical conditions – electron density ($n_e$), and electron temperature ($T_e$) – and chemical composition of the ionized gas.

The abundance of heavy elements in the interstellar medium (ISM), measured with respect to that of hydrogen (H), provides a fundamental tracer of galaxy evolution. Oxygen, the most abundant heavy element, is widely used as a proxy for measuring the overall heavy element content of the ISM (e.g. [27, 8, 35]). The strong collisionally excited lines (CELs) of O$^{++}$ are readily detected in the optical spectra of H II regions. Faint auroral CELs, historically hard to detect in many galactic and nearby H II regions, are now detected out to a red-shift, $z \sim10$, thanks to the James Web Space Telescope (JWST), while the extremely faint recombination lines (RLs) of O$^{++}$ can only be reliably measured in the brightest nearby regions ([16, 19, 38, 11, 55, 59]).

Gas-phase O$^{++}$ abundances are commonly derived using the so-called “direct method”. In which, the electron temperature is inferred from ratios of CELs with different excitation energies, while the electron density is obtained independently from ratios of CELs with close-lying upper levels, whose relative populations are set by the density-dependent competition between radiative decay and collisional de-excitation. Given $T_e$ and $n_e$, the emissivities of O$^{++}$ CELs can be calculated directly, allowing the O$^{++}$/H$^+$ abundance ratio to be determined by comparing their brightness to that of hydrogen RLs ([45]).

Because CEL emissivities depend exponentially on $T_e$, the direct method is intrinsically sensitive to the thermal structure of a nebula. Hotter regions contribute disproportionately to the observed emission, while cooler regions contribute comparatively little once collisional excitation becomes inefficient. As a result, if present, unresolved temperature variations along the line of sight can bias the inferred $T_e$ toward higher values, leading to systematically underestimated O$^{++}$ abundances from CELs ([49, 38]).

RLs of elements heavier than Helium (“metals” in astrophysical lingo) provide an alternative abundance diagnostic that is largely insensitive to temperature fluctuations. Because the emissivities of O$^{++}$ and H$^+$ RLs have a very similar temperature dependence, their ratio provides an almost temperature-independent measure of the O$^{++}$ ionic abundance ([45, 50]). However, O$^{++}$ RLs are intrinsically 3 to 4 orders of magnitude fainter than their CEL counterparts, restricting their detection to the brightest and nearest H II regions.

Since the early work of Wyse [65], it has been recognized that abundances derived from heavy-element RLs systematically exceed those obtained from CELs in H II regions, often by factors of up to 3 ([52, 16, 20, 19, 15, 7, 38, 59]). This long standing offset is quantified through the abundance discrepancy factor (ADF), defined as the logarithmic difference between RL- and CEL-based abundances. Despite decades of observational and theoretical effort, this discrepancy has persisted as a fundamental obstacle to the accurate determination of chemical abundances in ionized nebulae for more than 80 years.

Resolving this discrepancy is imperative. If heavy-element abundances in star-forming regions are underestimated by a factor of 2–3, the metal content of galaxies and the cosmic heavy-element budget would be systematically biased ([47]). Such an offset would propagate across a wide range of astrophysical inferences, directly impacting constraints on physical processes associated with metal production (e.g., nucleosynthesis, [28]), ISM enrichment (heavy element stellar yields, stellar winds and supernovae feedback, [44, 24]), the baryon cycle on galactic scales (gas inflows and outflows, [63, 53]) and heavy element depletion (e.g. dust formation and metal capture in stellar remnants [14]).

The metallicity of the ISM further regulates gas cooling and the resulting star formation, shaping the stellar initial mass function both in theoretical models ([12] (Hocuk & Spaans 2010)) and observations ([66, 32]), with further implications for models of the Galactic habitable zone, shifting our expectations of where and when a galaxy like the Milky Way, is most likely to form life-hosting planets like Earth ([33, 26, 37]). Furthermore, since direct-method (i.e. CEL) abundances underpin most of the strong-line diagnostics used to infer chemical abundances in galaxies across cosmic time [7, 27, 35, 57], any systematic bias in the local CEL abundance scale would propagate into our understanding of galaxy chemical evolution [42].

The physical origin of the abundance discrepancy remains debated. Several hypotheses have been proposed, including shocks and turbulence by stellar feedback [51, 6, 54]; chemically enriched high density inclusions embedded within the ionized gas [61, 60], and departures from a Maxwellian electron energy distribution in the ionized ISM [43]. Many of these scenarios invoke spatial variations in the thermal structure of the plasma, commonly described through the temperature fluctuation formalism of [49].

Studies based on slit-spectroscopy of H II regions have provided evidence supporting a connection between the ADF and temperature fluctuations [52, 19, 38], but limit these comparisons to the areas covered by the slits, and do not provide an integral view across the nebula [41]. Spatially resolved observations probing the interior of individual H II regions are therefore essential to test these claims and to place direct observational constraints on the origin of the ADF.

Using deep, ultra-wide-field integral field spectroscopy from the Sloan Digital Sky Survey V (SDSS-V, [29]) Local Volume Mapper (LVM, [13]), Singh et al. [59] mapped the Milky Way H II region M8, better known as the Lagoon Nebula, at a spatial resolution of 0.21 pc per spaxel, enabling the first spatially resolved investigation of the O$^{++}$ abundance discrepancy across an entire H II region. Singh et al. [59] used these data to produce the first spatially resolved IFU maps of O$^{++}$ RLs and the ADF across an H II region, revealing that the discrepancy peaks near the main ionizing source, Herschel 36 (Her36, the brightest component of a trinary system O7.5+ star [2]), and decreases with distance away from it. This behavior is consistent with the slight increase in the ADF reported by Mesa-Delgado et al. [41] in the central $\approx200$” region of the Orion Nebula using slit-spectroscopy.

In this work, we establish the physical origin of the observed abundance discrepancies observed across the Lagoon Nebula. In § 2, we describe the data. Sections 3 outlines the methodology involving measurements of different temperature diagnostics and a simple plasma temperature perturbation model. We present the results in § 4, followed by Discussion in § 5. We summarize the main conclusions in § 6.

Data

We use integral field spectroscopic observations of M8 obtained as part of SDSS-V LVM early science program. The observations and data reduction are described in detail in Singh et al. [59], who presented spatially resolved maps of the physical conditions, O$^{++}$ abundance, and the first spatially resolved map of the ADF, across the entire nebula. Here, we use these measurements to investigate the origin of the abundance discrepancy. To assess radial trends and the effects of spatial averaging, we perform the analysis at three spatial scales: (1) individual spaxel spectra providing fully resolved maps; (2) spectra co-added in 1′ (0.36 pc) wide annuli centered on Her36; and (3) the average M8 spectrum obtained by averaging over the full FOV.

The LVM observations cover the full extent of the nebula, as shown in Figure 1. However, the spaxels satisfying the signal-to-noise ($S/N$) requirements for the CEL and RL-based diagnostics are concentrated within the central part of the observed field. We therefore define the region enclosed by the gray outline in Figure 1 as the FOV shown in the subsequent analysis maps; this region contains all measurements used in this work. For the CEL-based measurements, we require $S/N > 3$ on [O III] $\lambda4363$. For RL-based diagnostics, we use the O II V1 multiplet, consisting of eight lines: 4638.86, 4641.81, 4649.13, 4650.84, 4661.63, 4673.73, 4676.23, and 4696.36 Å. The two faintest lines, $\lambda4673.73$ and $\lambda4696.36$ are not reliably detected in our data, and therefore we use the remaining six lines. For measurements based exclusively on O$^{++}$ RLs, we require $S/N > 3$ for the sum of two brightest V1 lines, $\lambda4649$ and $\lambda4650$. For quantities that combine CEL and RL measurements, both the [O III] $\lambda4363$ and O II V1 criteria are applied simultaneously [59]. All strong CELs and HI lines used in this study have a $S/N > 5$.

Methods

Te and ADF estimation

In this study we make use of two independent electron temperatures diagnostics for O$^{++}$ derived in Singh et al. [59]. The CEL-based (i.e. direct-method) temperature, $T_{e}^{CEL}$(O$^{++}$), is derived from the [O III] $\lambda4363$/[O III] $\lambda(4959+5007)$ ratio, and the hybrid temperature, $T_{e}^{RLs/CEL}$(O$^{++}$), is derived from the ratio of the summed V1 RL flux to [O III] $\lambda5007$ CEL. Because the latter diagnostic combines RLs with a CEL, its temperature dependence is different than that of a purely CEL ratio based diagnostic [48, 59]. Both temperatures are calculated using electron densities from the $S^{+}$ diagnostic based on the [S II]$\lambda\lambda6717$, 6731 doublet: $n_{e}(S^{+})$.

The O$^{++}$/H$^{+}$ abundances from CELs and RLs are derived from the relative intensities of $\lambda4959+5007$ and the summed O II V1 multiplet, respectively, with respect to H$\beta$. For both abundance determinations, we adopt $T_{\mathrm{e}}^{\mathrm{CEL}}(\mathrm{O}^{++})$ and $n_{\mathrm{e}}(\mathrm{S}^{+})$. Using the same physical conditions in both methods ensures that the difference between the resulting abundances reflects the different sensitivities of CELs and RLs, rather than differences introduced by the adopted nebular conditions. In any case, the impact of $T_{\mathrm{e}}$ is negligible on the RL abundance (see Singh et al. [59], for details).

The $n_{\mathrm{e}}$, $T_{\mathrm{e}}$, and ionic abundances are calculated using PyNeb (Luridiana et al. [34]). The atomic data sources adopted are summarized in Table 3 of Singh et al. [59]. Uncertainties are propagated via 50 Monte Carlo realizations, after confirming that higher number of realizations do not significantly change the results. In each realization, the measured line fluxes are perturbed within their measurement errors, and all the derived quantities are recalculated. The final measurements and their uncertainties are the mean and standard deviation of the resulting distributions, respectively.

We quantify the temperature discrepancy as the difference between the temperature diagnostics:

and the abundance discrepancy factor as:

Plasma temperature perturbation modeling

To quantify how electron temperature inhomogeneities along the line of sight bias the observed line fluxes that are used in the different nebular diagnostics used in § 3.1, we calculate theoretical predictions using a simple temperature perturbation model based on the $t^{2}$ formalism of Peimbert [49].

The mean emissivity of a line whose volume-averaged emission arises from a plasma with spatial temperature fluctuations is expanded in a Taylor series around the emission-measure weighted mean temperature $T_{0}$:

where the normalized variance $t^{2} = \langle(T_{\mathrm{e}}-T_{0})^{2}\rangle/T_{0}^{2}$ parametrizes the fractional amplitude of the fluctuations, and the average temperature $T_{0}$ is a weighted average electron temperature over the volume under consideration (in this case the integral along the line of sight over the solid angle sustained by the areas over which the spectra is integrated), weighted by the densities of the particles involved in the emission mechanism with emissivity $\varepsilon$. As the CEL emissivity depends exponentially on $T_{\mathrm{e}}$, hotter regions contribute disproportionately to the observed [O III] $\lambda4363/\lambda(4959+5007)$ ratio, biasing $T_{\mathrm{e}}^{\mathrm{CEL}}(\mathrm{O}^{++})$ systematically above $T_{0}$. The hybrid RL/CEL ratio has a weaker effective temperature sensitivity, so $T_{\mathrm{e}}^{\mathrm{RLs/CEL}}(\mathrm{O}^{++})$ is less susceptible to the temperature fluctuations in the gas than $T_{\mathrm{e}}^{\mathrm{CEL}}(\mathrm{O}^{++})$ (see Figure 2).

Vector (PDF)

The emissivities required to evaluate Equation (3) are given by Peimbert & Peimbert [48]:

Where the normalization constant $C_{\lambda}$ depends only on atomic parameters, including collisional excitation and de-excitation parameters, radiative transition probabilities, and recombination coefficients. For the $\lambda4959$ transition the constant $C_{4959}=1.74\times10^{-17}\,\mathrm{erg\,s^{-1}\,cm^{3}}$ is derived from the collisional excitation equilibrium condition,

where $n_{\mathrm{u}}$ is the population density of $\mathrm{O}^{++}$ in the upper energy level of the $\lambda4959$ transition, $n(\mathrm{O}^{++})$ is the total ion density, and $A_{ul}$ is the Einstein spontaneous transition probability for transition from the upper (u) to lower level (l). The collision strengths and transition probabilities from Lennon & Burke [31] and Osterbrock & Ferland [45]. The ratio $C_{\mathrm{V1}}/C_{4959}=6.56 \times10^{-5}$ is taken from Peimbert & Peimbert [48], while $C_{\mathrm{H}\beta}=2.71 \times10^{-23}$ from Osterbrock & Ferland [45]. The [O III] $\lambda4959/\lambda5007$ ratio is a constant under typical H II regions conditions, and both these lines share the same $T_e$ dependence. All line emissivities are therefore evaluated at $T_e^{\mathrm{CEL}}$, which is calculated using equation (9).

Applying Equation (3) to each line, Peimbert & Peimbert [48] shows that the observed CEL and hybrid temperatures relate to $T_0$ and $t^2$ through:

where $T_0(\mathrm{O}^{++})$ is the average electron temperature weighted by the $n_e$ and the ion number density of $\mathrm{O}^{++}$. The ionic abundances derived from each line family are:

Where $I_{4959}$ and $I_{\mathrm{H}\beta}$ are the intensities of the [O III] $\lambda4959$ and $H\beta$ lines, respectively; and $I_{\mathrm{V1}}$ is the sum of individual O II V1 intensities.

Note the use of $T_e^{\mathrm{CEL}}$ in equation (12). As discussed in § 3.1, the impact of adopting either $T_e^{\mathrm{CEL}}$ or $T_e^{\mathrm{RLs/CEL}}$ when calculating the RL abundance is negligible. This is because the exponents in equations (4) and (7) are almost identical.

Equations 9 and 10 together constrain both $T_0$ and $t^2$ from the measured pair $(T_e^{\mathrm{CEL}}, T_e^{\mathrm{RLs/CEL}})$. We construct a model grid spanning $T_0=5000-12000\ \mathrm{K}$, $t^2=0-0.1$, and intrinsic abundance $12+\log(\mathrm{O}^{++}/\mathrm{H}^{+})=6.5-9.5$ in steps of 0.5 dex, and utilize the above equations to determine the predicted outcome of $T_e^{\mathrm{CEL}}(\mathrm{O}^{++})$, $T_e^{\mathrm{RLs/CEL}}(\mathrm{O}^{++})$, CEL and RL-based $\mathrm{O}^{++}$ ionic abundances, and finally $\Delta T_e$ and ADF at each grid point. Resulting figures from Figure 2 to Figure 4, are analyzed in § 4.

Vector (PDF)

Results

Plasma temperature perturbation model predictions

The upper panel of Figure 2 shows that temperature fluctuations preferentially drive $T_e^{\mathrm{CEL}}(\mathrm{O}^{++})$ toward higher values. Because CEL emissivities depend exponentially on the temperature, hotter regions contribute disproportionately to the [O III] $\lambda4363$ emission relative to the [O III] $\lambda4959,5007$ lines when averaging along the line of sight, biasing the [O III] $\lambda4363/\lambda4959$ ratio towards high values. In contrast, $T_e^{\mathrm{RLs/CEL}}(\mathrm{O}^{++})$ is much less sensitive to local temperature inhomogeneities (lower panel) because it is derived from the ratio of RL fluxes to the CEL flux rather than solely on CELs. Figure 3 shows that thermal fluctuations induce a systematic underestimation of CEL-based $\mathrm{O}^{++}$ abundance (mainly through their impact on the derivation of $T_{\mathrm{e}}$) while they leave the RL-based $O^{++}$ abundance largely unaffected.

Vector (PDF)

Importantly, Figure 4 show that the $\Delta T_{\mathrm{e}}$ increases linearly with $t^{2}$, with only a weak dependence on $T_{0}$. This is a no-table aspect of this work, as we have identified an observational quantity ($\Delta T_{\mathrm{e}}$), which by itself provides an excellent observational proxy of the electron temperature fluctuations amplitude in the ionized ISM.

As it has been long established in the literature, models predicts that the ADF increases with $t^{2}$. The relation is not linear, as it flattens for high $t^{2}$ values, and its normalization is highly dependent on the value of the average underlying electron temperature $T_{0}$ (see Figure 4). Previous ADF studies have simply used this relation to infer a particular value of $t^{2}$ to reproduce an observed ADF. The novel aspect of this work is the introduction of the observable quantity $\Delta T_{\mathrm{e}}$, which, when correlated with the ADF, allows us to test directly whether the observed behavior of a real H II region on spatially resolved scales is consistent with that expected in the presence of temperature fluctuations.

Spatially resolved Te and ADF

We examine the spatially resolved, 1′ (0.36 pc) wide annular-binned, and average measurements of $T_{\mathrm{e}}$ and the $O^{++}$ abundance across the M8 nebula. This allows us to characterize the spatial structure of the nebula and assess how the observed quantities relate to each other and change with spatial averaging. Figure 5 presents spatially resolved maps of the key quantities in this analysis. The RL-based O++ abundances are systematically higher than those derived from CELs across the entire nebula. Their logarithmic difference (the ADF) shows a clear spatial trend, with higher values near the ionizing center and declining outward. The two temperature diagnostics show markedly different spatial behavior. The CELs-based electron temperature $T_e^{\mathrm{CEL}}$(O++) increases toward the ionization center. As will be shown below, this behavior is most likely an artifact associated with enhanced temperature fluctuations, and their effect on the luminosity weighting of CELs toward the central region. In contrast, the hybrid electron temperature $T_e^{\mathrm{RLs/CEL}}$(O++), which combines RLs and CELs, is systematically lower than the CEL-based $T_e$ and remains spatially uniform across the nebula. Since all the diagnostics are based on O++ emission, these differences can not be attributed to the different ionic volumes traced by different ions. It is instead a signal that the gas near Her36 is more thermally inhomogeneous ([59]). The temperatures difference, $\Delta T_e$(O++), peaks in the central regions and decreases outward, closely tracing the spatial profile of the ADF, suggesting a direct connection between the temperature discrepancy and the abundance discrepancy.

The $\Delta T_e$–ADF correlation

Figure 6 presents the ADF measured in each annuli as a function of the temperature discrepancy ($\Delta T_e$). The spatially resolved measurements are also shown as density contours encompassing 50%, 75%, and 90% of the measurements. Regions with higher levels of abundance discrepancy between RL and CEL diagnostics, also show higher levels of electron temperature discrepancies.

The observed ADF and $\Delta T_e$, are highly correlated across every scale we probe. The spatially resolved measurements yield a Pearson correlation coefficient of 0.93 (p-value=$10^{-6}$), confirming a robust point-to-point correspondence throughout the full FOV. The 1′ (0.36 pc) wide annular binned measurements yield the same correlation coefficient of 0.93 (p-value=$4.4 \times10^{-8}$). Measurements from the average spectrum of the entire M8 lie along the same relation ([59]) demonstrating that the correspondence between ADF and $\Delta T_e$ is preserved across all spatial scales probed.

Temperature-fluctuation and the $\Delta T_e$ relation

The observed relation can be directly compared with the predictions of the temperature fluctuation model described in § 3.2. Figure 6 presents the theoretical predictions (pink curves) for the way in which the ADF and $\Delta T_e$ quantities should relate to each other in the presence of electron temperature fluctuations.

The upper axis of Figure 6 provides the mapping between $\Delta T_e$ and $t^2$ (at a fixed $T_0 = 6000$, although the effect of this parameter is very small). Comparing the data in Figure 6 to the model predictions, we observe that the M8 region behaves as expected for a plasma at roughly constant average electron temperature, in the presence of electron temperature fluctuations whose amplitude decreases as one moves away from the central ionizing source.

Although the hybrid $T_{\mathrm{e}}^{\mathrm{RLs/CEL}}(\mathrm{O}^{++})$ is substantially a better electron temperature diagnostic than $T_{\mathrm{e}}^{\mathrm{CEL}}$ alone, it remains subject to measurable biases (see Figure 2). By jointly modeling the impact of $T_{\mathrm{e}}$ fluctuations on both CELs and RLs simultaneously, we can constrain $T_0$ more accurately than either diagnostic in isolation. This approach provides a quantitative constraint on both the average electron temperature, $T_0$, and the amplitude of the temperature fluctuations, $t^2$.

Resolving the abundance discrepancy

The average M8 spectrum yields a best-fit average electron temperature of $T_0 = 6240^{+310}_{-360}\ \mathrm{K}$. This value is consistent, within the uncertainties, with the hybrid $T_{\mathrm{e}}^{\mathrm{RLs/CEL}}(\mathrm{O}^{++}) = 6364 \pm49$, but is $\sim1800\ \mathrm{K}$ lower than the CEL-based $T_{\mathrm{e}}^{\mathrm{CEL}}(\mathrm{O}^{++}) = 8025 \pm67\ \mathrm{K}$ ([59]). The annular binned spectra recover a similar characteristic temperature, with a value of $T_0 = 6020^{+30}_{-20}\ \mathrm{K}$, obtained by fitting all the points assuming a constant temperature model. While this a good representative estimate for the nebula, close inspection of Figure 6 does hint for a small increase of a few hundred K towards the outermost bins. The amplitude of the root-mean-square temperature fluctuations inferred from the average M8 spectrum is $22 \pm1\%$ ($t^2 = 0.050$), decreasing from a central value of $24 \pm0.4\%$ ($t^2 = 0.056$) down to $20 \pm3\%$ ($t^2 = 0.043$) for the outermost annular bin. These values are consistent with $t^{2}=0.045$ reported by García-Rojas et al. [20] for a localized region near the center of M8.

Adopting a value of $T_{0}=6240$ K to derive both RL- and CEL-based $\mathrm{O}^{++}$ abundances from the average spectrum of the nebula, reduces the abundance discrepancy in M8 from $0.47 \pm0.02$ dex to $-0.04 \pm0.02$ dex, effectively solving the abundance discrepancy problem. These findings establish that temperature inhomogeneities as the dominant origin of the abundance discrepancy in M8. They further demonstrate that jointly modeling the effects of temperature fluctuations across multiple diagnostics is the most reliable path to obtain unbiased temperature and abundance measurements.

Discussion

Despite decades of investigation, identifying a direct observational connection between the CEL vs RL abundance discrepancy and its physical origin in H II regions has remained challenging, especially because previous observations have generally lacked sufficient spatial coverage and resolution required to resolve the internal structure of individual nebulae ([16, 20, 48], among others). Mesa-Delgado et al. [41] and Mesa-Delgado & Esteban [40] investigated spatial variations in $\mathrm{ADF}(\mathrm{O}^{++})$ and their relation to the thermal structure using long-slit spectroscopy. Mesa-Delgado et al. [41] reported a weak increase in the ADF toward the ionizing center of Orion, while Mesa-Delgado & Esteban [40] observed the central $\approx200^{\prime\prime}$ region of M8 and found a largely homogeneous ADF distribution. These studies provide important spatial constraints, although the slit-based coverage limited the comparison to the areas sampled by the slit rather than providing an integral view of the nebula.

Using the SDSS-V Local Volume Mapper spatially resolved observations of M8, Singh et al. [59] argued that the abundance discrepancy in this nebula is not driven by density variations or by the presence of distinct chemically inhomogeneous phases (unlike scenarios proposed for PNe, [62, 18, 1, 21]). The authors suggest that $T_{\mathrm{e}}$ inhomogeneities could play a central role in producing the discrepancy.

The interpretation of the abundance discrepancy in terms of $T_{\mathrm{e}}$-fluctuations has been debated for decades ([16, 20, 41, 38, 10]). In particular, the role of $T_{\mathrm{e}}$ fluctuations in biasing direct-method CEL-based chemical abundances has recently been put into question. Chen et al. [10] tested the temperature fluctuation hypotheses in Mrk 71 using far-infrared CELs of $\mathrm{O}^{++}$, finding their infrared abundance consistent with the optical CEL abundance rather than the RL abundance, which they interpreted as ruling out temperature fluctuation in that object. This conclusion has been contested by Méndez-Delgado et al. (2024), who show that it is sensitive to the assumed flux calibration, reddening correction, and electron density, and that alternative treatment of the same data are consistent with the $t^{2}$ formalism.

Our spatially resolved analysis of the LVM data on M8 establishes that the sub-parsec spatially-resolved behavior of the nebular emission inside this H II region, is fully consistent with the temperature fluctuations picture. We present the most direct observational evidence to date for $T_{\mathrm{e}}$ inhomogeneities being the primary driver of the CEL vs RL abundance discrepancy in H II regions. This highlights the value of sub-pc integral-field spectroscopic observations of ionized nebulae, to resolve the micro-physics of interstellar ionized gas emission mechanisms.

The quantity $\Delta T_{\mathrm{e}}$, the difference between direct method and hybrid temperature diagnostics, which can be measured directly from deep H II region spectra, offers an excellent observational proxy for $t^{2}$. The methodology presented here, in which $\Delta T_{\mathrm{e}}$ is correlated with the observed ADFs across different areas inside an H II region, provides a new method to directly constrained not only $t^{2}$, but most importantly $T_{0}$, the average electron temperature of the emitting ionized gas. This is important because the ADF does not uniquely constrain the thermal structure: its relation to $t^{2}$ depends on $T_{0}$. Our results demonstrate that jointly considering ADF and $\Delta T_{\mathrm{e}}$ provides a more direct constraint on the underlying thermal structure, which can be used to reliably recover an abundance scale consistent between CELs and RLs.

As wide-field integral field spectroscopy with facilities such as SDSS-V LVM extends resolved coverage to larger samples of H II regions across a range of ionizing conditions and metallicities, $\Delta T_{\mathrm{e}}$ measurements can provide the statistical calibration of the $t^{2}$ correction the field has lacked, moving the Peimbert correction (Peimbert [49]) from a case-by-case assumption to an empirically grounded tool that could be applied generally.

A broader consequence of our result concerns the overall metal budget of galaxies, known as the “missing metals problem”, which arises from the discrepancy between the cosmic density of metals expected to have been produced by previous generations of star formation and that accounted for in stars, the ISM, the circumgalactic medium (CGM), and intergalactic medium (IGM). Current estimates suggesting 35-90 % of metals remain unaccounted for ([17, 9, 46]). Surveys such as COS-Halos (Peples et al. [47]) find that only about 20 − 25% of the metals produced over a galaxy’s lifetime are retained in stars, ISM gas, and dust, with the remainder presumed to reside in the CGM and IGM. The ISM component of this census rests on gas-phase metallicities derived from CELs. If temperature fluctuations of the kind characterized here cause CEL-based methods to systematically underestimate the true ISM oxygen-abundance, then the metal budget of galaxies would likewise be underestimated. Accounting for these effects may partially alleviate this missing metal/oxygen problem. The consequences of resolving the abundance discrepancy extend from individual galactic H II regions to the cosmological accounting of metal production and retention in galaxies.

Along the same lines, a direct implications of this work concerns extragalactic metallicity measurements. The JWST has enabled auroral line detections at redshifts approaching $z \sim10$, allowing the direct method to be applied to galaxies across all cosmic times, and the empirical calibration of strong-line metallicity diagnostics that are applicable at different redshifts ([55, 30, 57, 58]). These galaxies are, however, entirely unresolved; their integrated spectra sum together numerous H II regions spanning a range of ionization parameters, radiation field hardnesses, densities, and temperatures, none of which are individually separable from the integrated light ([56]). In principle, the resulting metallicities from these methods trace the luminosity weighted metallicity of the individual H II regions. But the biases associated with the thermal structure of these regions are not accounted for and should affect these diagnostics, biasing our estimates of the early chemical evolution of the universe.

The thermal structure of H II regions depends on the clumpiness of the neutral ISM, the hardness of the ionizing radiation field, the gas metallicity, and the local pressure, all of which evolve substantially with redshift. While the magnitudes of the temperature fluctuations we measure in M8 may differ substantially from those in galaxies at cosmic noon or in the epoch of reionization, their impact on the emission mechanisms constrained here is the same one that is currently ignored in galaxy-scale CEL-based metallicity measurements. Establishing how the thermal structure of H II regions affects chemical abundance diagnostics in a resolved nearby systems such as M8, is a necessary first step before corrections for these effects can be properly propagated into high-$z$ chemical abundance diagnostic calibrations.

Conclusion

Using spatially resolved SDSS-V Local Volume Mapper observations of M8, we investigate the origin of the longstanding abundance discrepancy between CELs and RL. Our main conclusions are:

1. We find a tight spatial correlation between the $O^{++}$ ADF and $\Delta T_e$. The agreement of this observed relation with the theoretical predictions of the temperature-perturbation formalism, shows that $T_e$ fluctuations are the primary drivers of the CEL–RL abundance discrepancy in M8. The electron temperatures inferred from the direct method are biased toward higher temperatures, leading to a systematic underestimation of the oxygen abundance when temperature fluctuations are ignored.

2. The observed $\Delta T_e$–ADF relation, together with the temperature-perturbation model, provides an empirical way to constrain the thermal structure of an H II region. In particular, $\Delta T_e$ increases linearly with the square of the amplitude of the temperature fluctuations ($t^2$), with only a weak dependence on the underlying average electron temperature ($T_0$). It therefore serves as a direct observational proxy of $t^2$ and a means of constraining the thermal structure of the nebula.

3. The mean temperature, $T_0$, constrained from the $\Delta T_e$–ADF relation, reconciles the CEL- and RL-based O$^{++}$ abundances in M8, effectively solving the abundance discrepancy problem. This demonstrates that accounting for the thermal structure of the plasma in ionized nebulae is essential for obtaining unbiased CEL-based abundances.

The broader significance of these results is that they establish the thermal structure of H II regions as an important factor in assessing the reliability of nebular abundance measurements. Determining whether the temperature fluctuations identified in M8 are common across different ISM environments is therefore a critical next step. The After-Sloan 5 (AS5) Local Volume Mapper 2 (LVM2) survey will extend this analysis to hundreds of H II regions across the Milky Way and nearby galaxies, providing the faint-line measurements needed to systematically measure temperature fluctuations and RL-based abundances across a range of physical conditions and metallicities. M8 thus provides a resolved benchmark for testing how the thermal structure of ionized gas shapes abundance measurements across the local Universe.

Acknowledgments

A.S. acknowledges the support from the ANID Basal projects FB21250718. G.A.B. and E.J. acknowledge the support from the ANID Basal project FB210003. C.M. acknowledge the support of grant UNAM/DGAPA/PAPIIT IG101224 and IN113226. J.E. M-D. and C. M. thank the support by SECIHTI CBF-2025-I-2048 project “Resolviendo la Física Interna de las Galaxias: De las Escalas Locales a la Estructura Global con el SDSS-V Local Volume Mapper” (PI: Méndez Delgado). KK and EE gratefully acknowledge funding from the European Research Council’s starting grant ERC StG-101077573 (“ISM-METALS”). OE acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project-ID 541068876. J.G.F-T gratefully acknowledges the support provided by ANID Fondecyt Regular No. 1260371.

Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, the National Science Foundation, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. SDSS telescopes are located at Apache Point Observatory, funded by the Astrophysical Research Consortium and operated by New Mexico State University, and at Las Campanas Observatory, operated by the Carnegie Institution for Science. The SDSS web site is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Carnegie Institution for Science, Chilean National Time Allocation Committee (CNTAC) ratified researchers, Caltech, the Gotham Participation Group, Harvard University, Heidelberg University, The Flatiron Institute, The Johns Hopkins University, L’Ecole polytechnique fédérale de Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania State University, Smithsonian Astrophysical Observatory, Space Telescope Science Institute (STScI), the Stellar Astrophysics Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Illinois at Urbana-Champaign, University of Toronto, University of Utah, University of Virginia, Yale University, and Yunnan University.

Facilities: SDSS-V Local Volume Mapper

Software: astropy ([3, 4, 5]), Scipy ([64]), Numpy ([22]), Matplotlib ([25]),