A forward-modelling approach to overcome PSF smearing and fit flexible models to the chemical structure of galaxies

Historically, metallicity profiles of galaxies have been modelled using a radially symmetric, two-parameter linear model, which reveals that most galaxies are more metal-rich in their central regions than their outskirts. However, this model is known to yield inaccurate results when the point-spread function (PSF) of a telescope is large. Furthermore, a radially symmetric model cannot capture asymmetric structures within a galaxy. In this work, we present an extension of the popular forward-modelling python package LENSTRONOMY, which allows the user to overcome both of these obstacles. We demonstrate the new features of this code base through two illustrative examples on simulated data. First, we show that through forward modelling, LENSTRONOMY is able to recover accurately the metallicity gradients of galaxies, even when the PSF is comparable to the size of a galaxy, as long as the data is observed with a sufficient number of pixels. Additionally, we demonstrate how LENSTRONOMY is able to fit irregular metallicity profiles to galaxies that are not well-described by a simple surface brightness profile. This opens up pathways for detailed investigations into the connections between morphology and chemical structure for galaxies at cosmological distances using the transformative capabilities of JWST. Our code is publicly available and open source, and can also be used to model spatial distributions of other galaxy properties that are traced by its surface brightness profile.


INTRODUCTION
Chemical enrichment of the interstellar medium is a powerful diagnostic of galaxy formation and evolution.The abundance of elements such as Carbon, Nitrogen, and Oxygen carries information about the star formation and merger history and about the accretion of pristine gas from the surrounding intergalactic medium (e.g.Maiolino & Mannucci 2019).
Studying spatially resolved chemical abundances as a function of cosmic time has long being hailed as one of the keys to understanding the processes that govern star formation (e.g.feedback and superwinds) as well as disentangling galactic growth mechanisms such as in-situ star formation, ex-situ star formation, major and minor mergers, cold-flow accretion, and environmental effects (Tremonti et al. 2004;Zahid et al. 2012Zahid et al. , 2014;;Gibson et al. 2013;Jones et al. 2013;Ho et al. 2015;Belfiore et al. 2019;Franchetto et al. 2021).
In the local Universe, much progress has been achieved in the past few years through the wholesale study of large samples of galaxies with integral field spectrography (IFS) (e.g.Sánchez et al. 2012; ★ E-mail: methab@student.unimelb.edu.auCroom et al. 2012;Bundy et al. 2015).A picture of the chemical structure of local spiral galaxies has emerged.In general, chemical abundance declines with radius, described by a so-called negative chemical gradient (e.g.Searle 1971;Vila-Costas & Edmunds 1992;Sánchez et al. 2014;Bresolin & Kennicutt 2015;Berg et al. 2013Berg et al. , 2020;;Sánchez-Menguiano et al. 2018;Poetrodjojo et al. 2021), which is interpreted as the result of inside-out star formation (the buildup of metals in the central regions over successive generations of star formation; e.g.Boissier & Prantzos 1999;Kobayashi & Nakasato 2011;Pilkington et al. 2012).
However, not all galaxies in the local Universe have metallicity profiles that follow this trend.Many studies of the abundance distributions of local galaxies have found that the metallicity does not decrease constantly throughout the entirety of a galaxy disc, showing flatter slopes in its inner regions (e.g.Belley & Roy 1992;Zinchenko et al. 2016), its outskirts (e.g.Martin & Roy 1995;Worthey et al. 2005;Bresolin et al. 2009;Bresolin 2011;Goddard et al. 2011;Lépine et al. 2011;Marino et al. 2012;Yong et al. 2012;Sánchez et al. 2014;Marino et al. 2016), or both (e.g.Rosales-Ortega et al. 2011;Sánchez et al. 2015;Sánchez-Menguiano et al. 2018).The most extreme deviations from this inside-out model show "inverted" gradients, in which metallicity increases with radius rather than decreasing (Sánchez et al. 2014;Pérez-Montero et al. 2016).
Modelling the chemical structure of galaxies at  ≳ 1 would allow us to deepen our understanding of how the interstellar media (ISM) of galaxies evolve, revealing how the processes that govern star formation have been changing over time.However, there are several problems that prevent the models that have been discussed thus far from being able to provide an insightful description of the chemical profile of galaxies at cosmic noon.
First of all, at these redshifts, spatially resolved metallicity measurements for galaxies are extremely difficult to obtain.Observational challenges include the faintness of distant galaxies, the need for infra-red wavelength coverage, and the limited angular resolution compared to the typical apparent size of galaxies.Having a poor angular resolution has been shown to flatten the inferred metallicity gradients of galaxies (e.g.Yuan et al. 2013;Mast et al. 2014;Acharyya et al. 2020) especially in their inner regions (Belfiore et al. 2017), which may lead to erroneous results if it is not corrected for.
Secondly, the models discussed thus far that deviate from a singlysloped negative metallicity gradient are not a priori physically motivated -rather, they are motivated by the fact that a simple negative linear gradient does not seem to fit the data.For this reason, many different theories have been proposed post-hoc to explain how the modelled radial profiles may have emerged.The presence of an inverted metallicity gradient may be interpreted as being due to strong tidal forces resulting from interactions (Kewley et al. 2006(Kewley et al. , 2010;;Rupke et al. 2010;Torrey et al. 2012), or from enriched gas being blown away by powerful winds driven by stupendous starbursts, while cold pristine gas is being accreted directly to the centre via cold-flows (Dekel et al. 2009).Flatter metallicity profiles in the outskirts of galaxies could be caused by a lack of sufficient material for star formation at large radii, preventing enrichment (e.g.Rosales-Ortega et al. 2011); galactic fountains wherein outflowing material is re-accreted onto the galaxy at large galactocentric distances (Belfiore et al. 2017); the accretion of metal poor satellite galaxies (Bird et al. 2012); radially symmetric gas inflow (e.g.Portinari & Chiosi 2000); or angular momentum transport (e.g.Lacey & Fall 1985), possibly driven by spiral arms, bars, or dynamical resonances (e.g.Martin & Roy 1995;Minchev et al. 2011).In the inner regions, proposed explanations for flattening of the abundance profile include accumulation of gas due to the tidal influence of bars or spiral arms (Rosales-Ortega et al. 2011;Sánchez et al. 2011) or at the inner Lindblad resonance (Fathi et al. 2007); but this could also simply be the effect of beam smearing in poorly-resolved data (Belfiore et al. 2017).Overall, the presence of breaks in a radial metallicity profile indicate that new processes beyond simple inside-out formation must be at play in the formation history of a galaxy, but they alone do not provide us with enough information to determine what those other processes are.
Thirdly, all of the models mentioned above assume that the metallicity profile of a galaxy is radially symmetric.There are many astrophysical reasons why this is not expected to be the case at high redshifts: (i) The frequency of minor and major mergers increases with redshift (Conselice & Arnold 2009).Mergers between galaxies of different metallicities ought to leave visible spatially correlated chemical signatures that will persist until the two galaxies are well-mixed.
(ii) Star forming galaxies at cosmic noon and beyond are generally not symmetric, with the fraction of clumpy and irregular galaxies increasing with redshift (Conselice 2014).The circularly-symmetric metallicity gradient model was developed to describe stable disc galaxies, is probably only apt to describe galaxies with Sa/Sb/Sc morphologies, and does not naturally apply to irregular systems.
(iii) Pristine gas is expected to be accreted along filaments, not in a spherically symmetric fashion (Kereš et al. 2005), causing asymmetric metal poor regions inside a galaxy.
(iv) Winds/outflows generally produce asymmetric regions of enriched gas (Martin et al. 2002;Cameron et al. 2021), and those gasrich regions are not expected to fall back in a spherically symmetric way.
(v) Spiral arms and bars are expected to produce azimuthal variations in the chemical structure of galaxies, either due to enhanced star formation along the spiral arms producing more metals (Spitoni et al. 2019) or from secular processes driving radial gas flows that bring metal poor gas in towards the galaxy centre and drive enriched gas out (Grand et al. 2016).
All of these astrophysical processes will leave observable signatures in the spatial distribution of metallicity within the galaxy that cannot be distinguished by a metallicity gradient (we demonstrate this in Figure 1).
However, given the observational challenges, most studies on the metallicity of the high-redshift Universe to date have focused on deriving integrated gas phase metallicities, (e.g.Sanders et al. 2020Sanders et al. , 2021Sanders et al. , 2023) ) or fitting simplified description of the chemical structure, such as a singly-sloped metallicity gradient model (e.g.Leethochawalit et al. 2016;Wang et al. 2022).Galaxies with strongly "inverted" gradients that do not agree with a simple inside-out formation model have recently been discovered at high redshift (Cresci et al. 2010;Queyrel et al. 2012;Swinbank et al. 2012;Stott et al. 2014;Troncoso et al. 2014;Wuyts et al. 2016;Wang et al. 2019) at a frequency that is much higher than what is generally predicted by state-of-the-art theoretical models, challenging current implementations of feedback (Ma et al. 2017;Hemler et al. 2021;Tissera et al. 2022).One reason may be because these azimuthally-asymmetric chemical features are not able to be detected by a metallicity gradient (Tissera et al. 2022).Indeed, recent simulations have shown that galaxies with very large positive or negative metallicity gradients show signs of morphological disturbances (Tissera et al. 2016).
Attempts have been made to find azimuthal variations in the abundance profiles of galaxies in the local Universe with a variety of statistics.Zinchenko et al. (2016) defined the "global azimuthal abundance asymmetry" as the difference between the average deviations from a singly-sloped metallicity profile between two half-planes of a galaxy, and used this statistic to search for large scale asymmetries in galaxies from the CALIFA DR2 catalogue.They found that this statistic was highly degenerate with uncertainties in the central position of the galaxy used for a radial fit, and hence could not find definitive evidence of any asymmetries in the data with this method.A related approach is to split galaxies into half-planes (e.g.Li et al. 2013) or quadrants (e.g.Rosales-Ortega et al. 2011), fit a linear profile to each sector, and search for statistically significant differences in the metallicity gradient between different regions.However, such models implicitly predict large discontinuities at a fixed angle of azimuth, which is not physically motivated.While they could potentially identify large scale azimuthal variations in the abundance profiles of galaxies, they have limited power in explaining where such variations come from.
Many studies have searched for evidence of azimuthal variation by measuring the scatter in the metallicity around a radially-symmetric abundance profile (e.g.Martin & Roy 1995;Rosolowsky & Simon 2008;Bresolin 2011;Croxall et al. 2015Croxall et al. , 2016;;Grasha et al. 2022, to name a few).In addition to not being able to explain the origin of any metallicity variations that are found, this method is confounded by the large uncertainties in metallicity diagnostics, especially present when strong emission line calibrations are used to determine the metallicity (Kewley & Ellison 2008).As such, metallicity variations found with this approach are treated as upper limits to any true azimuthally varying component.Metha et al. (2021) demonstrate how uncorrelated noise in metallicity estimates can be separated from small-scale azimuthal variations using the semivariogram, a mathematical tool adopted from geostatistics.However, this geostatistical approach requires high-resolution galaxy data, with sampling on scales of a few 100 pc (Metha et al. 2024 in prep.), making them unsuitable for the majority of high-redshift galaxy data currently available.
Another approach is to search for correlations between deviations from a radially symmetric metallicity profile and asymmetric features in a galaxy's surface brightness profile.This method has found some success in identifying asymmetric metallicity profiles in the local Universe and providing physically justified explanations for their origins.Several studies have found enrichment associated with the spiral arms of well-resolved galaxies (e.g.Martin & Roy 1995;Cedrés et al. 2012;Sánchez-Menguiano et al. 2016, 2018, 2020;Ho et al. 2017Ho et al. , 2018;;Vogt et al. 2017;Sakhibov et al. 2018).For some of these studies, these variations show an enrichment on the trailing edge of the spiral arm and a dilution on the leading edge, which implies that spiral arms drive large scale radial gas flows that bring metal-poor gas into the galaxy and drive metal-rich gas out, in agreement with zoomin galaxy simulations (Grand et al. 2015(Grand et al. , 2016)).Other studies show enrichment on the spiral arms, which can be attributed to enhanced star formation (Spitoni et al. 2019;Mollá et al. 2019).Ho et al. (2017) note a slightly different trend wherein spiral arms appear more metal rich with a rapid drop off in the metallicity at the leading edge and a gradual decline at the trailing edge, which they explain through a cycle of self-enrichment of Hii regions as they approach a spiral arm, followed by a rapid dilution after the density wave has been crossed.In all of these cases, morphological information from the light profiles of these galaxies is used to not only identify chemical substructures, but also to provide explanations on the origins of these features, providing a more comprehensive picture of galaxy evolution for these targets.
In this study, we present a method that enables physicallymotivated models to be fit for high redshift galaxies, wherein deviations from a one-dimensional metallicity profile are correlated with asymmetrical features in a galaxy's surface brightness profile.We implement this method in lenstronomy1 (Birrer & Amara 2018;Birrer et al. 2021), an open-source public general usage astronomy Python package that is widely used in the community.This extension of lenstronomy allows two-dimensional models to be fit to a galaxy's chemical structure in a general and rigorous way.It is rigorous in that it accounts for the effects of pixelisation and PSF smearing using lenstronomy's forward-modelling approach, allowing us to accurately recover metallicity gradients of distant galaxies at  ∼ 2, whereas traditional approaches would be biased; and it is general in that the code we present is entirely flexible, allowing the user to describe the distribution of metallicity with multiple components, as required by the data.This allows for large-scale deviations from a radially-symmetric trend to be captured, and allows these deviations to be associated to different morphological components of galaxies.
We organise this paper as follows.After a brief demonstration that motivates moving beyond a radial metallicity gradient (Section 2), in Section 3 we describe how we implement pixellisation and PSF smearing for light-weighted data in lenstronomy.In Section 4, we demonstrate the code by applying it to a simple gradient model, in the presence of PSF smearing and pixelisation.In Section 5 we illustrate the ability of the code to describe clumps in metallicity.Section 6 discusses the broader applications of the code.Section 7 provides a brief summary.

MOTIVATION
In this Section, we present a very simple illustration on how the assumption of circular symmetry can be misleading in understanding the chemical distributions of galaxies.In Figure 1, we show 2D metallicity profiles for four different simulated galaxies.Each of these profiles are meant to represent idealised versions of a different astrophysical scenario.Panel (a) shows a galaxy with a metallicity gradient that is efficiently mixed in the azimuthal direction.Panel (b) has a large asymmetry in its metallicity profile, as may be expected from the recent accretion of a metal-poor satellite.Panel (c) shows a galaxy with a chemical enrichment of 0.2 dex along its spiral arms, motivated by the observations of Sánchez-Menguiano et al. (2016); Ho et al. (2017), andVogt et al. (2017), reflecting a possible consequence of the enhanced star formation within spiral arms.In panel (d), we show a galaxy with correlated local fluctuations in its metallicity profile on small-scales, as found by recent studies of the local Universe (e.g.Metha et al. 2021Metha et al. , 2022;;Li et al. 2021Li et al. , 2023)), representing inefficient mixing of metals in the ISM.
Below their 2D profiles, we plot the radial metallicity profiles that would be recovered for each of these toy models, using two different common methods from the literature.In red, we show the best-fit relation between metallicity and radius computed using an ordinary least-squares approach, using every pixel as a data point.When this method is used, these four model galaxies have identical radial metallicity profiles.Although they are clearly distinct in two dimensions, reflecting the different physical processes at play, this information is lost when a one-dimensional model is fit.
In black, we show the median metallicity (with the 16 th and 84 th percentiles shown as errorbars) computed in 5 radial annuli of equal width for these galaxies.Grey dashed lines show the best-fit metallicity profiles to this azimuthally-averaged data.We find that for scenarios (b) and (c), where the metallicity profile is asymmetric on large scales, the metallicity profiles computed using all pixels and radially-averaged annuli do not agree.This disagreement is not a consequence of either method being less accurate than the other; rather, it is due to a fundamental failure of the radially-symmetric metallicity gradient model to describe the true chemical distributions of these galaxies.
From this demonstration, we can see how using a model that assumes circular symmetry could confound the interpretation of high redshift data, where metallicity maps are often clumpy and irregular (e.g.Cresci et al. 2010;Förster Schreiber et al. 2018;Wang et al. 2020;Curti et al. 2020b).For example, the accretion of a high metallicity satellite galaxy could mimic an inverted gradient, even though the formation pathways for these two scenarios are clearly distinct.

METHODOLOGY
Metallicity, as we observe it, is a light-weighted quantity.The metallicity measured within an aperture is not simply the average mass of elements heavier than Helium divided by the total mass of elements The radial metallicity profile of these four galaxies.Red solid lines show the best fit metallicity profile from an ordinary least-squares fit using all pixels.Despite the clear differences in their 2D profiles, this is exactly the same for each galaxy.Black points with error bars show the median, 16 th and 84 th percentiles of the metallicity in five annuli of equal width for each galaxy.Grey dashed lines are the best-fit metallicity profile using the median metallicity in these annuli instead.For cases (b) and (c), where there are large asymmetrical structures in the metallicity profiles of the galaxies, these two commonly-used methods do not recover the same radial metallicity profiles.contained within that region.Instead (assuming unbiased metallicity diagnostics), it is the average metallicity of the spectral energy density of all photons that fall within that aperture.
Let  ( ì ) and  ( ì ) be the intrinsic metallicity (reported in linear, not logarithmic units) 2 and surface profile of a galaxy, respectively.Then, the metallicity observed at each location Z (ì ) after being smeared by a point spread function with kernel  (ì ) is given by: To find the metallicity measured within an aperture , we must then take a light-weighted average of this quantity over a region: Here, Ĩ ( ì ) is the PSF-smeared surface brightness at each location ì  within the aperture .
Our new tracer module implemented in lenstronomy accounts 2 We note that these methods are not sensitive to the absolute value of the metallicity reported, and will work equally well in computing the effects of PSF smearing and blending between sources for any quantity that is linearly proportional to the metallicity.Because of this, our methods are not affected by the large systematic offsets between different metallicity calibrators that are known in the literature (Kewley & Ellison 2008).
for both of these effects, in order to properly integrate the effects of beam smearing and pixelisation of light-weighted quantities into a forward-modelling framework.Additionally, we also allow the light-weighted metallicity from different components to be blended together, to model galaxies that are not well described by a single light profile.For a model of a galaxy with multiple surface brightness profiles  1 (ì ), . . .,   ( ì ), each with their own unique metallicity profiles  1 (ì ), . . .,   ( ì ), we may compute the intrinsic metallicity profile produced by a sum of these components as: . (3) This metallicity profile can then be smeared by the point spread function (Equation 1) and pixelised (Equation 2) to model the metallicity that we would measure within each pixel.This process allows us to model families of metallicity distributions that are not radially symmetric, which may be more suitable for describing peculiar or merging systems.
We show examples of how these equations affect observed metallicity distributions of galaxies in Figure 2. Modelling the light distribution of all of these galaxies as an exponential, Ĩ ( ì ) ∝ exp(−||ì ||/  ), we show the effects of a light-weighted convolution of our four toy model galaxies introduced above with a Gaussian PSF kernel of two different widths.The first of these (middle row of Figure 2) has a PSF FWHM of ∼ 0.6  , and the second (bottom row) has a PSF FWHM of ∼ 1.3  corresponding to cases where the resolution is moderate or poor, respectively.We see that at moderate resolution, fine details on the small-scale structure of galaxies are buried, but other non-axisymmetric spatial trends associated with these galaxies can still be recovered.At poor resolution, spatial features such as spiral arms become washed out, but largescale variations in the metallicity profile of the galaxy can still be seen, such as the recent merger with a metal-poor companion shown in the second column.This demonstrates how (i) the effects of PSF smearing must be accounted for to properly model a galaxy's metallicity distribution (we make this argument quantitatively in Section 4), and (ii) even in the presence of large PSF smear, there are opportunities to fit multi-component 2-dimensional models to observed metallicity distributions in order to capture astrophysical information about the evolution history of these galaxies.
While this methodology was designed for resolved metallicity data, it can be equally well applied for any light-weighted quantity that can be spatially resolved over a galaxy, such as its velocity distribution, or the colour observed in each spaxel.This module was released as a part of lenstronomy package in Version 1.11.6, available through PyPI and the anaconda distribution.A tutorial on its use for fitting metallicity profiles is also available online. 3

USE CASE #1: UNBIASING METALLICITY GRADIENTS
COMPUTED FROM PSF-SMEARED DATA Acharyya et al. (2020) showed that when the size of a PSF becomes comparable to the scale radius of a galaxy, the metallicity gradients recovered using standard techniques were biased to recover flatter values for the metallicity gradient.In this Section, we demonstrate that the forward-modelling approach presented in this work and implemented inlenstronomy is able to account for the effect of PSF smearing and accurately recover the underlying metallicity gradients of galaxies, even when the effects of PSF smearing grow very large.
To demonstrate this, we simulated observations of a simulated galaxy with an exponential light profile, with an effective radius of   = 1", typical of a spiral galaxy (e.g.van der Wel et al. 2024).We did not account for any inclination effects, assuming this galaxy is observed face on.The metallicity profile of this galaxy was chosen to have a central metallicity of 12 + log([O/H]) = 8 and a strong metallicity gradient of −0.2 dex per arcsecond.
We simulated observations of this galaxy using Gaussian PSFs of 8 different widths, with full widths at half maxima (FWHMs) ∈ {0.1, 0.2, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0} arcsec.Using our light-weighed forward modelling approach implemented in lenstronomy, we simulated observations of this galaxy with a field of view of 8 ′′ × 8 ′′ , to capture the outskirts of this galaxy.The size of each pixel was taken to be 0.066 ′′ , to match the pixel size of JWST's NIRISS instrument.In this experiment, the size of a pixel was kept constant as the PSF was changed.
Uncertainty in the metallicity observed for each pixel is set by the signal to noise of the different emission lines in each pixel.We model this as following  2  =  1 +  2  , where  is the intensity of light in each pixel (see Appendix A for the physical motivation behind this noise model).We fit the constants  1 and  2 such that  2  = 0.01 at the centre of the galaxy, and  2  = 0.02 at a distance of 1  from the galaxy's centre.These values were chosen based on typical uncertainties found in resolved metallicity maps from gravitationally lensed observations (e.g.Wang et al. 2020).
3 github.com/astrobenji/lenstronomy-metals-notebooksUsing this collection of simulated galaxy observations, we then tested our ability to recover the galaxy's metallicity gradient in two different ways.First, we used a standard weighted least-squares (WLS) approach, ignoring the presence of a PSF and weighting the metallicity recovered in each spaxel by its inverse variance, as implemented through the Python package scikit-learn.In this procedure, the central metallicity,   , and the metallicity gradient, ∇, of a galaxy are estimated from a collection of metallicities {  } that are measured at physical distances {  } from the galaxy's centre.Error on measured values of {  } are treated as negligible, while each metallicity measurement has its own measurement uncertainty {  }.Then, best-fit values for the two parameters   and ∇ are computed using chi-squared minimisation: Secondly, we used our updated version of lenstronomy, incorporating the effects of PSF smearing to recover the metallicity gradients via a forward-modelling approach.Specifically, we used a particle swarm optimisation algorithm (Kennedy & Eberhart 1995) to find the initial conditions for a MCMC approach using emcee (Foreman-Mackey et al. 2013), then ran emcee using 100 walkers and 250 steps, discarding 50 steps as a burn-in stage.We compute the median value of ∇ and a 68% credible interval by looking at the distribution of samples fit by emcee, following the recommendations of the package's developers (Hogg & Foreman-Mackey 2018).
In Figure 3, we plot the best-fit values of ∇ computed for each simulated galaxy observation, with each of these two data pipelines.We see that, when a WLS approach is used with no PSF correction, the recovered metallicity gradient becomes flat when the with of the PSF becomes comparable to the radius of the galaxy (FWHM/  ≥ 0.7).This bias is significant, and the effect grows stronger when the beam is made wider.On the other hand, when lenstronomy is used to forward model the effects of PSF smearing, no such bias is seen, even when the FWHM of the PSF becomes much larger than the scale radius of the galaxy.Six out of eight MCMC trials contain the true value of ∇ within their 68% credible intervals, and all contain the true value of ∇ within a 95% credible interval (not shown), which is the statistically expected result.This indicates that a forward-modelling approach that accounts for the effects of PSF smearing is required to get unbiased results, such as the method implemented in lenstronomy.
While forward-modelling is a powerful approach for recovering spatial information from PSF-smeared data, its performance is still sensitive to resolution.Forward modelling cannot be used to recover any information on spatial scales that are smaller than a single pixel.Whereas the effects of PSF smearing can be modelled and corrected for, the data will not contain any information on processes that occur on spatial scales smaller than the scale at which the data is sampled (Cressie 1993).This represents a fundamental limit on the quality of data for which models of the metallicity distribution can be constructed.
To demonstrate this fundamental limitation, we explored how the metallicity gradients recovered using lenstronomy are affected by the size of the galaxy in question.For this exploration, we fix the pixel size to be half the FWHM of the PSF, simulating Nyquist sampling, and explore galaxies with effective radii that extend over 1, 1.4, 2, 3, 5, and 8 pixels, corresponding to beam widths of 2, 1.4, 1, 0.66, 0.4, and 0.25 FWHM/  .In each case, we set the field of view of each observation to be up to 4  from the galaxy's  1, with no PSF smear.Below these models, we show the light-weighted properties of these galaxies smeared by a Gaussian PSF with moderate (middle row) and poor resolution (bottom row).As the resolution becomes increasingly coarse, certain details of galaxy structure are lost.We show this effect quantitatively in Figure 3. Recovered values of the metallicity gradient as a function of the size of the point spread function, using the classical weighted least-squares approach without forward-modelling (blue), and using lenstronomy's implementation of a forward-modelling approach (red).When the size of the PSF, measured using its FWHM, becomes comparable to the size of the galaxy, the metallicity gradient recovered by the traditional approach becomes biased to values closer than zero.However, forward modelling can be used to eliminate this biasing effect.centre, to ensure the entirety of the galaxy lies within our field of view.
To assess whether any bias is seen in the recovered metallicity gradients using lenstronomy, we fit the metallicity profile using The effect of the galaxy size on the performance of metallicity gradient fitting using two methods: the standard approach that neglects PSF smearing (blue), and our forward-modelling approach (red).For marginallyresolved galaxies with ∼ 1 − 2 pixels per   , a forward-modelling approach yields more reliable metallicity gradient estimates.For small galaxies with   ≲ 1 pixels, there is no longer enough information in the data to precisely fit the metallicity gradient using either method.
each experimental setup 100 times, regenerating the noise from the simulated heteroskedastic metallicity variance maps for each trial.As a control, we repeat this using the WLS method without accounting for PSF smearing.We show the 16th, 50th (median) and 84th percentiles in the metallicities recovered at each resolution in Figure 4. We find that both methods become unreliable when the size of a pixel becomes comparable to the   of the galaxy, but lenstronomy produces results that are more consistent with the true metallicity gradient at all resolutions.In particular, when the galaxy extends over only a few pixels, the uncertainty on the metallicity gradient recovered by lenstronomy becomes large too, successfully capturing the lack of information in the data.
Such an effect of the uncertainty on the metallicity profile has been noted before, e.g. by Belfiore et al. (2017).For this reason, longitudinal studies that seek to capture radial metallicity profiles for a wide variety of galaxies in the local Universe are limited to galaxies of a certain size, set by the beam width of the telescopes in use.In Sánchez et al. (2021), in order to ensure the quality of their data, they restrict their analysis to galaxies with PSF FWHM/  < 0.5.Barrera-Ballesteros et al. ( 2023) take a similar approach, but use a different limiting galaxy size to ensure that PSF FWHM/  < 1.0.Our Figure 4 shows that, when a forward modelling approach is used to correct the metallicity profiles for the effects of PSF smear, a looser criterion of FWHM/  < 1.5 is required for metallicity gradients to be recovered accurately.This demonstrates that the use of the methods included in our newly updated version of lenstronomy hold the power to extend the analysis of local galaxies to a larger sample size, allowing the radial profiles of more compact systems that are more reminiscent of high-redshift galaxies to be explored.

USE CASE #2: FITTING BLENDED METALLICITY PROFILES TO CLUMPY GALAXIES
Most galaxies at cosmic noon are not disk galaxies (e.g.Conselice et al. 2008).Most galaxies at cosmic noon do not show a wellestablished metallicity gradient (Simons et al. 2021).To understand the chemical structure of these galaxies, more complex models of the chemical morphology of galaxies must be considered.
To demonstrate the ability to use lenstronomy to fit metallicity profiles other than a linear metallicity gradient, we construct a simulated setup of a clumpy galaxy, using three identical exponential light profiles with effective radii of   = 0.25 ′′ (shown in the left panel of Figure 5).The intrinsic metallicity of the central clump is chosen to be 12 + log(O/H) = 7.6 throughout.The second clump, offset by 0.25 ′′ in the positive x-direction from the central clump, has an intrinsic metallicity of 12 + log(O/H) = 7.8, and the final clump, above and to the left of the central clump separated by 0.35 ′′ is chosen to have a metallicity of 12 + log(O/H) = 7.9.Such a metallicity structure was chosen as a plausible setup involving three interacting dwarf galaxies coming together at cosmic noon, which would appear to have a positive metallicity gradient when observed using a radially-symmetric model.
We simulate an observation of this galaxy using the JWST instrument NIRISS, with a resolution of 0.066 ′′ per pixel, using a simulated PSF of the F150W filter computed using the Python package webbpsf 4 .As in the previous Section, we model the uncertainty in the metallicity in a heteroskedastic way, letting  2  = 0.01 at the centres of each clump, increasing to  2  = 0.02 at a distance of 1  from each clump.
We show the metallicity profile of this galaxy, without noise, in

7.9
Table 1.Median, 16th, and 84th percentiles of the recovered metallicities of each component of the simulated galaxy.We find that lenstronomy is very successful at accurately recovering the signal, even when the amount of noise is large.
the right-hand panel of Figure 5.We note that we have not simulated any mixing of metallicities between the three components that are interacting to form this system, as the turbulent processes that govern metal mixing in galaxies and their associated timescales are still not well understood.Instead, all gradients in metallicity between components are simply an artefact of the light blending between components and PSF smear.We attempt to fit the metallicity profile of this galaxy in three different ways.Firstly, we use the WLS method with no PSF correction to fit a radially symmetric metallicity profile to this galaxy.Secondly, we fit a linear gradient using lenstronomy, accounting for PSF smearing and pixelisation, but ignoring the structure of the galaxy.Finally, we use lenstronomy to attempt to estimate the intrinsic metallicity of each of the three components that makes up this galaxy.
We find that when a linear gradient is fit, without accounting for pixelisation or PSF smearing, a central metallicity of 12+log(O/H) = 7.73 ± 0.02 with an inverted metallicity gradient of 0.10 ± 0.03 dex/arcsec is preferred.This was consistent with the results found by lenstronomy when a linear gradient model was fit, with a central metallicity of 12 + log(O/H) = 7.73 ± 0.02 and a metallicity gradient of 0.12 ± 0.04 dex/arcsec being favoured.However, neither of these solutions are successful at describing the internal metallicity structure of this galaxy to high accuracy.In Figure 6, we plot the differences between the metallicities predicted by the linear gradient model fit using a least-squares approach and the true simulated metallicity of each spaxel.We find that there exist large, asymmetric residuals between the modelled and observed metallicity profiles with this galaxy, showing deviations of up to 0.1 dex, even within the brightest regions where the signal-to-noise ratio of the metallicity data was highest.
Finally, we use lenstronomy to attempt to recover the metallicity profile of this galaxy by modelling each clump as having a constant metallicity throughout.We find that, even with noisy metallicity data, if the light profile of the galaxy can be modelled accurately then the metallicity of each of the three clumps could be recovered with a high degree of accuracy, with errors smaller than 0.01 dex (see Table 1).The largest deviation between this model and the data was found to be only 0.0084 dex, representing a 10-fold improvement in prediction accuracy over models in which metallicity is assumed to be radially symmetric throughout a galaxy.More formally, using the Bayesian information criterion, the 3-component model is highly preferred to a linear gradient model for this system, with ΔBIC= 35.

DISCUSSION
Our approach of using forward-modelling to correct for the effect that PSF-smearing has on flattening metallicity gradients is not unique.A similar approach has been discussed in Carton et al. (2017), and applied to MUSE data in Carton et al. (2018).However, this forwardmodelling approach differs from our own in several ways.Firstly, their model incorporates a photoionisation model, and so requires as input Figure 6.The difference between the metallicity of the galaxy measured within each pixel, and the metallicities predicted using a linear-gradient model.Because this galaxy is not radially symmetric, the metallicity gradient model leaves large residuals of the order 0.1 dex between the observed and recovered metallicities.several emission-line images.On the other hand, the lenstronomy approach takes as input a metallicity map that can be constructed a priori through any method that the observer has access to, making it more flexible.Secondly, the model of Carton et al. (2017) allows resolved metallicity measurements to be binned (i.e.into a Voronoi tessellation), whereas our code does not have this feature.Finally, the model of Carton et al. (2017) assumes that the true metallicity distribution of a galaxy is always a metallicity gradient, while our code allows the user to specify any functional form for the underlying metallicity structure and surface brightness distributions, allowing for asymmetries in a galaxy's chemical profile to be captured.Furthermore, lenstronomy also contains the capability to work on gravitationally-lensed sources, allowing the intrinsic metallicity distributions of such systems to be uncovered the source plane.
The quality of our data is improving, and the quality of our models ought to improve with it.Thanks to JWST, we are already collecting IFS data of gravitationally-lensed galaxies at  ∼ 3 with sub-kpc resolution (Wang et al. 2022).Such exquisite data quality allows models that include local variations in metallicity that depart from radial symmetry to be fit.These non-axisymmetric models allow signatures of gas-driven galaxy evolution, such as minor mergers, recent starbursts, asymmetric outflows, or cold flows of pristine gas into a galaxy's outskirts along cosmic filaments to be recognised in observational data.Tissera et al. (2016) examined abundance profiles for galaxies within a cosmological simulation, and found that the galaxies with the steepest positive and negative metallicity gradients showed evidence of substructures such as rings, bars, and close companions, demonstrating a link between abundance profiles and morphology.The techniques that we explore in this work allow us to delve deeper into this relationship in a quantitative way.
Such irregular metallicity profiles have been seen in observational data for high redshift galaxies.Förster Schreiber et al. ( 2018) examine the two-dimensional distributions of the line ratio [N ii]/H, which is approximately proportional to metallicity (Pettini & Pagel 2004).They find that often, large asymmetries can be seen in the [N ii]/H profiles, reflecting asymmetries in the surface brightness profiles of these systems.A similar effect is commented on in Curti et al. (2020b).In their Figure 11, they show four galaxies observed at high redshift with large asymmetries in their metallicity profiles, revealing clumpy substructures which could be a sign of ongoing interactions within a merging system, or different phenomena such as gas flows that act on timescales shorter than secular processes.
When a radial gradient alone is used to characterise all metallicity variation within high-redshift galaxies, radially-symmetric physical models are employed to explain this variation.Of the three highredshift galaxies exhibiting inverted metallicity gradients presented in Cresci et al. (2010), two show visible departures from radial symmetry in their 2D metallicity profiles (their Figure 1).However, as this variation was not captured, all three radial metallicity profiles were explained using the same model of cold gas flows being preferentially directed into a galaxy's centre (Kereš et al. 2005).Schönrich & McMillan (2017) argue that this model may not be sufficient to explain the inverted gradients captured within this data, as the effects of enhanced star formation efficiency in a galaxy's central regions may increase the metallicity of the central regions more than could be compensated for by a pristine gas flow.Fitting more flexible models to the data has the potential to reveal alternative explanations for such inverted, irregular metallicity profiles.
The systems with large asymmetries in their chemical profiles reported in Curti et al. (2020b) also show large deviations from radial symmetry in their photometric and kinematic data.This new version of lenstronomy allows us to account for information present in these ancillary data product, to inform us of the presence of substructures within the galaxy when fitting more complex metallicity profiles.This gives us a more holistic approach to describe the chemical structure of galaxies, combining data from different instruments and surveys to produce a more complete description of the dominant physical processes that govern galaxy evolution at high redshift.

Potential applications to other galaxy properties
Many different publically-available code bases exist that use forward modelling to accurately fit light profiles to galaxy data, including the effects of pixelisation and PSF smearing, such as GALFIT (Peng et al. 2002) and galight (Ding et al. 2022).This work presents an extension of lenstronomy that uses forward-modelling to compute the effects of gravitational lensing, PSF smearing, and pixelisation on the spatial distribution of light-weighted properties of a galaxy.This paper investigates the performance of this code base on accurately resolving the metallicity distributions of galaxies.However, the algorithms we have developed will also work on any other lightweighted property of galaxies.For example, forward-modelling of elemental abundance ratios such as [C/O] would allow us to see signatures of recent star formation, as  elements such as Oxygen are primarily released in core-collapse supernovae with delay times of ∼ 10 Myr, whereas a significant amount of Carbon is additionally released during the end of life of intermediate mass stars, with delay times closer to ∼ 100 Myr−10 Gyr (Kobayashi & Nakasato 2011;Kobayashi et al. 2020;Jones et al. 2023).However, such observations would be complicated by details of turbulent metal-mixing in these high-redshift galaxies, which is not well understood (Metha et al. 2021).
Alternatively, the age distribution of stars throughout a galaxy could be constrained via SED fitting.Such an analysis has already been attempted by Giménez-Arteaga et al. (2023).By fitting a SED independently to each pixel of a galaxy observed by NIRCam, several properties, including the age, dust attenuation, and star formation rate density of the galaxy could be recovered for each pixel.However, the effects of PSF smearing through pixels of different light sensitivity were not accounted for in this analysis, limiting their interpretation of resolved variables to scales larger than a PSF beam.Using lenstronomy, these final age maps could be combined with light profile data in order to recover the intrinsic age distribution of stars within these galaxies at native resolution, perhaps allowing a detailed analysis on the properties of small star-forming clumps to be performed on high redshift galaxies.
Another natural application of this forward modelling approach would be to analyse kinematic structures in galaxies captured by IFS data.Like metallicity, velocity dispersion is a light-weighted quantity.A similar approach has already been used by de Graaff et al. (2023) to uncover kinematic profiles of six galaxies at 5.5 <  < 7.6 observed with JWST NIRSpec data, using a bespoke forward-modelling pipeline that is specifically designed to account for the effects of the microshutter arrays used for multi-object spectroscopy (MOS).The new module implemented in lenstronomy is less suitable for NIR-Spec MOS data than the code that they present, but would be suitable for observation campaigns using more typical IFS instruments, such as NIRISS or MUSE.

SUMMARY
Spatially-resolved observations of light-weighted properties such as metallicity can be heavily impacted by instrumental effects, such as PSF smearing and pixelisation.In this work, we present an added module in the Python package lenstronomy, and demonstrate its ability to correct for these effects, as well as fit flexible models of the chemical structure of galaxies, through forward-modelling.We summarise our main conclusions below: • We explored the ability of lenstronomy to recover the metallicity profile of a galaxy under the effects of PSF smearing.We found that, when the size of a pixel was small compared to the size of a galaxy, the new lenstronomy code was able to recover the true metallicity gradients of a simulated galaxy (observed with realistic noise) without bias, even when the size of the PSF became comparable to the size of the galaxy itself.This makes this code very useful for situations where an instrument has a large PSF that extends over many pixels.
• We additionally tested the performance of lenstronomy over Nyquist-sampled mock-IFS data.We found that accounting for PSF smear and pixelisation using forward-modelling yielded more accurate metallicity gradients than when these effects were ignored for marginally resolved data (∼ 1.5 − 3 pixels per   ), allowing accurate metallicity profiles to be recovered for smaller systems.Data with a resolution of ∼ 1 pixel per   was found to not contain enough information for a metallicity gradient to be fit precisely even when PSF smear is accounted for, leaving large uncertainties on the recovered metallicity profile parameters.
• To test the ability of lenstronomy to fit non-linear metallicity profiles of galaxies, we simulated a galaxy composed of three "clumps" of identical brightness and different metallicities.We found that, while traditional model-fitting pipelines would characterise this galaxy as having an inverted metallicity gradient, lenstronomywas able to successfully resolve the metallicities of the three components with high accuracy.
These tests demonstrate the utility of lenstronomy to capture accurately details of a galaxy's chemical profile in two-dimensions, even when the seeing is poor.This is particularly relevant at highredshift, where a multitude of processes are expected to cause large deviations from a circularly-symmetric galaxy profile.In future works, we will apply these techniques to resolved maps of gravitationally-lensed galaxies captured by NIRISS in the GLASS-JWST survey (PI:Treu).As with the rest of the lenstronomy package, all of our code is actively maintained and documented, is available for public use under the MIT license and is distributed through the python packaging index.

Figure 1 .
Figure1.Top: The four different simulated metallicity profiles, with no noise, representing four different astronomical scenarios: (a) a smooth profile with a negative metallicity gradient, representing inside-out star formation; (b) an asymmetric metallicity profile, representing a galaxy that has recently undergone a major merger with a metal-poor companion; (c) a galaxy with metal-enrichment along its spiral arms; (d) a galaxy with small-scale correlated metallicity fluctuations, representing inefficient metal mixing.Bottom: The radial metallicity profile of these four galaxies.Red solid lines show the best fit metallicity profile from an ordinary least-squares fit using all pixels.Despite the clear differences in their 2D profiles, this is exactly the same for each galaxy.Black points with error bars show the median, 16 th and 84 th percentiles of the metallicity in five annuli of equal width for each galaxy.Grey dashed lines are the best-fit metallicity profile using the median metallicity in these annuli instead.For cases (b) and (c), where there are large asymmetrical structures in the metallicity profiles of the galaxies, these two commonly-used methods do not recover the same radial metallicity profiles.

Figure 2 .
Figure 2.In the top row, we show the four toy metallicity distributions shown in Figure1, with no PSF smear.Below these models, we show the light-weighted properties of these galaxies smeared by a Gaussian PSF with moderate (middle row) and poor resolution (bottom row).As the resolution becomes increasingly coarse, certain details of galaxy structure are lost.We show this effect quantitatively in Figure3.

Figure 3 .
Figure 3. Recovered values of the metallicity gradient as a function of the size of the point spread function, using the classical weighted least-squares approach without forward-modelling (blue), and using lenstronomy's implementation of a forward-modelling approach (red).When the size of the PSF, measured using its FWHM, becomes comparable to the size of the galaxy, the metallicity gradient recovered by the traditional approach becomes biased to values closer than zero.However, forward modelling can be used to eliminate this biasing effect.

Figure 4 .
Figure 4.The effect of the galaxy size on the performance of metallicity gradient fitting using two methods: the standard approach that neglects PSF smearing (blue), and our forward-modelling approach (red).For marginallyresolved galaxies with ∼ 1 − 2 pixels per   , a forward-modelling approach yields more reliable metallicity gradient estimates.For small galaxies with   ≲ 1 pixels, there is no longer enough information in the data to precisely fit the metallicity gradient using either method.

Figure 5 .
Figure5.The light profile (left) and metallicity profile (right) for a simulated galaxy with an irregular morphology and a non-linear metallicity profile.This galaxy is made of three different components, all of which are modelled to have exponential light profiles with the same brightness and the same   .Each of the three components of this galaxy is modelled as having a constant, but different, metallicity, with the lowest metallicity in the centre, and the highest metallicity in the top-left corner.