# A geostatistical analysis of multiscale metallicity variations in galaxies -
III. Spatial resolution and data quality limits

Benjamin Metha1,2,3, Michele Trenti1,2, Andrew Battisti2,5, Tingjin Chu4
1School of Physics, The University of Melbourne, VIC 3010, Australia
2Australian Research Council Centre of Excellence for All-Sky Astrophysics in
3-Dimensions, Australia
3Department of Physics and Astronomy, University of California, Los Angeles,
430 Portola Plaza, Los Angeles, CA 90095, USA
4School of Mathematics and Statistics, The University of Melbourne, VIC 3010,
Australia
5Research School of Astronomy and Astrophysics, Australian National
University, Cotter Road, Weston Creek, ACT 2611, Australia
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Geostatistical methods are powerful tools for understanding the spatial
structure of the metallicity distribution of galaxies, and enable construction
of accurate predictive models of the 2D metallicity distribution. However, so
far these methods have only been applied to very high spatial resolution
metallicity maps, leaving it uncertain if they will work on lower quality
data. In this study, we apply geostatistical techniques to high-resolution
spectroscopic maps of three local galaxies convolved to eight different
spatial resolutions ranging from $\sim 40$pc to $\sim 1$ kpc per pixel. We fit
a geostatistical model to the data at all resolutions, and find that for
metallicity maps where small scale structure is visible by eye (with $\gtrsim
10$ resolution elements per $R_{e}$), all parameters, including the
metallicity correlation scale, can be recovered accurately. At all resolutions
tested, we find that point metallicity predictions from such a geostatistical
model outperform a circularly symmetric metallicity gradient model. We also
explore dependence on the number of data points, and find that $N\gtrsim 100$
spatially resolved metallicity values are sufficient to train a geostatistical
model that yields more accurate metallicity predictions than a radial gradient
model. Finally, we investigate the potential detrimental effects of having
spaxels smaller than an individual Hii region by repeating our analysis with
metallicities integrated over Hii regions. We see that spaxel-based
measurements have more noise, as expected, but the underlying spatial
metallicity distribution can be recovered regardless of whether spaxels or
integrated regions are used.

###### keywords:

ISM:structure, galaxies:ISM, galaxies:abundances, methods:statistical,
software: data analysis, ISM: abundances

††pubyear: 2023††pagerange: A geostatistical analysis of multiscale
metallicity variations in galaxies - III. Spatial resolution and data quality
limits–C

##  1 Introduction

The distribution of metals (chemical elements heavier than Helium) throughout
a galaxy is sensitive to many details of their formation and evolution that
are not yet fully understood. It is well established that metals are released
into the interstellar medium (ISM) of galaxies by supernovae (Andrews et al.,
2020), but we do not know how efficiently they can mix into their local
environment (Li et al., 2021). We know that a fraction of the metals released
in these supernovae ejecta escape the galaxy via metal enriched outflows (e.g
Heckman et al., 1990; Reichardt Chu et al., 2022), but we do not know how
large this escape fraction is, nor how much of the ejected material falls back
onto the galaxy later in a galactic fountain (Christensen et al., 2018). We
know that galaxies accrete metal-poor gas to fuel star formation (Larson,
1972; Trapp et al., 2022), but there is still some debate as to whether this
happens preferentially on the central regions or the outskirts of galactic
discs (Schönrich & McMillan, 2017). Detailed analysis of the internal chemical
structure of galaxies can shed light on the relative strengths and timescales
of all of these processes, making it a key element for addressing many
theoretical challenges in understanding galaxy formation and evolution (Naab &
Ostriker, 2017).

Traditionally, the internal metallicity structure of galaxies has been
understood through the lens of a radial linear gradient model (Searle, 1971).
Through a large collection of observational campaigns (e.g. Vila-Costas &
Edmunds, 1992; Zaritsky et al., 1994; Croom et al., 2012; Sánchez et al.,
2012; Berg et al., 2013; Ho et al., 2015; Bryant et al., 2015; Bundy et al.,
2015; Belfiore et al., 2017, to name a few), it has been established that most
spiral galaxies in the local universe follow a negative metallicity gradient,
with more metal enrichment in their central regions than in the outskirts of
their disks. This can be explained through an “inside-out" galaxy formation
model, in which star formation begins first in the central regions of a
galaxy, allowing those regions to be more metal-enriched through successive
generations of star formation than the outskirts, where star-formation has
only recently begun (e.g Boissier & Prantzos, 1999).

Recent technological advancements in high resolution integral field
spectroscopy (IFS) has led to acquisition of detailed metallicity maps for
galaxies in the local universe, with resolutions finer than $\sim 100$ parsec
(Erroz-Ferrer et al., 2019; López-Cobá et al., 2020; Emsellem et al., 2022).
Using data from MUSE, and psuedo-IFS data obtained through the TYPHOON survey
(Seibert et al. in prep.), large azimuthal variations in metallicity have been
seen to exist in many local spiral galaxies in addition to a metallicity
gradient (e.g. Sánchez-Menguiano et al., 2016, 2018; Ho et al., 2017, 2018,
2019; Kreckel et al., 2019). This data has also been used to show that small-
scale variations in metallicity exist in addition to a radial trend, and to
analyse the spatial scale of these variations, using a variety of statistical
techniques (Kreckel et al., 2020; Li et al., 2021, 2022; Williams et al.,
2022).

In the first two papers of this series (Metha et al. 2021, hereafter M21; and
Metha et al. 2022, hereafter M22), we have explored how tools and techniques
from geostatistics, the mathematical study of stochastic processes that occur
over a continuous spatial domain, can be applied to analyse data products from
high-resolution studies of metallicity variations within star-forming galaxies
in the local universe. In M21, seven galaxies from the PHANGS-MUSE survey were
analysed using a geostatistical approach to recover the statistical properties
of their small-scale metallicity fluctuations, allowing predictions of an
analytical model to be tested. In M22, geostatistical models were fit to high-
resolution metallicity maps of eight large spiral galaxies observed as a part
of the TYPHOON survey, producing more accurate predictive models of the
metallicity distribution throughout these galaxies that can be applied to
regions obscured by diffuse ionised gas (DIG) contamination where the
metallicity cannot be directly observed.

Despite these successes, a wealth of open questions remain. What are the main
astrophysical sources responsible for these small-scale metallicity
variations? Are they present in all galaxies, or only in large star-forming
spirals in the local universe? Is the size of these small-scale variations set
by the size ($R_{e}$; effective radius) of the galaxies? Do these small-scale
metallicity variations also exist in passive galaxies, or dwarf galaxies? Are
star-forming galaxies with strong spiral arms expected to be more or less well
mixed than those without (i.e. grand design vs flocculent)? Are these small-
scale variations stronger or weaker at cosmic noon, where galaxies are more
clumpy and irregular and often show flat or inverted metallicity gradients?
Finally, does the two-dimensional chemical structure of these galaxies match
what is predicted by hydrodynamical simulations such as IllustrisTNG
(Pillepich et al., 2018) or FIRE (Wetzel et al., 2023)?

To answer these questions and others, it is important to extend these spatial
statistical analyses to a larger sample of real and simulated galaxies, with a
variety of sizes, morphologies, star formation rates, and redshifts. A large
wealth of legacy IFS data exists (Sánchez, 2020) for which geostatistical
analysis could potentially be applied. However, the application of
geostatistical techniques to understand the ISM of galaxies has at present
only been explored for high-resolution IFS (or IFS-like) data products.

Resolution has been identified as the most important figure of merit in the
design of an IFS survey, particularly when the goal is to capture detailed
information about sharp structures that vary over small scales (Mast et al.,
2014). Therefore, it is important to understand how a geostatistical analysis
performs on data at a range of resolutions. This will provide useful
information for the design of future surveys and numerical simulations if
their aims include detailed analysis of the small-scale processes that
determine ISM metallicity.

Similar studies have been done in the past exploring resolution limits in IFS
surveys in the context of fitting metallicity gradients to galaxies. By
rebinning IFS data and simulating observations with poor angular resolution,
Yuan et al. (2013) found that metallicity gradients appear flatter when
observed with low-resolution IFS (FWHM $\gtrsim 340$ pc for a system at
z=1.49). These effects were further investigated and quantified by Acharyya et
al. (2020) using simulated high-resolution IFS datacubes. They found that
input metallicity gradients are generally recovered accurately when the FWHM
of the instrumental PSF is smaller than $\sim 0.25R_{e}$.

In this work, we investigate the impact of spatial resolution by following the
methodology of Yuan et al. (2013) and Mast et al. (2014). Specifically, we
rebin high-resolution IFS data of three galaxies observed by the TYPHOON
survey to increasingly coarse resolutions, in order to assess how the results
of geostatistical analysis are sensitive to resolution. Further, we explore
how the quality of model fits depends on the number of metallicity data
points. Finally, we provide a first investigation into the quantitative
effects of aperture bias by comparing models fit using unbinned spaxels
associated with Hii regions to models fit using integrated Hii region data.

The structure of this paper is as follows. We discuss the data and how it is
processed to produce low-resolution maps in Section 2. In Section 3, we
discuss our metallicity diagnostics, our methods for finding Hii -dominated
spaxels and Hii regions, and the geostatistical model and fitting methods that
we use. We present our results for NGC 5236 in detail in Section 4. We show
how the ability of our geostatistical pipeline to precisely and accurately
recover details on the multiscale metallicity structure of these galaxies is
affected by the number of available data points in Section 4.1, and by
resolution in Section 4.2. In Section 4.3, we explore the trade-off between
using unbinned Hii spaxels, which may be affected by aperture bias, and binned
Hii regions, which are fewer and necessarily at a lower spatial resolution
than individual spaxels. To ensure our findings are not simply peculiarities
of one galaxy, we repeat our analysis for two other spiral galaxies in Section
5. We discuss overall results in Section 6, and summarise the main conclusions
in Section 7, including general advice on when and how these geostatistical
methods can be effectively applied to astronomical datacubes. Finally, we
present our analysis based on alternate metallicity diagnostics in Appendix A,
a primer on some geostatistical techniques and concepts introduced earlier in
this series in Appendix B, and a geostatistical explanation for how the
variance of the metallicity measured in an IFS datacube depends on its
resolution in Appendix C.

##  2 Data

We select NGC 5236 (M83) as our fiducial galaxy in this study, using high
resolution spectroscopic data from the TYPHOON survey Poetrodjojo et al.
(2019). To ensure these results are valid for a range of galaxies, we also
include data from two other galaxies in the TYPHOON survey, NGC 6744 and NGC
7793 (Seibert et al. in prep.). We summarise the properties of these galaxies
in Table 1.

Name | NGC 5236 | NGC 6744 | NGC 7793
---|---|---|---
RA | 13:39:55.96 | 19:09:46.10 | 23:57:49.83
DEC | -29:51:55.5 | -63:51:27.1 | -32:35:27.7
Distance (Mpc) | 4.9 | 11.6 | 3.6
$R_{e}$ (kpc) | 3.3 | 7.2 | 1.8
$R_{25}$ (kpc) | 9.18 | 33.7 | 4.89
Position angle | $54.0^{\circ}$ | $13.7^{\circ}$ | $94.4^{\circ}$
Inclination | $12.5^{\circ}$ | $53.5^{\circ}$ | $55.6^{\circ}$
Morphology | Sc | Sbc | Scd
Stellar Mass $(\log_{10}M/M_{\odot})$ | 10.4 | 10.9 | 9.2
Star formation rate $(M_{\odot}/$yr) | 4.0 | 4.0 | 0.25
Table 1: Properties of the three galaxies explored in this paper. RA and Dec
values are taken from the NASA/IPAC Extragalactic Database. Morphological
classification and position angle data are taken from HyperLEDA. Inclinations
are calculated based on the axis ratios reported by HyperLEDA, assuming a disc
scale height of $q_{z}=0.2$ for all galaxies, following Battisti et al.
(2017). Distances, stellar mass estimates, and SFR estimates are taken from
Leroy et al. (2019). Half-light radii, also known as effective radii $(R_{e})$
are taken from the $K_{S}$ band effective radii measured from the 2MASS large
galaxy atlas (Jarrett et al., 2003). B-band isophotal radii ($R_{25}$) are
taken from de Vaucouleurs et al. (1991). Both radii have been converted to
physical units based on the adopted distance values.

Datacubes for all three of these galaxies were constructed using a Progressive
Integral Step Method (PrISM) using the 2.5m du Pont telescope located at the
Las Campanas Observatory in Chile (Grasha et al. 2022, Seibert et al. in
prep.). A series of adjacent exposures were taken with a long slit of width
$1.^{\prime\prime}65$, which were re-binned to psuedo-IFS datacubes with a
native angular resolution of $1.^{\prime\prime}65$. TYPHOON observations were
only obtained when the seeing was smaller than the width of the long slit. We
note that Las Campanas Observatory has a median seeing of $0.^{\prime\prime}6$
(Persson et al., 1990).

For NGC 5236 at a distance of $4.89$ Mpc (Leroy et al., 2019), an angular size
of $1.^{\prime\prime}65$ corresponds to a physical pixel size of $39$ pc. Each
datacube was rebinned to several additional resolutions by binning blocks of
$f\times f$ adjacent spaxels together. We label each map according to its bin
factor $f$. In this study, bin factors of $1,2,4,6,8,10,12$ and $24$ are
tested, corresponding to physical resolutions from $39$ pc to $940$ pc for NGC
5236. This selection of resolutions covers the range of physical resolutions
from MUSE observations of local galaxies, to high resolution zoom-in
cosmological simulations such as auriga ($125$ pc, Grand et al. 2017), to IFS
observations of gravitationally lensed galaxies with the NIRISS instrument of
JWST ($\sim 200$ pc, Wang et al. 2022), up to large IFS surveys of galaxies
with coarse spatial resolution such as CALIFA ($\sim 800$ pc, Espinosa-Ponce
et al. 2020). Therefore, the range of spatial resolution considered allows us
to explore for which, if any, of these kinds of data products and
observational surveys a geostatistical analysis is feasible.

To illustrate what features of the galaxy are visible, in Figure 1 we show
H$\alpha$ images for NGC 5236 at each resolution explored. At the native
resolution, individual Hii regions can be seen as “bubbles" along the spiral
arms of this galaxy, and regions of low H$\alpha$ emission in the inter-arm
regions are visibly resolvable. As $f$ is increased, structural details of
this galaxy are lost. At $f=4$ (with a spatial resolution of $\sim 150$ pc),
individual Hii regions begin to blur together, but variations can still be
seen along spiral arms. These small-scale variations become difficult to see
at $f=8$, but the azimuthally-varying spiral structure is still visible. At
the coarsest resolution ($f=24$), no azimuthal structure is apparent, and a
radial trend of decreasing H$\alpha$ brightness is all that can be seen.

Figure 1: Surface brightness of H$\alpha$ for NGC 5236 at its native
resolution (top left), and for spatially binned datacubes containing the total
emission from $f\times f$ spaxels in each pixel. For the very high resolution
maps ($f=1-2)$, individual Hii regions may be resolved clearly. On larger
scales ($f=4-6$), bright regions are blurred together, but small features are
still visible. For intermediate resolution maps ($f=8-12$), small scale
features become harder to see, but large-scale features such as spiral arms
are still visible. With low-resolution IFU observations $(f=24)$, all features
apart from a radial trend in H$\alpha$ brightness are lost.

##  3 Methods

For each rebinned datacube, lzifu (Ho et al., 2016) was used to find the
intensities of 13 strong emission lines, including H$\alpha$, H$\beta$,
[Oiii]$\lambda 5007$, [Sii]$\lambda\lambda 6716,31$, and [Nii]$\lambda 6583$.
The strength of all emission lines are corrected for dust attenuation assuming
Case B recombination, fixing the Balmer decrement to be
H$\alpha/$H$\beta=2.86$, using the extinction curve of Cardelli et al. (1989),
assuming $R_{V}=3.1$ to match the Milky Way value.

###  3.1 Hii spaxels selection

Strong emission line ratio based metallicity diagnostics are only valid for
Hii regions – regions recently ionised by hot, bright, young O/B stars. We
ensure that all spaxels used in this analysis are predominately ionised by Hii
region emission at each resolution by imposing data quality cuts designed to
isolate Hii regions from diffuse ionised gas (DIG), a significant source of
contamination for metallicity studies (e.g. Reynolds, 1990; Poetrodjojo et
al., 2019; Kumari et al., 2019, M22).

We use two BPT diagram-based cuts (Baldwin et al., 1981), restricting our
analysis to all spaxels that lie below both the Kewley et al. (2001)
demarcation line on the [Sii]-based BPT diagram and the Kauffmann et al.
(2003) demarcation line on the [Nii]-based BPT diagram. To ensure these cuts
are robust, only spaxels with measurements of H$\alpha$, H$\beta$,
[Oiii]$\lambda 5007$ [Sii]$\lambda\lambda 6716,31$, and [Nii]$\lambda 6583$
fluxes with S/N$>3$ are considered. Following Zhang et al. (2017), we
additionally impose a surface brightness cut, including only spaxels with a
H$\alpha$ surface brightness greater than $10^{39}$ erg s-1 kpc-2. After these
cuts, for NGC 5236, 7488 spaxels are classified as being predominately ionised
by Hii emission at native resolution. We hereafter refer to such spaxels as
Hii spaxels. As $f$ is increased to simulate coarser resolution observations,
the number of Hii spaxels recovered by this methodology decreases. We plot the
number of Hii spaxels available for analysis as a function of the resolution
of the binned datacubes for all three galaxies analysed in this work in Figure
2.

Figure 2: Number of Hii spaxels as a function of resolution. For each galaxy,
as the resolution of the data is decreased, the number of spaxels classified
as Hii dominated also decreases.

From observations of Hii regions in the Milky Way, it has been established
that the spatial extent of a Hii region can range from tens of parsecs to kpc
scales (Cosens et al., 2018), with typical sizes ranging from 80-800 pc
(Espinosa-Ponce et al., 2020). Therefore, the high resolution metallicity maps
examined in this work contain Hii spaxels that are significantly smaller than
an individual Hii region. This may introduce bias into measurements based on
individual high-resolution spaxels, as emission line ratios are not constant
throughout a Hii region (e.g. Ferland et al., 2013), and strong emission line
based metallicity diagnostics are calibrated to work on integrated Hii regions
(e.g. Mannucci et al., 2021).

To assess whether this bias is significant, we additionally extract
metallicities for integrated Hii regions detected using the publically-
available code pyHIIextractor (Lugo-Aranda et al., 2022). We choose this
program over other Hii region detection software because it does not require a
license, it is actively maintained, and it was designed for optimal
performance on high-resolution IFS data. Based on the H$\alpha$ intensity map
and its uncertainty, positions and radii of hundreds of Hii region candidates
are extracted from each galaxy at native resolution. The strengths of the
emission lines for each candidate Hii region are then computed by summing all
(extinction-corrected) fluxes from spaxels contained within each region.

One limitation of pyHIIextractor for our application is that the same spaxel
can be contained within two (or more) different overlapping extracted regions.
In these cases, we assign spaxels to the Hii region whose centre is closest to
them. This step is necessary to prevent spurious correlations arising between
neighbouring Hii regions that are not physical, but merely statistical
artefacts from double-counting.

We plot all Hii region candidates on the two BPT diagrams, discarding any
regions lying above the Kewley et al. (2001) demarcation line on the
[Sii]-based BPT diagram and the Kauffmann et al. (2003) demarcation line on
the [Nii]-based BPT diagram, leaving a final sample of Hii regions for each
galaxy. We note that this processing step decreases both the number of data
points available for a metallicity model fit and the effective resolution of
the data, but may be appropriate in some cases to mitigate against aperture
bias. We explore how including this processing step affects a geostatistical
model fitting procedure in Section 4.3.

###  3.2 Metallicity diagnostics

To ensure that our results are robust, we investigate using four different
metallicity diagnostics – the N2S2H$\alpha$ diagnostic with the calibration of
Dopita et al. (2016), the Scal diagnostic of Pilyugin & Grebel (2016), and the
O3N2 and RS32 diagnostics of Curti et al. (2020). When using each diagnostic,
we discard any Hii spaxels (or regions) that have S/N$<3$ in any of the lines
used. This final data quality cut excludes $<1\%$ of the data. For all of the
diagnostics we explore, we see similar trends for how the geostatistical
methods tested are affected by resolution and data quality. We choose
RS32$=\log($[O iii]$\lambda 5007/$H$\beta+$[S ii]$\lambda\lambda
6717,6731$/H$\alpha$) as our fiducial diagnostic, presenting results from this
diagnostic in the main text. This diagnostic is calibrated using direct
electron temperature-based metallicity estimates of stacked spectra of SDSS
galaxies with a 4th order polynomial fit for $Z$ as a function of RS32. We
refer the reader to Appendix A for details on the other metallicity
diagnostics investigated in this work, as well as for a presentation of the
main results for NGC 5236 when these alternative metallicity diagnostics are
used.

###  3.3 Geostatistical model fitting

In this Section, we introduce the hierarchical geostatistical model that we
fit to our data, and the parameter determination procedure. The model is
identical to the one introduced in M22. However, the fitting procedure has
been updated and improved. For a summary on key background on other tools from
geostatistics that are used in this work and have beenintroduced in the
earlier papers of this series, we refer the reader to Appendix B.

For each data set that we test, we consider the observed metallicity $Z_{\rm
obs}$ at each location $\vec{x}$ to be a combination of the true underlying
metallicity at that point ($Z(\vec{x})$), plus some uncorrelated, stationary,
zero-mean observational noise $\epsilon(\vec{x})$:

$Z_{\rm obs}(\vec{x})=Z(\vec{x})+\epsilon(\vec{x}).$ (1)

We further separate $Z(\vec{x})$ into two components: the process mean
$\mu(\vec{x})$, which captures the large-scale metallicity trends of the
galaxy, plus a random component $\eta(\vec{x})$ with zero-mean, which
represents real, small-scale metallicity fluctuations around this mean:

$Z(\vec{x})=\mu(\vec{x})+\eta(\vec{x}).$ (2)

As in M21, M22, $\mu(\vec{x})$ is modelled as being linearly dependant on the
distance between $\vec{x}$ and the centre of the galaxy. To fully specify this
model, two parameters must be fit: $\nabla Z$, the radial metallicity
gradient, and the characteristic metallicity $Z_{\text{char}}$ at a user-
defined radius $r_{\text{char}}$:

$\mu(\vec{x})=Z_{\text{char}}+\nabla
Z\cdot\left(r(\vec{x})-r_{\text{char}}\right).$ (3)

In M21, M22, and in most studies when metallicity gradients are fit,
$r_{\text{char}}$ is chosen to be zero. $Z_{\text{char}}$ therefore becomes
the central metallicity of the galaxy. However, we found that such a choice of
parameterisation led to strong anticorrelation between $Z_{\text{char}}$ and
$\nabla Z$, with $\rho<-0.9$ between these two parameters irrespective of the
galaxy under investigation, the metallicity diagnostic used, or the resolution
explored. Such a strong anticorrelation has been noted before in the
literature, e.g. by Zaritsky et al. (1994).

To mitigate this anticorrelation, and in turn improve the efficiency of
convergence using Monte Carlo methods, we parameterise the linear function by
fitting the characteristic metallicity at a radius of
$r_{\text{char}}=0.4R_{25}$ (see Zaritsky et al. 1994). The metallicity at
this radius has been observed to approximately match the integrated
metallicity of face-on spiral galaxies using a variety of metallicity
calibrations (e.g Moustakas & Kennicutt, 2006; Rosales-Ortega et al., 2011).
We found that such a parameterisation significantly lowered the degree of
correlation between $Z_{\text{char}}$ and $\nabla Z$, for all fits of our
galaxy maps.

We parameterise the covariance of the small-scale, random metallicity
fluctuations $\eta(\vec{x})$ as an exponential function, following M22. Such a
functional form was chosen because in the small-scale limit, the structure of
this field will mirror that of a passive scalar driven by Burger’s turbulence
(Lee & Gammie, 2021), which is a good model for shock-induced structures (Bec
& Khanin, 2007; Falceta-Gonçalves et al., 2014), and may be the most suitable
paradigm of turbulence for describing the supernova shock-driven mixing of
passive scalars around Hii regions. This function has two parameters: $\phi$,
which sets the distance at which random metallicity fluctuations are
correlated, and $\sigma^{2}$, the total variance about the process mean. For
two data points $\vec{x},\vec{y}$ separated by a distance of $h$, we model the
covariance between the small scale metallicity fluctuations to be:

$\text{Cov}(Z(\vec{x}),Z(\vec{y}))=\sigma^{2}\exp\left(-\frac{h}{\phi}\right).$
(4)

All four parameters of this hierarchical model ($Z_{\text{char}},\nabla
Z,\phi$ and $\sigma^{2}$) are fit simultaneously using a maximum likelihood
method, optimised using the Python package emcee (Foreman-Mackey et al.,
2013). We refer the reader to Equation 6 of M22 for the full form of the 4D
likelihood function that is maximised. To help distinguish between values of
$\sigma^{2}$ that are close to zero, and to ensure that the value of
$\sigma^{2}$ was always positive, $\log\left(\sigma^{2}\right)$ was fit
instead of $\sigma^{2}$, with flat priors.

For every data set explored, we estimate posterior distributions for each
parameter by running emcee for 480 steps using 72 walkers. Exploratory work
showed that the autocorrelation time, $\tau_{AC}$, for each parameter is $\sim
40$ steps. Based on this, 480 steps was chosen such that the program could run
for $10\tau_{AC}$ steps, after discarding $2\tau_{AC}$ steps for burn-in. This
ensures there are effectively $\sim 720$ independent samples from the
posterior distribution of each parameter, which is sufficient for estimating
parameters with realistic uncertainties for this 4 dimensional problem. Each
Monte Carlo chain was examined by hand to ensure that this burn-in length was
long enough; for cases where it was not, an additional 40 steps were ran so
that $3\tau_{AC}$ steps could be discarded.

For each run, walkers were initially started at a wide range of initial values
of $\phi$ (ranging from $0.1$ to $1.6$ kpc) and $\sigma^{2}$, in order to
ensure solutions were not confined to local minima, and a uniform prior was
enforced on $\phi$ to keep it bounded between $50$pc and $5$kpc. The lower
limit was chosen as it is similar to the size of a single pixel at native
resolution for this dataset, and it is not possible to observe fluctuations
that are smaller than a single pixel (Cressie, 1993). For initial guesses of
$Z_{\text{char}}$ and $\nabla Z$ we use the results of a simple weighted-least
squares (WLS) fit, with an uncertainty of $0.05$ dex on $Z_{\text{char}}$ and
an uncertainty of $0.01$ dex/kpc on $\nabla Z$, as central metallicities and
metallicity gradients recovered for galaxies using the WLS method have been
found to agree with those found with a geostatistical method to these levels
(M22). After a test run, the chains of each Monte Carlo run were examined for
convergence. We found that for some runs, all chains converged to low values
of $\phi$ except for those that were started at values of $\phi$ that were
much larger than the best-fit value. In these cases, we re-ran an MCMC fit
using a set of initial values for $\phi$ ranging from $0.1$ kpc to $0.45$ kpc.
In a few cases, the upper limit on the prior on $\phi$ was also lowered to $3$
kpc, to ensure all chains remained close to the best-fit value, as informed by
previous runs. We note that adjusting parameter ranges based on information
inferred from previous fits is a standard practice in a Bayesian framework,
where the aim is to achieve quicker convergence on the global minimum while
efficiently using the (limited) computational resources available.

Figure 3: Corner plot, showing 1D and 2D marginalised posterior distributions
for the four parameters of our geostatistical model, using all Hii spaxels
from our fiducial galaxy (NGC 5236), with our fiducial metallicity diagnostic
(RS32), at the native resolution ($39$ parsec per spaxel). A strong
correlation can be seen between $\phi$ and $\log(\sigma^{2})$, the two
parameters that describe small-scale metallicity fluctuations within the
galaxy, as raising $\phi$ will naturally lower the variability of metallicity
within the data. Such a correlation is seen in all emcee fits performed in
this work.

To illustrate the quality of our fits, in Figure 3 we show the approximate
posterior generated by emcee as a corner plot for NGC 5236, using all Hii
spaxels, our fiducial metallicity diagnostic of RS32, and a datacube without
any additional binning. From this plot, we see that the two small-scale
parameters, $\phi$ and $\log(\sigma^{2})$, are strongly correlated with each
other ($\rho=0.94$). As $\phi$ is increased, nearby data points are expected
to become more and more correlated with each other, reducing the overall
variability of the model. To match the amount of variation that is seen in the
data set that is being fit, $\sigma^{2}$ must be increased to compensate for
this effect. Therefore, models that fit the data and have a high value of
$\phi$ must have a lower value of $\sigma^{2}$, and vice versa. Such behaviour
is expected, and does not impact on the results presented in the following
Sections.

##  4 Results

Figure 4: Recovered semivariograms for NGC 5236, when only a sub-sample of
spaxels are used. In all cases, the semivariogram has the same shape –
however, when the number of data points is $\lesssim 100$, small number
statistics lead to unreliable estimates of the variance between pairs of
spaxels at each separation, increasing the noise in the semivariogram, making
its interpretation less clear.

In this Section, we assess how the performance of the geostatistical methods
explored for galactic metallicity modelling in M21 and M22 depends on three
facets of the datacube used: the number of data points available, the
resolution of the data, and whether the fit is performed on individual Hii
spaxels, or over integrated Hii regions. However, separating out these effects
is difficult. When spaxels are binned to make coarser datacubes, the number of
data points available for statistical inference and model fitting decreases
(see Figure 2). Furthermore, integrating Hii spaxels into regions both
decreases the number of available data points for model fitting, and the
resolution of the available data points.

With this in mind, in Section 4.1, we first assess how these methods are
affected the number of data points with measured metallicity values available.
We test this by randomly selecting subsamples of a limited number of Hii
spaxels from our flagship galaxy NGC 5236 at native resolution, allowing us to
change the number of available data points for analysis without introducing
any additional spatial effects.

Then, in Section 4.2, we report on how the model parameters recovered and the
accuracy of metallicity predictions depend on the spatial resolution of the
data by comparing results for our binned datacubes for NGC 5236 to the results
found at native resolution. By comparing the results for the rebinned
datacubes to the results on the native resolution data sets of Section 4.1, we
can tell if there are any additional effects caused by having coarsely
resolved data that would not be expected by merely changing the number of data
points available.

Finally, in Section 4.3, we see if the results are the same for Hii regions as
they are for Hii dominated spaxels, by comparing spaxels that pass the BPT-
diagnostic cuts of Kewley et al. (2001) and Kauffmann et al. (2003) to regions
found using pyHIIextractor that also pass these cuts (for further details, see
Section 3.1). We interpret the results of this experiment in light of what we
have learned from Sections 4.1 and 4.2, to see if binning the data into Hii
regions has any impact on the output of the model that cannot be explained by
the other two effects investigated in this work.

###  4.1 The effects of missing data

To test how geostatistical methods perform when the number of available Hii
spaxels is small, we generate random subsamples of 50, 100, 300, and 500 Hii
spaxels from NGC 5236. We compare the performance of our geostatistical
methods on these subsampled datacubes to the full datacube, which has 7454 Hii
spaxels for which metallicities could be determined using RS32.

Figure 5: Best-fitting values, with uncertainty, of the metallicity
correlation scale $\phi$, the (log) total variance $\sigma^{2}$, the
characteristic metallicity $Z_{\text{char}}$, and the metallicity gradient
$\nabla Z$ for NGC 5236 as a function of the number of spaxels at native
resolution used for the fit. Generally, we find that the best-fitting
parameter values are consistent with the values recovered when all spaxels are
used (far right data point), albeit with large scatter when the number of
available data points is low.

We create a semivariogram using each of the subsamples of the native
resolution datacube in Figure 4, and plot them against the original
semivariogram computed for this data. Each semivariogram shows how the
variance in the metallicity between spaxels increases as the distance between
spaxels is increased, after accounting for the large-scale metallicity
gradient.111We provide a formal mathematical definition of the semivariogram,
an introduction to this method, and its relation to some other tools used in
astronomy to extract the covariance structure of random fields in Appendix
B.1. A bin size of $40$pc, chosen to match the resolution of the data and the
minimum separation between spaxels, is used for all semivariograms.

We see that all semivariograms have the same general shape, indicating that
the same small-scale structure is being detected in all five datasets.
However, we also see that as the number of available data points decreases,
the noise in the recovered semivariogram increases. When only $50$ data points
are used, some separations have missing values for the semivariance, as there
are not enough pairs of Hii spaxels separated by approximately that distance
for the variance between them to be estimated robustly. This issue could be
partially overcome by creating a semivariogram with a coarser binning,
allowing a more clear view of larger features, but doing so would result in
less retention of information about the small-scale variability of the data.
Since we are primarily investigating how many data points are necessary to
recover these small-scale properties, we do not do this.

In Figure 5, we compute estimates for all four parameters of our
geostatistical model for each subsample of the native resolution data set, and
investigate the convergence of these parameters. We see that even with as few
as 50 data points, the properties of the large-scale variability ($\nabla Z$
and $Z_{\text{char}}$) are recovered accurately. In fact, $\nabla Z$ is
consistent between all models, and $Z_{\text{char}}$ varies by only $0.02$
dex. Similarly, all models converge to the same value for the amount of
scatter around the mean trend ($\sigma^{2}$), although the uncertainty on this
parameter increases as the number of available data points decreases.

Recovered values on the metallicity correlation scale $\phi$ are consistent
when $300$ or more data points are available for fitting. When only $100$ data
points are used, the best-fit value of this parameter is larger, with larger
uncertainties ($\phi=518^{+271}_{-166}$ pc), but it is still consistent with
the value found at native resolution ($\phi=333^{+34}_{-28}$ pc) within the
1$\sigma$ error bars. For the model trained on only $50$ spatially resolved
metallicity measurements, a smaller estimate of $\phi$ is obtained, with
smaller error. By examining the posterior distribution of the Monte Carlo fit
for this trial, we find that this is due to estimates clustering tightly to
the lower limit of $\phi$ at $50$ pc, indicating that no evidence for spatial
correlations is seen in this limited data set.

Next, we evaluate how the goodness-of-fit of our geostatistical model is
affected by having a limited number of data points. We judge the quality of
our model fits over each subsampled datacube in terms of their predictive
power. Following M22, we assess this using ten-fold cross validation. This
method is a standard technique to estimate the ability of a model to
generalise to unknown data points (Field et al. 2012, Chapter 7.7.2.), and is
recommended as a reliable general method for model selection (Hogg & Villar,
2021). For each subsampled datacube, the Hii spaxels are split into ten evenly
sized groups. Data from all but one group are used to fit parameters for the
geostatistical model described in Section 3.3, and generate point predictions
for the metallicities of the Hii spaxels in the final group via the process of
universal kriging, outlined in Appendix B.2. This is repeated for all ten
groups to produce a predicted metallicity for every spaxel.

We summarise the results of this exercise by looking at the median absolute
deviation (MAD) of the predicted metallicity of each spaxel to the metallicity
observed for that spaxel, which we use as a metric for the model’s
performance. We then normalise this value to the variance of the data by
dividing by one standard deviation, $\sigma$, to make this metric easier to
interpret.222For reference, data from a Gaussian noise distribution around a
known mean is expected to have MAD$/\sigma=\sqrt{\frac{2}{\pi}}\approx 0.8$.

As a comparison point, we also repeat this cross-validation process for a
metallicity gradient model with no small-scale component, fit for each set of
data points using the standard weighted least squares (WLS) method,
implemented in the Python package statsmodels. In the language of our model,
this is equivalent to fixing $\eta(\vec{x})=0$.

We plot MAD/$\sigma$ for the linear gradient model and the geostatistical
model as a function of the number of data points available in Figure 6. For
the linear model with no small-scale component, we do not see a trend with the
prediction accuracy model changing with the number of data points over the
range we have chosen. This implies that $50$ data points is sufficient to
construct a metallicity gradient model; and using more data points in the fit
will not make this model any better. On the other hand, we see the accuracy of
the geostatistical model scales with the number of data points available for
training. With only $50$ spatially-resolved metallicities, the metallicity
predictions of a geostatistical model have comparable accuracy to the
predictions from a metallicity-gradient model. With $100-300$ data points, the
geostatistical model begins to outperform the azimuthally-symmetric model that
only captures large-scale metallicity variations within a galaxy. With a
larger sample of data points, the improvement increases dramatically.

This result is expected from geostatistical theory. Kriging works by taking
advantage of the fact that nearby data points are expected to have similar
properties (Tobler 1970, Appendix B.2). When more data points are used, it is
more likely that each data point lies within $\lesssim\phi$ parsecs of another
spaxel with known metallicity, allowing the metallicity to be inferred using
local galaxy properties as well as global ones. Furthermore, when more data
points are used, the recovered values of the small-scale structure parameters
$\phi$ and $\sigma^{2}$ are more accurate.

Figure 6: The normalised prediction error (median absolute deviation per
standard deviation) of the predictions for metallicity computed using the
cross-validation method for NGC 5236 using a simple linear gradient (blue
line) and a geostatistical model (red line). As long as the number of data
points $\gtrsim 100$, predictions are more accurate when a geostatistical
model including small-scale metallicity variations was included.

###  4.2 The effects of spatial resolution

To determine the resolution limit at which small-scale metallicity variations
are visible, we construct a semivariogram for the small-scale metallicity
variation using data convolved to each binning factor $f$. We plot these
results in Figure 7.

Figure 7: The semivariogram for the TYPHOON metallicity map of NGC 5236,
binned to eight different spatial resolutions. As the resolution becomes
coarser, the total variance in metallicity between spaxels decreases, as small
scale inhomogeneities are binned together. When the spatial resolution is
finer than $\sim 250$pc ($f\leq 6$), the shape of the semivariogram matches
that of the native resolution map, indicating that the resolution is fine
enough for small-scale metallicity fluctuations to be accurately resolved. For
resolutions of $300-400$ pc per pixel, the general shape of the semivariogram
can still be seen, and an estimate of the scale at which metallicities become
locally uncorrelated can still be made. At resolutions coarser than this, it
is difficult to infer the presence of small-scale fluctuations, and accurately
describe their properties.

The most striking trend that can be seen in this series of plots is that the
height of the semivariograms at large separations decreases as $f$ is
increased. This implies that the total amount of variance about a linear
gradient is lower when the data is of a coarser quality. Qualitatively, this
is to be expected. When the spatial resolution of IFU spectroscopy is coarser
than the size of metallicity variations, different features are averaged
together over each individual spaxel, decreasing the variance around a galaxy-
scale metallicity model. Quantitatively, we explore how this spatial averaging
would impact the variance seen at different resolutions if the true small-
scale metal variations follow an exponential correlation function with a
correlation scale of $\phi=330$ pc, as is suggested by our highest-resolution
data, in Appendix C. We find that the reduction of variance seen at different
resolutions matches the reduction predicted with this model of spatially-
correlated metallicity variations very tightly.

The separation at which the semivariance stops increasing (called the range of
the data in geostatistics literature; e.g. Cressie 1990) is not seen to
increase with $f$ for all semivariograms computed with $f\leq 6$
(corresponding to physical resolutions finer than 250 pc per pixel, plotted as
blue lines in Figure 7). Visually, when $f\leq 6$, all semivariograms show the
same shape, beginning to curve downwards at a separation of $\sim 0.4$ kpc,
before flattening out completely by $\sim 1$kpc. For coarser data ($f>6$,
plotted as brown lines in Figure 7), this shape is no longer clearly visible.
At $f=8$, flattening is still seen at $\sim 1$kpc, but the shape of the
semivariogram for sub-kpc data is no longer accurate. At $f>8$, we see the
apparent range of the data increase with $f$, as trends associated with
physical small-scale variations become washed out, and the size of the pixels
becomes the most important scale for setting the observed structure of the
metallicity map.

To show this result in a more quantitative way, we present the best-fit value
of $\phi$, the metallicity correlation scale, with uncertainty, as a function
of the resolution of available data in Figure 8. Using the native resolution
data, the most likely value of $\phi$ is $332^{+33}_{-28}$ pc. As $f$ is
increased from $1$ to $10$, the best-fit value of $\phi$ increases slightly
and has slightly larger uncertainties, but remains consistent with the value
found for native resolution data. At $f=12$, the size of a single spaxel
becomes $470$ pc – larger than the value of $\phi$ recovered for the high-
resolution maps. At this resolution, our model fitting procedure returns a
value of $\phi=527^{+135}_{-90}$ pc, inconsistent with the value of $\phi$
suggested by the higher-resolution maps. At our coarsest resolution, $\phi$ is
inferred to be $443^{+308}_{-244}$ pc. While this value is consistent with the
$\phi$ recovered at the finest resolutions, the uncertainties are very large,
reflecting the fact that the value of $\phi$ cannot be well-constrained by
this data, as a single pixel is expected to be much larger than the typical
size of small-scale variations.

Figure 8: Best-fitting values, with uncertainty, of the metallicity
correlation scale $\phi$, the characteristic metallicity $Z_{\text{char}}$,
and the metallicity gradient $\nabla Z$ for NGC 5236 as a function of the
spatial resolution of metallicity maps, reported in parsec (bottom axis) and
effective radii ($R_{e}$, top axis). In the top panel, the grey shaded region
shows the area of parameter space where the correlation scale $\phi$ is
smaller than the size of a pixel. We find that kpc-scale resolved metallicity
surveys are able to accurately recover large-scale trends of the metallicity
profile ($Z_{c},\nabla Z$), while the correlation scale can only be accurately
constrained by observations with a resolution $\lesssim 400$pc, or
$0.12R_{e}$.

Figure 8 also shows that the recovered values of $Z_{\text{char}}$ and $\nabla
Z$ are relatively consistent throughout all data maps explored, excepting the
results found for the datacube at native resolution, where a flatter
metallicity gradient is preferred. We note that the fitted parameters for this
model at native resolution are the same as those shown in Figure 5. On its
own, this result is not significant – while the values of $Z_{\text{char}}$
and $\nabla Z$ recovered using our emcee-based parameter fitting pipeline do
not agree with any other recovered values of $Z_{\text{char}}$ and $\nabla Z$
at other resolutions to $1\sigma$, we would not expect the $1\sigma$ credible
intervals of eight different model fits to all overlap. However, we also find
a similar result of the recovered value of $\nabla Z$ being flatter when two
out of our three alternative metallicity diagnostics are used (N2S2H$\alpha$
and Scal, with results inconclusive for O3N2; see Appendix A), and for the
galaxy NGC 7793 when the size of a spaxel is $\lesssim 100$ pc (see Section
5). While formally significant, the size of this deviation is small – smaller
than the difference between metallicity gradients recovered for this galaxy
when different metallicity diagnostics are used (see Figure 9 of Poetrodjojo
et al. 2019).

This result is unexpected from previous studies of the effects of resolution
on recovering the metallicity gradient. In fact, previous research has shown
that a flattening of the metallicity gradient is expected for poorly
spatially-resolved data due to the effects of pixelisation and beam smearing
(e.g. Yuan et al., 2013; Acharyya et al., 2020). However, in this case, it is
not a degradation of the data that causes an apparent flattening of the
observed gradient, as a tendency for $Z$ to decrease as $r$ increases is still
seen even in the native resolution datacube.

To test the validity of our statistical pipeline, we analysed the values of
$\nabla Z$ recovered using a simple WLS approach without incorporating any
small-scale variation. We present our results in Figure 9. At each resolution,
we plot the measured metallicity of each Hii spaxel, with error bars, at its
distance from the galaxy centre. We overplot the best-fit linear relations
describing the large-scale metallicity trends recovered using both our
geostatistical approach, and using least-squares. For all of the binned data
sets ($f>1$), the profiles recovered using WLS and our geostatistical method
show good agreement. However, when $f=1$, the geostatistical approach prefers
to describe the data with a flatter metallicity gradient. This difference is
significant, with $p<5\times 10^{-5}$ for the marginalised posterior value of
$\nabla Z$ being consistent with the least-squares value. Visually, at $f=1$,
both lines appear to describe the data fairly well. While there is slightly
more variance in the data points around the flat profile predicted by the
geostatistical model, this is not unexpected given the different ways that
these two models were optimised – in the WLS approach, parameters for the line
of best fit are chosen to minimise the scatter of all data points around the
radial metallicity profile, whereas in our geostatistical model, data points
are naturally expected to have some scatter around the large scale trend (for
$f=1$, the most likely value for this scatter is $\sigma=0.05$ dex).

Figure 9: The radial metallicity profile of NGC 5236, computed using (i) our
geostatistical methodology (red dashed lines) and (ii) the weighted least
squares (WLS) algorithm (orange dashed lines) at different rebinned
resolutions. When the WLS method is used, consistent values of the metallicity
gradient are recovered for all binning factors $f$, with a slight tendency for
the metallicity gradient to become flatter at the coarsest resolutions. For
most of the resolutions tested, the radial metallicity profile recovered using
our geostatistical approach agrees with the profile recovered using WLS.
However, for the unbinned data ($f=1$), the geostatistical method prefers a
flatter metallicity gradient.

When the WLS method was used, the metallicity gradient recovered for the
highest resolution map ($\nabla Z=-0.0126$ dex/kpc) was within the range of
metallicity gradients recovered for the other rebinned datacubes
($-0.013\leq\nabla Z\leq-0.010$ dex/kpc), with a slight trend of the
metallicity gradient becoming flatter at coarser resolutions, as is expected
from a coarser spatial sampling. However, for the N2S2H$\alpha$ and Scal
diagnostics, an increase in the WLS-reported metallicity gradient was seen at
the finest resolutions. We discuss a possible interpretation for this effect
in Section 4.2.1.

In Figure 10 we plot the normalised prediction error (MAD$/\sigma$, as
calculated using ten-fold cross-validation) of the geostatistical model, and
the linear model with no small-scale component, against the resolution of the
datacubes. We find that when a linear gradient model is used, a similar
prediction error value of MAD$/\sigma=0.75-0.85$ is seen at all resolutions,
with a slight trend of the predictive accuracy decreasing as the resolution is
made coarser. At all resolutions investigated, predictions from a
geostatistical model were found to be more accurate than those from a
metallicity-gradient model (although this improvement is only marginal for the
coarsest resolution datacube with $f=24$). For resolutions finer than $\sim
330$ pc (the scale at which adjacent pixels are separated by a distance less
than the metallicity correlation scale, $\phi$), the prediction accuracy of
the geostatistical method improves rapidly as the resolution is improved. This
highlights the power of this geostatistical modelling framework for being able
to produce accurate predictive maps of the 2D metallicity distributions of
galaxies, especially for high-resolution datacubes.

We note that while the geostatistical model fit to the highest resolution data
has a value of $\nabla Z$ that does not agree with the other models to
$1\sigma$, this model is still able to predict the metallicity of Hii regions
with more accuracy than any other model. This indicates that the predictive
accuracy of our geostatistical model is not severely impacted by the different
values of $\nabla Z$ recovered – the model with a lower central metallicity
and a flatter metallicity gradient is just as suitable for describing the
observed data as a model with a stronger negative metallicity gradient.

Figure 10: Normalised prediction error (median absolute deviation per standard
deviation) for metallicity computed using the cross-validation method for NGC
5236 using a simple linear gradient (blue line) and a geostatistical model
(red line). For all resolutions tested, predictions were more accurate when a
geostatistical model including small-scale metallicity variations was
included, with the greatest improvement seen at fine resolutions.

####  4.2.1 Impact of spatial resolution on Hii spaxel selection

There are many possible reasons why the metallicity gradients recovered with
the finest resolution spaxels exhibit small, but statistically significant
deviations compared to metallicity gradients measured at simulated coarser
spatial resolutions.

First, it could be that the intrinsic radial trend in metallicity for these
galaxies is not linear. High resolution metallicity maps have revealed breaks
in a single metallicity profile that may not be visible in coarsely resolved
data, potentially due to the influence of galactic bars (e.g Chen et al.,
2023). Flattening in metallicity profiles have been reported in the extended
discs of spiral galaxies beyond $R_{25}$ (e.g. Rosales-Ortega et al., 2011;
Goddard et al., 2011), including for this galaxy (Bresolin et al., 2009).
However, the data used in this study is limited to a maximum galactocentric
distance of $\sim R_{25}$, where flattening has not been found in previous
studies of this galaxy’s metallicity profile. Other studies have reported
flatter metallicity profiles in the inner regions of galaxies (e.g. Belley &
Roy, 1992; Zinchenko et al., 2016), but we do not see this in NGC 5236 at any
resolution (see Figure 9).

The differences in inferred gradients could also be related to the selection
of spaxels at each resolution. While we use a consistent Hii spaxel selection
method for all of our rebinned datacubes, the set of spaxels that are
identified as being Hii-dominated is not the same between the emission line
maps with different resolution. This effect is visible in Figure 9 – at $f=2$,
several spaxels are classified as Hii-dominated at radii of $\sim 11$ kpc that
have no analogues in the $f=1$ map.

Additionally, (1) the Hii region selection includes an empirically calibrated
threshold from Kauffmann et al. (2003) based on Sloan Digital Sky Survey data
that may not apply to apply to highly-resolved datacubes with resolutions
finer than 100pc, and (2) the BPT diagnostic lines of Kewley et al. (2001) are
computed to be the upper limits of the locations on the BPT diagrams that can
be inhabited by Hii regions – however, no guarantees are made that spaxels
that are not ionised by Hii emission must lie outside of this region.
Therefore, if a statistically appreciable number of spaxels outside Hii
regions are instead flagged by the DIG diagnostics as Hii dominated for
finely-resolved data, this could lead to a bias in the metallicity gradient
computed for high-resolution data.

To explore this hypothesis, we re-ran the analysis of Section 4.2, ensuring
that the same selection of spaxels were used at each resolution. To do this,
we followed a procedure similar to the one outlined in Poetrodjojo et al.
(2019). First, using the native resolution data, for each metallicity
diagnostic we define a mask containing only the spaxels that passed all data
quality control cuts. When binning the data to coarser resolutions, only the
spaxels that passed all data quality cuts at the finest resolutions were
aggregated. For this test, we did not recompute line ratios for each rebinned
spaxel with LZIFU. Instead, we simply summed the emission line fluxes of each
spaxel within the mask contained within each binned spaxel, summing the errors
in quadrature.

Similar to Poetrodjojo et al. (2019), and following M22, a harsher quality
control cut for defining Hii spaxels was used than the BPT diagram based
methods. Based on the methods of Kaplan et al. (2016), we compute the fraction
of light in each spaxel that is coming from Hii regions ($C_{\text{H{ii} }}$)
using the line ratio [S ii]/H$\alpha$. In more detail, using the native
resolution maps, we assumed that the 100 spaxels with the brightest H$\alpha$
flux were entirely ionised by Hii regions ($C_{\text{H{ii} }}=1$), and the 100
spaxels with the faintest H$\alpha$ flux (that still have a S/N ratio greater
than $3$ for both [S ii] and H$\alpha$) were entirely ionised by DIG
($C_{\text{H{ii} }}=0$). We take the median values of [S ii]/H$\alpha$ of each
of these populations to be the typical values of [S ii]/H$\alpha$ for Hii
regions and the DIG, respectively, and computed values of $C_{\text{H{ii} }}$
for all spaxels by linearly interpolating where the measured value of [S
ii]/H$\alpha$ falls between these two extremes. Spaxels with a value of
$C_{\text{H{ii} }}>0.95$ are considered to be Hii dominated and are included
in our mask.333This threshold value of $0.95$ is (i) intrinsically arbitrary,
and (ii) stricter than the threshold value of $0.9$ used in Poetrodjojo et al.
(2019). Such a conservative choice is consistent with the analysis presented
in M22, where a 0.9 threshold was applied to $\log($[S ii]/H$\alpha)$, which
for the TYPHOON emission line map of NGC 5236 is equivalent to a 0.95
threshold in [S ii]/H$\alpha$.  Additionally, any spaxels with S/N < 3 in
H$\alpha$, H$\beta$, or any emission lines used in the metallicity diagnostic
were discarded.

Figure 11 shows how the number of spaxels available for our flagship galaxy
(NGC 5236) with our fiducial metallicity diagnostic (RS32) changes as a
function of resolution with this new Hii classification scheme. We compare
this with the number of spaxels retained using the BPT diagrams and surface
brightness cut as well. We see that, for the native resolution map, the [S
ii]/H$\alpha$-based cut is much stricter than the BPT diagram-based cut: while
$7488$ spaxels are classified as Hii dominated at native resolution using the
BPT diagram-based diagnostics, only $2351$ pass the [S ii]/H$\alpha$-based cut
and have S/N > 3 for H$\alpha$, H$\beta$, [S ii]$\lambda\lambda 6717,31$, and
[O iii]$\lambda 5007$. As the resolution is made coarser, the number of
spaxels available for fitting decreases much less rapidly when the same set of
spaxels are used at all resolutions. This may be due to an excess of spaxels
being included at fine resolutions when the BPT-based spaxel selection
criteria are used; or it could be due to spaxels being lost to DIG
contamination at coarse resolutions.

Figure 11: Number of data points available for NGC 5236 as a function of
resolution, when Hii spaxels are selected (i) using the BPT diagnostics of
Kewley et al. (2001) and Kauffmann et al. (2003) together with the H$\alpha$
surface brightness cut of (Zhang et al., 2017); and (ii) by selecting a sample
of Hii spaxels at native resolution using the line ratio [S ii]/H$\alpha$, and
rebinning this set of spaxels only to coarser resolutions. At the finest
resolutions, the [S ii]/H$\alpha$-based cut removes many more spaxels than the
BPT diagram based cuts. At coarser resolutions, however, more spaxels are
retained, as binned spaxels that contain even one spaxel at native resolution
that passes the data control cuts at native resolution are kept.

Using a consistent selection of spaxels, but varying the resolution of the
data, we then repeat our experiment of Section 4.2, exploring how our set of
recovered parameters vary with resolution. We plot our results in Figure 12.
We find that in this case, the metallicity gradient varies much less than was
found for the case where Hii spaxels were selected using a BPT-based approach.

Figure 12: Parameters of our geostatistical model recovered with emcee after
ensuring that the same set of spaxels are used at all resolutions. For this
data, the value of $\nabla Z$ recovered does not become significantly flatter
as the resolution is changed.

A similar result was seen in the investigations of Poetrodjojo et al. (2019).
When DIG was excluded from the data at native resolution, metallicity
gradients were found to be mostly constant with resolution. However, large
variations in the metallicity gradient were found when different metallicity
diagnostics were used, suggesting that systematic uncertainties in diagnostics
are the most important limiting factor to achieve precise gradient
measurements, rather than any effects resulting from having data with very
fine spatial resolution.

###  4.3 The effects of aperture bias

In Sections 4.1 and 4.2, we assumed that all metallicities measured for Hii
spaxels are unbiased tracers of the underlying metallicities of the regions of
the galaxy that they are drawn from. However, there are reasons to consider
this may not be the case. Spherically-symmetric models of the interior
structure of Hii regions such as Cloudy (Ferland et al., 2013) indicate that
the ratio of emission lines within a Hii region is expected to change as a
function of the distance from the centre (see e.g. Figure 1 of Mannucci et al.
2021). For this reason, when the size of an IFS spaxel is smaller than the
size of a Hii region, the relative strengths of its emission lines are
dependent not only on the metallicity within the spaxel, but also on the
distance from the Hii region’s centre to the spaxel’s centre, and on the
relative sizes between the two. This may introduce an additional source of
error in the metallicity estimated within each spaxel. By comparing the
metallicities inferred from integrated DIG-free spectra to those that would be
inferred from narrow slit spectroscopy, Mannucci et al. (2021) showed that in
some circumstances aperture bias may have a stronger effect on recovered
metallicities than DIG contamination. However, other studies with different
methodologies have found the effects of aperture bias to be mild (e.g.
Arellano-Córdova et al., 2022).

Figure 13: Comparing semivariograms calculated using Hii dominated spaxels vs
integrated Hii regions found with pyHIIextractor. We see that individual
spaxels carry more noise than Hii regions, and show a small but noticeable
increase in variance from $\sim 0.5-1$ kpc that is not seen in the region-
based semivariogram. The additional fluctuations in the region-based
semivariogram are likely caused by a low number of data points compared to the
spaxel-based semivariogram.

One step that can be added to an IFS data processing pipeline as a
prophylactic measure against aperture bias is to first identify Hii regions
within a galaxy, and then integrate the emission from each spaxel within each
region to get emission line ratios for each Hii region. The line ratios for
each region may then be converted into metallicities for each region. This
ensures that no line ratio is artificially increased or decreased due to
spatial sampling effects.

However, this method comes at some costs. Firstly, by binning multiple
adjacent spaxels into single regions, the spatial resolution of the data is
decreased, and the number of available data points that can be used to fit a
2D metallicity model is lowered. Secondly, many Hii region-finding algorithms
make assumptions on the size, brightness, and geometry of Hii regions that may
not be valid (see e.g. Jin et al. 2022), yield different results at different
spatial resolutions (Mast et al., 2014), and usually depend on parameters that
must be hand-tuned by the user, making the results of a Hii region search
difficult to replicate (Lugo-Aranda et al., 2022).

In this context, a key question to ask is whether it is better to bin data
into Hii regions before beginning a geostatistical analysis, or to use higher-
resolution, possibly biased data from individual Hii spaxels.

In Figure 13, we show semivariograms for the small scale metallicity
variations within our fiducial galaxy (NGC 5236) using our fiducial
metallicity diagnostic (RS32), comparing results found when all spaxels are
used (blue line) to results using data binned to 485 Hii regions (red line),
found using pyHIIextractor that pass both BPT diagnostics (Section 3.1). We
see that both semivariograms show the same shape, with the same flattening at
$\sim 0.4$ kpc. Two other effects are also visible. First, the semivariogram
produced using Hii regions is more noisy. This is expected, because there are
fewer data points available for the fitting – see Section 4.1 and Figure 4.
Secondly, the height of the semivariogram in the large-scale limit computed
using Hii regions is shorter than the semivariogram computed using Hii
spaxels. This, too, is expected, because the size of a Hii region is larger
than a Hii spaxel, and so averaging over a larger region reduces the amount of
variance in metallicity (see Figure 7, and further discussion in Appendix C).

By comparing the height of the two semivariograms near zero separation, we see
that the Hii regions have far less uncorrelated variance than what is seen for
Hii spaxels. This observation has two possible interpretations. The first is
that integrating over Hii regions washes out small-scale metallicity
variations that exist within an individual Hii region. The second
interpretation (which we consider more likely) is that metallicities estimated
for integrated Hii regions are more accurate than the metallicities estimated
for individual Hii spaxels. If this interpretation is correct, then by
correcting for aperture bias, the error on the metallicity measurement of each
spaxel is reduced. By comparing the height of these semivariograms near zero
separation, the amount of variation reduction when this formulation is used
could be estimated, potentially allowing the effects of aperture bias to be
quantified for a range of commonly-used metallicity diagnostics. We postpone
such an analysis for a future study.

Overall, integrating spaxels into Hii regions yields a metallicity map with
fewer data points, coarser resolution, and more accurate metallicity
measurements. In Table 2, we show what effect this integration step has on the
best-fit values of $Z_{c},\nabla Z,\phi,$ and $\sigma^{2}$ recovered for this
galaxy. We find that most of these parameters agree within their $68\%$
credible intervals. The variance is found to be significantly lower for Hii
regions, which is an understandable consequence of the coarser resolution of
the binned Hii region data compared to the unbinned Hii spaxel data (see
Appendix C). A $3.6\sigma$ difference is seen between the best-fit values of
$Z_{\text{char}}$ for binned regions and individual spaxels, possibly pointing
to a systematic error in metallicity measurements when this diagnostic is used
on regions smaller than an individual Hii region. However, the actual
difference is very small, with the median values computed for regions and
spaxels differing by only $0.02$ dex. This implies that the effect of such a
systematic bias, if it exists, is small enough to be ignored, particularly
given the larger (0.7 dex) uncertainties associated with strong emission line
ratio based metallicity diagnostics (Kewley & Ellison, 2008).

For the other two parameters, there are some slight differences within the
bounds of uncertainty. For the Hii region data, $\nabla Z=-3.7\pm 2.7\times
10^{-3}$ dex kpc-1. This is slightly more negative than the gradient found
using native resolution Hii spaxel data, but flatter than the gradients
recovered for binned data, indicating that this data may be close to the limit
of resolution at which the preferred models become those with no (or very
weak) metallicity gradients, and variations within the galaxy is
preferentially described through small-scale trends. Interestingly, the value
of $\phi$ found for the Hii region based data is lower than the value
recovered when Hii spaxels are used – but the size of this effect is not
strong enough to be significant.

| Hii spaxel fit | Hii region fit
---|---|---
$Z_{\text{char}}$ | $8.76\pm 0.004$ | $8.78\pm 0.004$
$\nabla Z$ ( 10-3 dex kpc-1) | $-0.2^{+2.5}_{-2.4}$ | $-3.7\pm 2.7$
$\phi$ (pc) | $332^{+33}_{-28}$ | $275^{+38}_{-32}$
$\log_{10}\left(\sigma^{2}\right)$ | $-2.61\pm 0.04$ | $-2.50\pm 0.04$
Table 2: Best fit parameter results for the metallicity model for NGC 5236
when results are computed using data from each spaxel associated with a Hii
region, compared to results when using integrated emission line ratios from
Hii regions found using pyHIIextractor. In both cases, all spaxels/regions lie
below the Kewley et al. (2001) and Kauffmann et al. (2003) demarcation lines
in the S2- and N2-BPT diagrams. We find significant differences in the
recovered values of two parameters for this galaxy with this diagnostic: the
variance, which is understandable as Hii regions are binned over a larger area
than individual Hii spaxels, naturally decreasing the variance in the
metallicity between Hii regions (see Appendix C), and the characteristic
metallicity, $Z_{\text{char}}$, possibly reflecting a small systematic bias
when this diagnostic is used on spaxels smaller than an individual Hii region.

Based on this analysis, we conclude that binning Hii spaxels into Hii regions
does reduce the uncertainty of the metallicity estimated at each data point,
and does not significantly impact the results of a geostatistical model fit.
As long as the data remains at a resolution finer than the expected size of
small-scale variations within a galaxy, and the number of data points remains
high ($\gtrsim 100$), we advise binning Hii spaxels into Hii regions before
fitting a geostatistical model. However, if this is not possible, there is no
cause for concern, as a model fit using unbinned Hii spaxels will likely yield
similar results.

A similar conclusion was reached by Vogt et al. (2017). They analysed the two-
dimensional abundance distribution of HCG 91c captured by MUSE in two ways:
(i) spaxel-by spaxel, and (ii) binning the data into Hii regions using the
open-source Python fitting code contained in
brutus444https://fpavogt.github.io/brutus/. They found that the metallicity
gradient extracted was consistent when either spaxels or Hii regions were
used. Additionally, Vogt et al. (2017) noticed coherent, spatially-localised
variations in metallicity on sub-kpc scales when both methods were used,
representing an amplitude fluctuation of $0.2$ dex when computed from a
spaxel-based analysis, and $0.15$ dex when the data was binned into Hii
regions. Such fluctuations agree with what we find in our analysis, with the
spatially correlated variation remaining after spaxels were binned, albeit
with a slightly lower variance.

This result also shows that the fluctuations that we reveal with a
semivariogram analysis are not caused by the varying ionisation conditions
within Hii regions. Rather, they capture intrinsic variations in the
metallicity of the ISM between Hii regions.

##  5 Results for NGC 6744 and NGC 7793

In Section 4.2, we find that, based on the galaxy NGC 5236 and the RS32
diagnostic, the limiting resolution at which our geostatistical pipeline can
accurately recover the statistical properties of small scale metallicity
fluctuations is $\sim 350$ pc. This value is comparable to the best-fit value
of $\phi$, the scale of these fluctuations. Because of this, it is unclear as
to whether this resolution limit depends on the size of the intrinsic
fluctuations in NGC 5236, or is constant over a wider galaxy sample.

We aim to make a qualitative assessment of the robustness of this resolution
limit. For this, we expand our analysis to two other galaxies – NGC 6744, an
intermediate spiral galaxy at a distance of 11.6 Mpc, and NGC 7793, a
flocculent spiral galaxy 3.6 Mpc away. Each of these galaxies was binned to
increasingly coarser resolutions, using values of $f=1,2,4,6,8,10,$ and $12$.
This corresponds to physical resolutions ranging from $93$ to $1113$ pc
($0.012-0.14R_{e}$) for NGC 6744, and $28$ to $346$ pc ($0.016-0.19R_{e}$) for
NGC 7793. For NGC 7793, a binning with $f=24$ (corresponding to a physical
resolution of $691$ pc) was also attempted, but only 31 spaxels satisfied our
data cuts at this resolution, preventing a meaningful geostatistical model
from being fit (c.f. Section 4.1).

At each resolution, a geostatistical model was fit to each galaxy using emcee.
We show the median, 16th and 84th percentiles from our posterior distributions
of $\phi$, $Z_{\text{char}}$, and $\nabla Z$ for NGC 6744 at each resolution
in Figure 14, and for NGC 7793 in Figure 16, using RS32 as our metallicity
diagnostic. We also show the predictive accuracy of our geostatistical models
using ten-fold cross-validation for NGC 6744 in Figure 15 and for NGC 7793 in
Figure 17. We compare the metallicities predicted for each spaxel via
universal kriging to the measured metallicities of each spaxel, comparing it
against the predictive accuracy a simple linear gradient model fitted via the
weighted least squares (WLS) method.

For NGC 6744, at native resolution, the most likely value of $\phi$ is
$1.30^{+0.15}_{-0.12}$ kpc – much larger than the value of $\phi$ found for
NGC 5236, even when NGC 5236 was binned to a resolution coarser than NGC
6744’s resolution. This indicates that correlations between metallicities
exist on much larger scales in this galaxy than in NGC 5236. Here, we see that
the correlation scale $\phi$ is fairly stable to increasing values of $f$
until a resolution of $742$ pc per pixel is reached, comparable to half the
value of $\phi$ found at native resolution. After this point, the best fit
values of $\phi$ become significantly higher, and the size of their error bars
increases dramatically. This in turn affects the predictive accuracy of this
model, which we plot in Figure 15. After the point where the fit to $\phi$
stops being accurate, the geostatistical model becomes markedly worse, yet its
predictive performance remains similar to a metallicity gradient model with no
small-scale component.

For NGC 6744, the best fit values of the metallicity gradient and
characteristic metallicity are stable for all resolutions tested, with all
results consistent with each other within $1\sigma$ error bars.

Figure 14: Best-fitting values, with uncertainty, of the metallicity
correlation scale $\phi$, the characteristic metallicity $Z_{\text{char}}$,
and the metallicity gradient $\nabla Z$ for NGC 6744 as a function of the
spatial resolution of metallicity maps. The grey-shaded region in the top
panel represents the region where $\phi$ is less than or equal to the spaxel
size. Figure 15: Normalised prediction error (median absolute deviation per
standard deviation) of the predictions for metallicity computed using the
cross-validation method for NGC 6744 using a simple linear gradient (blue
line) and a geostatistical model (red line). When the small-scale structure
can be resolved accurately, predictions from a geostatistical model are more
accurate. However, at resolutions coarser than $\sim 750$ pc ($\sim
0.1R_{e}$), kriging predictions are no more accurate than predictions from a
linear model with no small scale component.

For NGC 7793, at native resolution, the most likely value of $\phi$ is
$254^{+24}_{-21}$ pc. Best fit values of $\phi$ are seen to deviate
significantly from this value with large error bars when data with a
resolution of $\gtrsim 170$ pc were used. This implies again that the result
was stable to increases in resolution until a pixel becomes larger than
$\sim\frac{1}{2}\phi$.

Figure 16: As in 14, but for NGC 7793. Figure 17: Normalised prediction error
(median absolute deviation per standard deviation) of the predictions for
metallicity computed using the cross-validation method for NGC 7793 using a
simple linear gradient (blue line) and a geostatistical model (red line). When
the small-scale structure can be resolved accurately, predictions from a
geostatistical model are more accurate. However, at resolutions coarser than
$\sim 200$ pc ($\sim 0.1R_{e}$), kriging predictions are no more accurate than
predictions from a linear model with no small scale component.

We show how the predictive accuracy of our model is affected by resolution for
NGC 7793 in Figure 17. At the two sharpest resolutions tested, the
geostatistical model is able to predict metallicities at locations within the
galaxy at a much greater accuracy than is possible with a simple linear
gradient. At these resolutions, individual regions of enhanced star formation
are well-resolved, and modelling small-scale deviations in metallicity
captures much of the details of metallicity variability within this galaxy. At
resolutions between $100$ and $200$ pc per spaxel, small-scale structure
within this galaxy can still be marginally resolved and identified by our
geostatistical pipeline, leading to models that, while not as successful at
identifying local deviations from a large-scale metallicity trend, still
exhibit a significant improvement in predictive power compared to a linear
gradient model. When this galaxy is observed with resolutions coarser than
$\sim 200$ pc per spaxel, the size of a pixel becomes comparable to the
typical size of a small-scale metallicity fluctuation, and any significant
benefit in predictive power from fitting a geostatistical model is lost.

Overall, this analysis shows that the resolution requirements for a
geostatistical analysis to be valuable does indeed depend on the spatial
extent ($\phi$) of the small-scale variations within a galaxy. This makes
coming up with general guidelines for the resolution requirements needed in
order to study the internal metallicity distribution of galaxies using these
geostatistical techniques difficult, especially when the typical size of
small-scale fluctuations is not known in advance. However, we have noticed
that the resolution at which geostatistical models begin to perform poorly
often coincides with the resolution at which small-scale structure within a
galaxy stops being apparent on a visual inspection.

This leads us to suggest the following “Golden Rule” for when a geostatistical
analysis is appropriate for IFS datacubes: when small-scale deviations from a
linear metallicity gradient appear to exist _on visual inspection_ , the size
(in terms of both $\phi$ and $\sigma^{2}$) of these small-scale fluctuations
can be captured accurately, and used to train predictive geostatistical models
which can improve our understanding of small-scale metal mixing processes
within galaxies. Our early results indicate that this is often achieved when
the size of a resolution element is $\lesssim 0.1R_{e}$. However, we caution
the reader that this guideline is only based on observations of three large
local star-forming spiral galaxies, and it is not yet known if this relation
will hold true for dwarf galaxies, or galaxies with other morphologies. At
resolutions where small-scale fluctuations appear to exist but are not well-
resolved (where $\phi$ is comparable to but smaller than the size of a pixel),
geostatistical models offer some improvement in predictive power over a simple
linear gradient model. Still, the best estimates of any parameters of the fit
should be taken with some healthy scepticism, as poor resolution has been
shown to bias $\phi$ to larger values, and lower $\sigma^{2}$ (see Appendix
C). Finally, when the data is so coarse that no structure aside from a radial
trend can be seen, there will not be enough data to fit a geostatistical model
that is more informative than a metallicity gradient would be on its own.

##  6 Discussion

From the results presented in Section 4.2 and Section 5, we find that a
resolution finer than $\sim\frac{1}{2}\phi$ is required to accurately recover
the statistical properties of small-scale metallicity variations. From the
literature (e.g. Kreckel et al. 2019, M21, M22), small-scale fluctuations in
ISM metallicity are seldom seen on scales larger than 1kpc (however, see
Williams et al. 2022)

For this reason, we do not consider it valuable to use this geostatistical
methodology to analyse data products of $\sim 1$kpc IFS surveys, such as SAMI
(Croom et al., 2012), MaNGA (Bundy et al., 2015), or CALIFA (Sánchez et al.,
2012). This is unfortunate, as these surveys contain a large number of
galaxies spanning multiple orders of magnitude in mass and star formation
rate, which would allow investigations into how the properties of small-scale
metallicity variations are set by a galaxy’s global properties. However, there
still exists a number of IFS surveys containing dozens of galaxies with
spatial resolutions finer than $\sim 250$ pc. The combined catalogue AMUSING++
(López-Cobá et al., 2020) contains IFS observations of 635 galaxies observed
with spaxel sizes of $0^{\prime\prime}.2$ and a typical seeing of
$1^{\prime\prime}.0$. Approximately $40\%$ of galaxies in this sample are
close enough to be resolved with a seeing finer than $\sim 300$ pc, making
this a statistically significant legacy dataset of IFS observations for which
these geostatistical techniques could be applied. This dataset also contains
galaxies of different morphologies. Analysing such a diverse population of
galaxies using a geostatistical pipeline could provide insights into how the
different drivers of turbulence affect gas mixing in a galaxy throughout its
lifetime.

Similarly, geostatistical analysis could be applied to very high resolution
cosmological simulations of galaxy formation. Specifically, applications would
be ideal for zoom simulations such as FIRE (Wetzel et al., 2023) and AURIGA
(Grand et al., 2017), which simulate metallicity variations within galaxies on
spatial scales as fine as $\sim 2$ pc. Such simulated datasets could readily
be analysed with our geostatistical pipeline, in order to test the sub-grid
physical prescriptions adopted in these simulations against observations, and
thus both validate/falsify current choices and inform future theoretical
models.

Extending a geostatistical analysis to the population of galaxies at cosmic
noon and beyond will be possible with the next generation of extremely large
thirty meter class ground-based telescopes. The European Extremely Large
Telescope (E-ELT; Gilmozzi & Spyromilio 2007) and the Thirty Meter Telescope
(TMT; Skidmore et al. 2015) both include IFS instruments with resolutions of 4
milliarcseconds (HARMONI and IRIS, respectively) as part of their first
generation suite of instruments. Under the $\Lambda$CDM cosmology of Planck
Collaboration et al. (2020b), this angular resolution corresponds to datacubes
with a physical resolution finer than 35 parsec per pixel at _all redshifts_.
Similarly, the IFS instrument GMTIFS planned for the Giant Magellan Telescope
(GMT; Sharp et al. 2018) will have an angular resolution of 6 milliarcseconds,
which corresponds to a physical resolution of 52 parsec per pixel or finer at
all redshifts. Based on the resolution analysis presented here, any one of
these instruments will produce dataset suited for geostatistical analysis that
will characterise the two-dimensional chemical structure of galaxies at high
redshift to a high level of detail. Comparing such galaxies to local ones, it
should be possible to determine whether the process of metal mixing is the
same in these chaotic starbursting systems as they are in the Universe today.
In turn, this will contribute to understanding whether the processes that
govern stellar feedback are constant throughout cosmic time. With existing
high redshift metallicity diagnostics (e.g Sanders et al., 2023; Laseter et
al., 2023; Nakajima et al., 2022), based on the wavelength coverage of these
instruments, such analysis will be possible for galaxies up to $z\approx 4$.
Progressing beyond this redshift limit is possible if new rest-frame UV
metallicity diagnostics for high redshift galaxies are developed.

Overall, our analysis highlighted the applicability and benefits of a
geostatistical analysis. As with any novel diagnostic, the results obtained
presents interesting aspects for further research. In particular, at the
finest resolutions ($\lesssim 50$ pc per pixel), the best-fit geostatistical
models generally tended to fit weaker metallicity gradients than was seen for
data at coarser resolutions. This behaviour is not intuitive, as other studies
have shown that binning data to coarser resolutions leads to flatter recovered
gradients for galaxies (Yuan et al., 2013; Mast et al., 2014; Acharyya et al.,
2020). However, these studies have focused on simulating galaxy observations
with resolutions far coarser than the range in which this effect was seen to
be important.

One interpretation of this effect is that mean metallicity gradients of
galaxies are an emergent property that comes about from averaging over small-
scale fluctuations. This would imply that the global metallicity properties of
galaxies are set by their small-scale gas physics – however, Baker et al.
(2023) argue that kpc-scale features alone cannot be responsible for
completely determining a galaxy’s metallicity. Another interpretation is that
the metallicity gradient of a galaxy is not an intrinsic property of the
galaxy itself, and that higher resolution data reveals significant departures
from a linear trend that require a different model to produce a good fit.

We find that a flexible geostatistical model accounting for small-scale
metallicity fluctuations within the data produces very good predictions for
the metallicity of unknown regions within a galaxy, especially with high
resolution data. Obtaining accurate metallicity predictions at specific points
of a galaxy is important not only for studies of galaxy evolution and ISM
physics, but also for cosmology. Accurate estimates of the metallicity at the
site of many standard candles is required in order to calibrate their period-
luminosity relations (e.g. Ripepi et al., 2020; Bhardwaj et al., 2023).
Therefore, the predictive accuracy of metallicity models is directly linked to
the precision of local Universe measurements of the Hubble constant (Riess et
al., 2022).

The circularly-symmetric metallicity gradient model has been shown to be a
good descriptor of the internal metallicity structure of local disc galaxies
on the kpc-scale (e.g. Sánchez et al., 2014; Bresolin & Kennicutt, 2015). For
these galaxies, such a model is well-motivated by theories of inside-out
galaxy evolution. However, its applicability to merging systems, dwarf
galaxies, and galaxies that display large asymmetries in their surface
brightness profile is less well motivated. At high-redshift, galaxies are
known to be more clumpy and irregular than local galaxies (Conselice, 2014),
suggesting that models with more flexibility than a linear gradient model may
be important for understanding the baryon cycle at cosmic noon and earlier
epochs. Understanding the data quality and resolution required to make a
meaningful fit to a galaxy’s metallicity map using a geostatistical model is
an important first step to determining whether these models will be useful
tools to analyse data from future IFS surveys of irregular systems and
galaxies at high redshift.

In this work, the physical resolution of the data was defined to be the
angular size of spaxels in psuedo-IFS data, times the distance to the galaxy.
Distances to galaxies in the local universe can be highly uncertain. For
example, we follow Leroy et al. (2019) in quoting the distance to NGC 6744 as
$11.6$ Mpc, but other references have found the distance to this galaxy to be
as close as $7.6\pm 1.7$ Mpc (Sorce et al., 2014). Such large discrepancies in
the distances to these galaxies make the values of $\phi$ much more uncertain
than the values that we report – for example, adopting the distance from Sorce
et al. (2014) would decrease the value of $\phi$ measured for NGC 6744 from
$1.3$ kpc to $852$ pc. However, all results we have found on the relationship
between the size of small-scale fluctuations and the resolution of the data
would remain valid, as these two quantities scale with distance in a perfectly
correlated way.

An alternative definition of the spatial resolution of an IFS datacube is
based on the FWHM of the PSF for that data. For TYPHOON data, the PSF is
several times smaller than the angular size of a spaxel. Because it is
negligible at native resolution, we do not simulate the effects of having a
large PSF that extends over many pixels at any coarser resolutions, for
consistency. When the PSF is much larger than the size of a spaxel in an IFS
observation, metallicity gradients appear flatter (Acharyya et al., 2020), and
small-scale features may be blurred or lost. Techniques such as forward
modelling can be used to model the effects of poor seeing on large-scale
galaxy properties and correct for them (Metha et al. in prep.), and help
distinguish between physical small-scale features and small-scale features
caused by a known PSF profile. Additionally, correlations expected between
nearby spaxels caused by instrumental effects could be directly incorporated
into a geostatistical model by changing the form of $\epsilon(\vec{x})$ to
capture this correlated structure.

Another important effect that was not directly controlled for in the series of
resolution experiments presented in this work is the effect of altering the
signal to noise ratio, S/N, of the data. Metallicity maps with coarser
resolution will have a greater S/N than those at finer resolution due to
binning. Presumably, this would affect the performance of a model-fitting
procedure. We note that for all resolutions tested, the same cut of S/N$>3$ on
all lines used for metallicity diagnostics was enforced, ensuring that even at
our finest resolution, no poor quality emission line fits are used to
determine metallicities. We further note that geostatistical methods are
suitable for distinguishing uncorrelated sources of noise, such as Poisson
noise in a detector, from correlated signals of small-scale metallicity
fluctuations. Thus, a geostatistical data analysis may be less impacted by
having low S/N data than some other traditional statistical analysis
techniques would be (see M21 for further details).

Finally, the inclination of a galaxy, $i$, is expected to make identification
and classification of galaxy small-scale features more difficult, by
effectively decreasing the resolution along one axis of a galaxy’s plane.
Understanding how the geostatistical models presented in this work are
affected by inclination is important to determine which systems these
techniques are best applied to, and for choosing targets for future high-
resolution IFS surveys. Because a galaxy can only be seen from one inclination
angle, and small-scale properties of the ISM are not uniform between galaxies,
mock-IFS observations of simulated galaxies are required to investigate this
effect. We leave investigations on this subject for future work.

##  7 Summary and conclusions

In this series of papers, we introduce a novel framework for interpretation
and modelling of astronomical datacubes based on spatial statistics
techniques. In M21, geostatistics was introduced as a useful tool for
extracting additional data from high-resolution metallicity maps, allowing
local deviations from large-scale metallicity profiles to be detected,
classified, and compared to theoretical models of metal mixing in a galaxy’s
ISM. In M22, we extended on this theme, showing how geostatistical models can
be used to make accurate predictive maps of the metallicity structure of a
galaxy in 2-dimensions, including within regions where DIG contamination
prevents the metallicity from being directly measurable. Here, we investigate
how the capability of geostatistical methods to observe, model, and predict
small-scale features in the metallicity maps of galaxies are affected by data
quality, to increase the sample size of (real and simulated) galaxies over
which these methods can be applied, and aid the design of future surveys that
seek to provide detailed descriptions of the ISM of a wide variety of
galaxies. In particular, we focus on how the best-fit physical parameters and
predictive power of a geostatistical model are affected by the resolution and
completeness of an IFS-like observation. We summarise our main results below:

  * •

We found that $\sim 100$ data points with known metallicity are required to
train a geostatistical model of metallicity variations within a galaxy’s ISM.
Below this value, noise due to small-number statistics makes computing a
semivariogram very difficult, making it difficult to infer the presence of
small-scale features. Furthermore, predictions from a geostatistical model
perform no better than a linear gradient model. Increasing the number of data
points available for a geostatistical analysis above this value increases the
accuracy of kriging predictions, and decreases the uncertainty on the model’s
parameters.

  * •

The resolution required to fit a geostatistical model that accurately captures
small-scale features in the metallicity distributions of galaxies was found to
depend on the size of those features. These features are typically of the
order of $200-300$ pc across, but vary from galaxy to galaxy. As a general
rule, we find that the geostatistical data pipeline presented in this work is
able to capture any small-scale features that are apparent on visual
inspection; and from the galaxies explored in this work, we find that this
tends to be achieved when the resolution is better than $\sim 0.1R_{e}$.

  * •

For coarsely resolved data, the best fit value of $\phi$, the typical spatial
extent of metallicity fluctuations, was found to depend on the resolution of
the data. Therefore, best-fit values of $\phi$ should be trusted if and only
if the size of a pixel is less than $\phi/2$. In these cases, predictions from
a geostatistical model vastly outperformed a metallicity gradient model.

  * •

Even when the resolution of the data was not sufficient to accurately
constrain $\phi$ to its true value, kriging predictions were still found to be
much more accurate than predictions from a linear gradient alone for datasets
with $\gtrsim 100$ data points, and pixel sizes $<\phi$ where small-scale
structure is marginally resolvable by eye. When fewer data points are
available, or the resolution is coarser, predictions from kriging offered only
a mild improvement over predictions from a linear metallicity gradient.

  * •

A tradeoff exists between using integrated Hii regions and individual Hii-
dominated spaxels for a spatially-resolved metallicity analysis. Integrating
spaxels into Hii regions decreases the number of data points available and the
resolution of the data, but also reduces the uncertainty of the metallicity of
each data point. Ultimately, we found that both methods yielded consistent
results. Therefore, we recommend binning data points into Hii regions before a
spatial analysis only when the average sizes of Hii regions stays below
$\sim\phi/2$, and the number of data points stays above $\sim 100$.

Based on these results, there already exists a large legacy data set of local
universe galaxy observations that have sufficient resolution and a sufficient
number of data points for geostatistical models to be applied (López-Cobá et
al., 2020). Key questions to answer include understanding how the metal mixing
scale of galaxies varies with their size, morphology, star formation rate, and
environment, to illuminate the interplay between small-scale gas physics and
galaxy evolution. Comparisons to zoom-in galaxy simulations will allow
detailed models of ISM mixing to be directly compared to observational data.
Future observations with high-resolution IFS instruments, such as the first
generation IFS instruments planned for E-ELT, GMT, and TMT, will allow
geostatistical methods to be applied to galaxies at redshifts up to $z\sim 4$,
revealing how turbulent mixing processes within galaxies have changed over
cosmic time.

## acknowledgements

We would like to thank the anonymous referee, whose comments helped to improve
the quality of this paper. BM acknowledges support from an Australian
Government Research Training Program (RTP) Scholarship. This research is
supported in part by the Australian Research Council Centre of Excellence for
All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number
CE170100013. BM would like to thank Prof. Tommaso Treu for his mentorship and
advice on using MCMC methods and for hosting a wonderful extended research
visit at UCLA. The authors extend their thanks to Alejandra Lugo-Aranda for
the advice on using pyHIIextractor, and to Dr. Kathryn Grasha for the advice
on using HIIphotand other helpful conversations. This research was conducted
on Wurundjeri land.

## data availability

The TYPHOON data will be publicly released in Seibert et al. (in prep). Prior
to this release, the TYPHOON data are available upon reasonable request to
Andrew Battisti. Any further data products created for this work are available
from the corresponding author (BM) upon reasonable request. A publically
available Python package containing implementations of the geostatistical
techniques explored in this series is currently in preparation (Metha et al.
in prep.).

## References

  * Acharyya et al. (2020) Acharyya A., Krumholz M. R., Federrath C., Kewley L. J., Goldbaum N. J., Sharp R., 2020, MNRAS, 495, 3819
  * Andrews et al. (2020) Andrews S., Fryer C., Even W., Jones S., Pignatari M., 2020, ApJ, 890, 35
  * Arellano-Córdova et al. (2022) Arellano-Córdova K. Z., et al., 2022, arXiv e-prints, p. arXiv:2206.04280
  * Baker et al. (2023) Baker W. M., et al., 2023, MNRAS, 519, 1149
  * Baldwin et al. (1981) Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
  * Battisti et al. (2017) Battisti A. J., Calzetti D., Chary R. R., 2017, ApJ, 851, 90
  * Bec & Khanin (2007) Bec J., Khanin K., 2007, Phys. Rep., 447, 1
  * Belfiore et al. (2017) Belfiore F., et al., 2017, MNRAS, 469, 151
  * Belley & Roy (1992) Belley J., Roy J.-R., 1992, ApJS, 78, 61
  * Berg et al. (2013) Berg D. A., Skillman E. D., Garnett D. R., Croxall K. V., Marble A. R., Smith J. D., Gordon K., Kennicutt Robert C. J., 2013, ApJ, 775, 128
  * Bhardwaj et al. (2023) Bhardwaj A., et al., 2023, ApJ, 944, L51
  * Boissier & Prantzos (1999) Boissier S., Prantzos N., 1999, MNRAS, 307, 857
  * Bresolin & Kennicutt (2015) Bresolin F., Kennicutt R. C., 2015, MNRAS, 454, 3664
  * Bresolin et al. (2009) Bresolin F., Ryan-Weber E., Kennicutt R. C., Goddard Q., 2009, ApJ, 695, 580
  * Bryant et al. (2015) Bryant J. J., et al., 2015, MNRAS, 447, 2857
  * Bundy et al. (2015) Bundy K., et al., 2015, ApJ, 798, 7
  * Cardelli et al. (1989) Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  * Chen et al. (2023) Chen Q.-H., Grasha K., Battisti A. J., Kewley L. J., Madore B. F., Seibert M., Rich J. A., Beaton R. L., 2023, MNRAS, 519, 4801
  * Christensen et al. (2018) Christensen C. R., Davé R., Brooks A., Quinn T., Shen S., 2018, ApJ, 867, 142
  * Conselice (2014) Conselice C. J., 2014, ARA&A, 52, 291
  * Cosens et al. (2018) Cosens M., et al., 2018, ApJ, 869, 11
  * Cressie (1990) Cressie N., 1990, Mathematical Geology, 22, 239
  * Cressie (1993) Cressie N., 1993, Statistics for Spatial Data, revised edn. Wiley, New York
  * Croom et al. (2012) Croom S. M., et al., 2012, MNRAS, 421, 872
  * Curti et al. (2020) Curti M., Mannucci F., Cresci G., Maiolino R., 2020, MNRAS, 491, 944
  * Dopita et al. (2016) Dopita M. A., Kewley L. J., Sutherland R. S., Nicholls D. C., 2016, Ap&SS, 361, 61
  * Emsellem et al. (2022) Emsellem E., et al., 2022, A&A, 659, A191
  * Erroz-Ferrer et al. (2019) Erroz-Ferrer S., et al., 2019, MNRAS, 484, 5009
  * Espinosa-Ponce et al. (2020) Espinosa-Ponce C., Sánchez S. F., Morisset C., Barrera-Ballesteros J. K., Galbany L., García-Benito R., Lacerda E. A. D., Mast D., 2020, MNRAS, 494, 1622
  * Falceta-Gonçalves et al. (2014) Falceta-Gonçalves D., Kowal G., Falgarone E., Chian A. C. L., 2014, Nonlinear Processes in Geophysics, 21, 587
  * Feigelson & Babu (2012) Feigelson E. D., Babu G. J., 2012, Modern Statistical Methods for Astronomy, doi:10.48550/arXiv.1205.2064.
  * Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Field et al. (2012) Field A., Miles J., Field Z., 2012], Discovering Statistics using R. SAGE Publications, London, England
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gelfand (2010) Gelfand A. E., 2010, in Gelfand A. E., Diggle P. J., Fuentes M., Guttorp P., eds, Handbook of Spatial Statistics. Taylor & Francis Group, 6000 Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487-2742, pp 517–539
  * Gilmozzi & Spyromilio (2007) Gilmozzi R., Spyromilio J., 2007, The Messenger, 127, 11
  * Goddard et al. (2011) Goddard Q. E., Bresolin F., Kennicutt R. C., Ryan-Weber E. V., Rosales-Ortega F. F., 2011, MNRAS, 412, 1246
  * Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
  * Grasha et al. (2022) Grasha K., et al., 2022, ApJ, 929, 118
  * Groves et al. (2023) Groves B., et al., 2023, MNRAS, 520, 4902
  * Haffner et al. (2009) Haffner L. M., et al., 2009, Reviews of Modern Physics, 81, 969
  * Heckman et al. (1990) Heckman T. M., Armus L., Miley G. K., 1990, ApJS, 74, 833
  * Ho et al. (2015) Ho I. T., Kudritzki R.-P., Kewley L. J., Zahid H. J., Dopita M. A., Bresolin F., Rupke D. S. N., 2015, MNRAS, 448, 2030
  * Ho et al. (2016) Ho I. T., et al., 2016, Ap&SS, 361, 280
  * Ho et al. (2017) Ho I. T., et al., 2017, ApJ, 846, 39
  * Ho et al. (2018) Ho I. T., et al., 2018, A&A, 618, A64
  * Ho et al. (2019) Ho I. T., et al., 2019, ApJ, 885, L31
  * Hogg & Villar (2021) Hogg D. W., Villar S., 2021, PASP, 133, 093001
  * Jarrett et al. (2003) Jarrett T. H., Chester T., Cutri R., Schneider S. E., Huchra J. P., 2003, AJ, 125, 525
  * Jin et al. (2022) Jin Y., Kewley L. J., Sutherland R. S., 2022, ApJ, 934, L8
  * Journel & Huijbregts (1978) Journel A., Huijbregts C. J., 1978, Mining Geostatistics. Academic Press
  * Kaplan et al. (2016) Kaplan K. F., et al., 2016, MNRAS, 462, 1642
  * Kauffmann et al. (2003) Kauffmann G., et al., 2003, MNRAS, 346, 1055
  * Kewley & Ellison (2008) Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183
  * Kewley et al. (2001) Kewley L. J., Heisler C. A., Dopita M. A., Lumsden S., 2001, ApJSupplement Series, 132, 37
  * Kewley et al. (2019) Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511
  * Kreckel et al. (2019) Kreckel K., et al., 2019, ApJ, 887, 80
  * Kreckel et al. (2020) Kreckel K., et al., 2020, MNRAS, 499, 193
  * Kumari et al. (2019) Kumari N., Maiolino R., Belfiore F., Curti M., 2019, MNRAS, 485, 367
  * Larson (1972) Larson R. B., 1972, Nature, 236, 21
  * Laseter et al. (2023) Laseter I. H., et al., 2023, arXiv e-prints, p. arXiv:2306.03120
  * Lee & Gammie (2021) Lee D., Gammie C. F., 2021, ApJ, 906, 39
  * Leroy et al. (2019) Leroy A. K., et al., 2019, ApJS, 244, 24
  * Li et al. (2021) Li Z., Krumholz M. R., Wisnioski E., Mendel J. T., Kewley L. J., Sánchez S. F., Galbany L., 2021, MNRAS, 504, 5496
  * Li et al. (2022) Li Z., et al., 2022, arXiv e-prints, p. arXiv:2206.08072
  * López-Cobá et al. (2020) López-Cobá C., et al., 2020, AJ, 159, 167
  * Lugo-Aranda et al. (2022) Lugo-Aranda A. Z., Sánchez S. F., Espinosa-Ponce C., López-Cobá C., Galbany L., Barrera-Ballesteros J. K., Sánchez-Menguiano L., Anderson J. P., 2022, arXiv e-prints, p. arXiv:2204.04252
  * Mannucci et al. (2021) Mannucci F., et al., 2021, MNRAS, 508, 1582
  * Mast et al. (2014) Mast D., et al., 2014, A&A, 561, A129
  * Matheron (1954) Matheron G., 1954, Note géostatistique, 1
  * Matheron (1963) Matheron G., 1963, Economic geology, 58, 1246
  * Metha et al. (2021) Metha B., Trenti M., Chu T., 2021, MNRAS, 508, 489
  * Metha et al. (2022) Metha B., Trenti M., Chu T., Battisti A., 2022, MNRAS, 514, 4465
  * Moustakas & Kennicutt (2006) Moustakas J., Kennicutt Robert C. J., 2006, ApJ, 651, 155
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Nakajima et al. (2022) Nakajima K., et al., 2022, ApJS, 262, 3
  * Persson et al. (1990) Persson S. E., Carr D. M., Jacobs J. H., 1990, Experimental Astronomy, 1, 195
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Pilyugin & Grebel (2016) Pilyugin L. S., Grebel E. K., 2016, MNRAS, 457, 3678
  * Planck Collaboration et al. (2020a) Planck Collaboration et al., 2020a, A&A, 641, A5
  * Planck Collaboration et al. (2020b) Planck Collaboration et al., 2020b, A&A, 641, A6
  * Poetrodjojo et al. (2019) Poetrodjojo H., D’Agostino J. J., Groves B., Kewley L., Ho I. T., Rich J., Madore B. F., Seibert M., 2019, MNRAS, 487, 79
  * Poetrodjojo et al. (2021) Poetrodjojo H., et al., 2021, MNRAS, 502, 3357
  * Reichardt Chu et al. (2022) Reichardt Chu B., et al., 2022, ApJ, 941, 163
  * Reynolds (1990) Reynolds R. J., 1990, in Bowyer S., Leinert C., eds, Proceedings of the 139th. Symposium of the International Astronomical Union Vol. 139, The Galactic and Extragalactic Background Radiation. p. 157
  * Riess et al. (2022) Riess A. G., et al., 2022, ApJ, 934, L7
  * Ripepi et al. (2020) Ripepi V., et al., 2020, A&A, 642, A230
  * Rosales-Ortega et al. (2011) Rosales-Ortega F. F., Díaz A. I., Kennicutt R. C., Sánchez S. F., 2011, MNRAS, 415, 2439
  * Sánchez (2020) Sánchez S. F., 2020, ARA&A, 58, 99
  * Sánchez-Menguiano et al. (2016) Sánchez-Menguiano L., et al., 2016, ApJ, 830, L40
  * Sánchez-Menguiano et al. (2018) Sánchez-Menguiano L., et al., 2018, A&A, 609, A119
  * Sánchez et al. (2012) Sánchez S. F., et al., 2012, Astronomy & Astrohpysics, 538, A8
  * Sánchez et al. (2014) Sánchez S. F., et al., 2014, A&A, 563, A49
  * Sanders et al. (2023) Sanders R. L., Shapley A. E., Topping M. W., Reddy N. A., Brammer G. B., 2023, arXiv e-prints, p. arXiv:2303.08149
  * Schönrich & McMillan (2017) Schönrich R., McMillan P. J., 2017, MNRAS, 467, 1154
  * Searle (1971) Searle L., 1971, ApJ, 168, 327
  * Sharp et al. (2018) Sharp R., et al., 2018, in Evans C. J., Simard L., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII. p. 107021V, doi:10.1117/12.2313010
  * Skidmore et al. (2015) Skidmore W., TMT International Science Development Teams Science Advisory Committee T., 2015, Research in Astronomy and Astrophysics, 15, 1945
  * Sorce et al. (2014) Sorce J. G., Tully R. B., Courtois H. M., Jarrett T. H., Neill J. D., Shaya E. J., 2014, MNRAS, 444, 527
  * Stein (1999) Stein M., 1999, Interpolation of spatial data: Some theory for kriging. Springer Verlag, New York
  * Sutherland et al. (2018) Sutherland R., Dopita M., Binette L., Groves B., 2018, MAPPINGS V: Astrophysical plasma modeling code, Astrophysics Source Code Library, record ascl:1807.005 (ascl:1807.005)
  * Tobler (1970) Tobler W. R., 1970, Economic Geography, 46, 234
  * Trapp et al. (2022) Trapp C. W., et al., 2022, MNRAS, 509, 4149
  * Vila-Costas & Edmunds (1992) Vila-Costas M. B., Edmunds M. G., 1992, MNRAS, 259, 121
  * Vogt et al. (2017) Vogt F. P. A., Pérez E., Dopita M. A., Verdes-Montenegro L., Borthakur S., 2017, A&A, 601, A61
  * Wang et al. (2022) Wang X., et al., 2022, arXiv e-prints, p. arXiv:2207.13113
  * Wetzel et al. (2023) Wetzel A., et al., 2023, ApJS, 265, 44
  * Wikle et al. (2019) Wikle C. K., Zammit-Mangion A., C. C. N. A., 2019, Spatio-temporal statistics with R. CRC Press, Taylor &; Francis Group
  * Williams et al. (2022) Williams T. G., et al., 2022, MNRAS, 509, 1303
  * Yuan et al. (2013) Yuan T. T., Kewley L. J., Rich J., 2013, ApJ, 767, 106
  * Zaritsky et al. (1994) Zaritsky D., Kennicutt Robert C. J., Huchra J. P., 1994, ApJ, 420, 87
  * Zhang et al. (2017) Zhang K., et al., 2017, MNRAS, 466, 3217
  * Zinchenko et al. (2016) Zinchenko I. A., Pilyugin L. S., Grebel E. K., Sánchez S. F., Vílchez J. M., 2016, MNRAS, 462, 2715
  * de Vaucouleurs et al. (1991) de Vaucouleurs G., de Vaucouleurs A., Corwin Herold G. J., Buta R. J., Paturel G., Fouque P., 1991, Third Reference Catalogue of Bright Galaxies

##  Appendix A Analysis with alternative metallicity diagnostics

Estimating metallicities from strong emission line diagnostics has always been
fraught with systematics. Different metallicity diagnostics can give
metallicity estimates that differ by as much as 0.7 dex (Kewley & Ellison,
2008). Furthermore, metallicity gradients measured in galaxies using different
metallicity diagnostics can be different as well (Poetrodjojo et al., 2021).

For this reason, care must be taken to ensure that the findings of this study
are not dependent on the specific metallicity diagnostic used. To ensure our
results are robust, we test three additional metallicity diagnostics, using
three different emission line ratios, calibrated in a variety of ways. The
first of these is the N2S2H$\alpha$ diagnostic, together with the theoretical
calibration of Dopita et al. (2016) devised using the plasma modelling code
Mappings 5.0 (Sutherland et al., 2018). This diagnostic uses three emission
lines: [Nii$]\lambda$6583, H$\alpha$, and the doublet
[Sii$]\lambda\lambda$6717,6731. This diagnostic has been shown to be largely
insensitive to both dust extinction and variations in the ionisation
parameter.

Secondly, we consider the S-calibration (Scal hereafter) of Pilyugin & Grebel
(2016). This empirical calibration was based on 313 reference Hii regions with
metallicities determined using the direct $T_{e}$-based method, and uses the
emission lines H$\beta$, [Nii$]\lambda$6583, [Oiii$]\lambda$5007, and the
doublet [Sii$]\lambda\lambda$6717,6731.

Finally, we examine how the results change when the O3N2 calibration of Curti
et al. (2020) is used. This diagnostic was calibrated using exactly the same
methods as our fiducial diagnostic, RS32. However, unlike the other
diagnostics discussed in this work, O3N2 has no dependence on sulphur, instead
depending only on H$\alpha$, H$\beta$, [Oiii$]\lambda$5007 and
[Nii$]\lambda$6853\. Thus, it provides an independent test against effects
that could artificially raise or lower the relative strength of sulphur lines,
such as DIG contamination (Haffner et al., 2009).

Figure 18: Medians, 16${}^{\text{th}}$, and 84${}^{\text{th}}$ percentiles of
the posterior distributions of the parameters of the geostatistical model
recovered for NGC 5236, when N2S2H$\alpha$ is chosen to be the metallicity
diagnostic. At the finest resolutions, a model with a metallicity gradient
close to zero is preferred. A value of $\phi\approx 200$ pc is preferred when
this diagnostic is used; and this value can be accurately recovered until the
pixel size becomes comparable to this value.

We show how the fitted values of $\phi,$ $Z_{\text{char}}$, and $\nabla Z$
vary with resolution for NGC 5236 when the N2S2H$\alpha$ diagnostic is used in
Figure 18. When this diagnostic is used, metallicity fluctuations are found to
be correlated on a slightly smaller scale than was found when RS32 was used,
with $\phi=204\pm 12$ pc at native resolution. A consistent value is found
when $f=2$, but for coarser resolutions with $f=4$ and $f=6$, this value
increases slightly, to $\sim 250$ pc. At resolutions coarser than this, the
recovered values of $\phi$ become increasingly high, as the size of a pixel
becomes larger that the size of these small-scale fluctuations. We note a
trend for the size of the fluctuations recovered to match the size of a pixel
for four of the models fitted. This makes sense – our binning process removes
all signal of correlation below the size of a single pixel, but some small-
scale correlations are still observable in the data, so the model assumes that
the scale of the correlation and the scale of a single pixel are the same.

It is unclear to us why the recovered size of metallicity fluctuations are
different when this diagnostic is used, as compared to RS32. Nonetheless, from
this sequence of model fits, we see that to accurately recover the value of
$\phi$, the size of a pixel must be smaller than $\sim\phi/2$ – the same
result that is seen for the two other galaxies discussed in Section 5.

We also notice the same effect seen in Section 4.1, where the recovered value
of $\nabla Z$ is close to zero at our two finest resolutions. At coarser
resolutions than this, the value of $\nabla Z$ recovered is similar to, yet
slightly steeper than, the metallicity gradient observed when our fiducial
RS32 metallicity diagnostic is used. Disparities in the metallicity gradient
recovered when different metallicity diagnostics are used have been seen
before for this galaxy (e.g Poetrodjojo et al., 2019), and the size of the
difference in our study ($\sim 0.005$ dex kpc-1) is consistent with such
expectations of variability across diagnostic tools.

Figure 19: As in Figure 10, but using N2S2H$\alpha$ as a metallicity
diagnostic. We see that a geostatistical model fit still outperforms a linear
gradient method at all resolutions tested.

In Figure 19, we show the comparative predictive accuracy of our
geostatistical model against a model that only contains a metallicity gradient
at all resolutions tested. We see that our main results found when RS32 was
chosen to be our metallicity diagnostic also hold for N2S2H$\alpha$: namely,
(i) the geostatistical model outperforms a standard metallicity gradient model
at all resolutions; (ii) the improvement in the fit is best seen for high-
resolution data; and (iii) even when the recovered value of $\nabla Z$ is
close to zero for finely-resolved data, the predictions of the geostatistical
model are still very accurate.

Figure 20: Medians, 16${}^{\text{th}}$, and 84${}^{\text{th}}$ percentiles of
the posterior distributions of the parameters of the geostatistical model
recovered for NGC 5236, when Scal is chosen to be the metallicity diagnostic.
Values of $\phi$ recovered become unreliable when the pixel size approaches
the size of the intrinsic fluctuations. At the finest resolution, the best-fit
value of $\nabla Z$ is close to zero, consistent with our other results.

Our results when the Scal diagnostic is used are very similar to the results
found when N2S2H$\alpha$ is used. This is to be expected, as these two
metallicity diagnostics have been found to be highly correlated (Groves et
al., 2023). We show how the fitted values change as we change the resolution
of our data in Figure 20. At native resolution, $\phi=222^{+18}_{-14}$ pc,
consistent with the value found when N2S2H$\alpha$ is used. Consistent values
of $\phi$ are observed until the size of a pixel becomes larger than $200$pc.
After this point, $\phi$ is consistently estimated to be equal to the size of
one pixel, until $f=24$ where no small-scale structure can be found, and the
best-fit value of $\phi$ is only recovered with very large uncertainties. As
is found for the other two diagnostics tested, at the finest resolution, the
geostatistical model prefers to fit a metallicity gradient that is close to
zero.

We plot the predictive accuracy of a geostatistical model for NGC 5236 when
this diagnostic is used compared to a metallicity-gradient only model in
Figure 21. We see the following results: (i) the geostatistical model always
outperforms a linear fit; (ii) when the resolution is so coarse that $\phi$
cannot be accurately recovered, the improvement over a linear gradient fit
does not improve as the resolution of the data improves; (iii) when the data
is of a fine enough resolution for $\phi$ to be recovered accurately, the
predictive accuracy of the model continues to increase as the resolution is
made finer, all the way down to native resolution; and (iv) having a recovered
gradient close to zero at the finest resolution does not affect the predictive
performance of this model.

Figure 21: As in Figure 10, but using Scal as a metallicity diagnostic.
Predictions from the geostatistical model always outperform those computed
from the linear gradient model. The improvement is most significant at
resolutions finer than $\sim 200$pc, comparable to the value of $\phi$
measured for this galaxy using this diagnostic. Figure 22: Medians,
16${}^{\text{th}}$, and 84${}^{\text{th}}$ percentiles of the posterior
distributions of the parameters of the geostatistical model recovered for NGC
5236, when O3N2 is chosen to be the metallicity diagnostic. When $f=6$, small-
scale structures stop becoming apparent in the data, as evidenced by the large
value of $\phi$. Around this point, the best-fit value of $\nabla Z$ decreases
and $Z_{\text{char}}$ increases. For smaller and larger resolutions, the trend
of $\phi$ increasing with resolution is similar to what is seen in other
diagnostics. Figure 23: As in Figure 10, but using O3N2 as a metallicity
diagnostic. We find that when $f=6$, a metallicity gradient alone provides an
extremely good predictive model for the data. Even so, at all resolutions
including when $f=6$, the linear gradient model is outperformed by a
geostatistical model.

Finally, we show the results when the O3N2 diagnostic is used. In Figure 22,
we plot the median, 16th, and 84th percentiles in the posteriors of each
parameter computed from a Monte-Carlo fit as a function of the resolution of
the data input. We find that at native resolution, $\phi=224^{+16}_{-15}$ pc,
and $\nabla Z$ is flatter than is found at coarser resolutions, consistent
with the other two diagnostics discussed in this Appendix. Results are similar
when $f=2$. However, an anomalous result is seen when $f=6$, and possibly at
the two resolutions closest to this resolution. For this particular bin
factor, using this diagnostic, the best-fit value of $\phi$ becomes very
large, with large error bars: $\phi=1.63^{+1.20}_{-0.51}$ kpc. This is larger
than the value seen at any other resolution, including when $f=24$. Non-
monotonic trends are also seen in the recovered values of $\nabla Z$ and
$Z_{\text{char}}$ versus resolution when this diagnostic is used: at $f=6$,
the value of $\nabla Z$ reaches a minimum, and the (negatively correlated;
$\rho=-0.57$) value of $Z_{\text{char}}$ reaches a maximum. This indicates
that for this particular resolution, with this particular diagnostic, the
variance in the data can be preferentially described by a strong negative
metallicity gradient, with no evidence for sub-kpc structures.

Further evidence for this interpretation is seen in Figure 23. When $f=6$, the
median absolute difference between predicted metallicities from a linear
gradient model and observed metallicities drops to $0.4\sigma$ – significantly
smaller than the WLS model error seen when any other diagnostic is used to
analyse NGC 5236 at any resolution, and half of the value of MAD/$\sigma$ that
is expected from Gaussian error. This indicates that for this combination of
resolution and metallicity diagnostic, a linear gradient model does an
excellent job of describing the metallicity variations within the data.
Interestingly, the best-fit metallicity gradient found using a WLS method does
not change significantly between $f=1$ and $f=10$ ($-0.012\leq\nabla Z\leq
0.09$ dex kpc-1), indicating that this improvement in the performance of the
model at this resolution is not caused by a more accurate recovery of the
parameters, but by some difference in the properties of the data when this
particular binning is used.

Even at $f=6$, when the data is very well-described by a simple linear trend,
we find that our geostatistical model still outperforms the linear fit, with a
lower MAD/$\sigma$ than is found even for the native resolution O3N2
metallicity map. We do not see any similar behaviour of the O3N2 diagnostic
producing highly linear metallicity maps at particular resolutions for the
other two galaxies investigated in this work. Because of this, we do not
speculate as to what aspects of the O3N2 diagnostic may be responsible for
this unusual behaviour.

Clearly, there is much work left to be done in understanding how strong
emission line based metallicity diagnostics calibrated on stacked, integrated
galaxy spectra may be applied to sub-kpc regions of the ISM. Extending on the
conclusions of Kewley & Ellison (2008), we caution the reader that relying on
individual strong emission line diagnostics may lead to erroneous conclusions
about the metallicity structure of a galaxy. Following Kewley et al. (2019),
we recommend using a combination of multiple (three or more) strong emission
line diagnostics using different line ratios to achieve a robust understanding
of a galaxy’s internal metallicity distribution.

##  Appendix B A primer on geostatistical methods

To capture details on the small-scale structure of the metallicity data around
a large-scale metallicity trend, this work resorts to tools and techniques
from _geostatistics_ , the field of mathematics and statistics that concerns
itself with analysing stochastic processes that occur over a continuous
spatial domain. This field is vast, having evolved over several decades since
the seminal works of Matheron (1954), and we do not attempt to cover a
comprehensive review here.555For those interested in geostatistics who want to
read further, we recommend Wikle et al. (2019) as an introductory text, or
Cressie (1993) for a comprehensive reference. Rather, we introduce two
important concepts that have been discussed extensively in the previous papers
in this series: the semivariogram, a tool for exploratory data analysis and
separating correlated signals from uncorrelated noise, introduced for
astronomical applications in M21; and universal kriging, an interpolation
method that takes as input a set of data points and a geostatistical model,
and provides as output predictions of the value of a random field at a set of
specified locations (discussed extensively in M22). Furthermore, we provide a
comparison between the semivariogram, the two point correlation function and
the power spectrum – two popular tools in the analysis of spatial structure in
astronomical datasets (e.g. Planck Collaboration et al., 2020a; Li et al.,
2021).

###  B.1 The semivariogram and related methods

To capture details on the small-scale structure of the metallicity data around
the metallicity trend, we use the semivariogram, a classic tool for
exploratory data analysis from the geostatistical literature (Matheron, 1963).
The semivariogram is defined as:

$\gamma(h)=\frac{1}{2}\text{Var}\left(Z(\vec{x})-Z(\vec{y})\right),$ (5)

where the variance is computed over all pairs of points $\vec{x}$ and
$\vec{y}$ for which $h-\delta_{h}/2\leq|\vec{x}-\vec{y}|\leq h+\delta_{h}/2$.
Here, $\delta_{h}$ is the size of bins used to compute the semivariogram,
which should be chosen to be large enough so that enough pairs are present at
each separation of interest that a variance can be reliably computed, but as
small as possible to retain high spatial resolution (e.g. Cressie, 1993;
Journel & Huijbregts, 1978).

Intuitively, the semivariogram shows how the variance between data points
increases as the separation between them increases. Three informative
parameters can be read from a semivariogram: the range, which is the
separation between data points at which the variance ceases to increase,
indicative of the sizes of the largest substructures within the data; the
sill, which is the height of the semivariogram after it has flattened off,
revealing the total variance in the data; and the nugget, which is the height
of the semivariogram at zero separation, equivalent to the amount of
uncorrelated noise in the data. By comparing the height of the sill to the
nugget, we can quantify the amount of correlated versus uncorrelated variance
in the data, allowing the amplitude, in dex, of metallicity fluctuations to be
inferred while ignoring the effects of uncorrelated instrumental and
calibration errors – see M21 for further discussion and an illustration.

The semivariogram is closely related to the two-point correlation function,
$\xi(h)$, defined to be the Pearson correlation coefficient for the value of
$Z$ at two points separated by a distance of $h$. Under the assumptions that
the random field under investigation is stationary and isotropic, the two
diagnostic functions are related by the following equation:

$\gamma(h)=\sigma^{2}\left(1-\xi(h)\right),$ (6)

where $\sigma^{2}$ is the overall variance of $Z$ throughout the field. Unlike
the two-point correlation function, the semivariogram is not normalised;
therefore, the semivariogram is sensitive to the amplitude of metallicity
fluctuations, whereas the two-point correlation function is not.

Another related metric often used in astronomy (and particularly cosmology)
for determining the covariance structure of a random field is the power
spectrum, which is calculated by taking the Fourier transform of a random
field and then determining the square of the magnitude of each frequency
component. By the Wiener–Khinchin theorem, for particularly well-behaved
random fields, this function captures the same details that the semivariogram
captures, but with units of spatial frequency rather than separation. However,
the formal requirements on the data quality to use Fourier methods to reveal
the covariance structure of the random field in this way are prohibitively
strict: the random field must have no edges and no missing data, with periodic
signals, and homoscedastic (constant variance) Gaussian error throughout
(Feigelson & Babu, 2012). For analysis of the cosmic microwave background
where these methods find their natural home, these requirements are readily
satisfied. However, metallicity maps of galaxies are always finite fields; are
heteroscedastic (not homoscedastic) due to galaxies being brighter (and
therefore having higher S/N) in their central regions; and have large portions
of missing data due to DIG contamination. We therefore recommend the