Measuring the dust attenuation law of galaxies using photometric data

Fitting model spectral energy distributions (SED) to galaxy photometric data is a widely used method to recover galaxy parameters from galaxy surveys. However, the parameter space used to describe galaxies is wide and interdependent, and distinctions between real and spurious correlations that are found between these parameters can be difficult to discern. In this work, we use the SED fitting code BAGPIPES to investigate degeneracies between galaxy parameters and the effect of the choice of different sets of photometric bands. In particular, we focus on optical to infrared wavelength coverage, and on two parameters describing the galaxies' dust attenuation law: $A_V$ and $\delta$, which characterize dust column density and the slope of a flexible dust attenuation law, respectively. We demonstrate that 1) a degeneracy between the residual (the difference between truth and recovered value) $A_V$ and star formation rate exists, but this is lifted when WISE bands are included; 2) BAGPIPES is able to accurately recover the input $A_V$ and $\delta$ distributions and relations (differences in slope of less than 1.7$\sigma$ for a flat relation, less than 1.2$\sigma$ for an observationally-motivated relation from Salim et al. 2018) and is not introducing spurious correlations between these parameters. Our findings suggest that the information needed to constrain $A_V$ and $\delta$ well enough individually exists in the data, especially when IR is added. This indicates that recent works finding a correlation between $A_V$ and $\delta$ are not being misled by fitting degeneracies from their SED fitting code.


INTRODUCTION
Galaxies are incredibly complicated systems, where their defining properties are highly dependent on their evolution and environment.In addition, the information needed to describe a galaxy is quite varied, typically a model needs to account for their stellar mass, rate of star formation, dust content, metallicity, star formation history, and redshift (see (Conroy 2013) for a review).All of this information about the galaxy is compressed into the observational signal we observe, its electromagnetic radiation in a form of a Spectral Energy Distribution (SED).Comparing the observed SED to synthetic models is an extremely widespread way of deriving the physical properties of galaxies.
Several works have explored the parameter degeneracies and systematics associated to SED fitting.For instance, studies have revealed the inherent relationships between choice of star formation history (SFH) parameterization (or lack thereof, non-parameteric models have proven promising, see ⋆ E-mail: palmese@cmu.edu† NASA Einstein Fellow Leja et al. 2019;Lower et al. 2020) and recovered stellar mass (Lower et al. 2020) or age (Simha et al. 2014), treatement of metallicity and stellar mass (Mitchell et al. 2013), and the degeneracy between age and metallicity (Worthey 1994).Other studies have investigated the effects of galaxy dynamics on SED fitting, considering phenomena such as quenching (Ciesla et al. 2016), morphology (Wuyts et al. 2009) or past mergers (Zine & Salim 2022).Systematic uncertainties have also been shown to enter results as a function of redshift when attempting to determine stellar masses (van der Wel et al. 2006;Paulino-Afonso et al. 2022).Despite this, little exploration has been performed with a varying dust attenuation law, though studies have been done which found that dust's effect on other galaxy parameters is non-negligible (Lo Faro et al. 2017;Leja et al. 2018;Lower et al. 2022).Due to the attenuation law parameters introducing extra degeneracies with other physical properties and with each other (Qin et al. 2022), most analyses tend to assume a fixed slope for it.
The assumption of a fixed slope is not the case for recent works that have attempted to recover the dust law for Supernovae (SN) host galaxies (Meldorf et al. 2022;Dixon et al. 2022).These works in particular sought to resolve previously unexplained correlations between the host galaxy mass and Hubble diagram residuals, called the mass step.As physical parameters such as dust and stellar mass can or are already used through the aforementioned correlations to improve the standardization of supernovae (Sullivan et al. 2011;Betoule et al. 2014), a biased estimate of these parameters could lead to systematics in the resulting cosmological measurements.For example, due to the dimming caused by dust, failing to account for dust can cause an overestimation of luminosity distances to these SN.Hence, systematic errors may be introduced into measurements of parameters calculated from the standard candle relationship of Type Ia SNe, such as the dark energy equation of state parameter, w, or the matter density of the universe Ωm, if galaxy properties are not properly estimated (Paulino-Afonso et al. 2022).
In this work, we focus specifically on the dust content of galaxies.The effect of dust on the light we observe is twofold.The first effect is extinction: light from the galaxy is absorbed and reddened through interactions with dust particles.The second effect is attenuation.Attenuation is a broader phenomenon which includes extinction effects, but also accounts for the redirection of light by scattering off dust.This has the dual effect of reflecting light that would have missed the observer had there been no dust into the path of the observer and vice versa.We describe these phenomena with attenuation laws or curves.These curves represent the ratio between luminosity emitted (no dust interaction) and received (affected by dust) as a function of wavelength.
While this attenuation curve is conveniently specified by a small number of parameters, fitting these parameters accurately can be challenging, as the dust parameters AV and δ tend to be degenerate with other properties such as star formation rate, and with each other by definition.Many recent works have found that a higher value for AV correlates with a flatter attenuation curve, i.e. a higher value of δ (Arnouts, S. et al. 2013;Kriek & Conroy 2013;Salmon et al. 2016;Leja et al. 2017;Salim et al. 2018;Decleir et al. 2019;Battisti et al. 2020;Boquien et al. 2022).However, some recent works (Qin et al. 2022) have claimed that said AV − δ relationship could be driven by a degeneracy between the fitting parameters.We attempt to reproduce the results of Qin et al. (2022) to determine if the AV − δ correlation is driven by fitting degeneracies.
In this work, we seek to explore the effect that the SED fitting process has on recovering specific parameters.We first test the reliability of SED fitting for parameters of interest by considering the distributions of residuals between truth values (the values the models were generated with) and recovered estimations of these parameters.We then specifically analyze the distributions of dust parameter residuals, searching for correlations between the main dust parameters AV and δ, other parameters, and their respective errors or bias distribution scatter.We finally turn to the analysis performed in Qin et al. (2022) and attempt to reproduce their results.We extend their analysis to a data-driven AV − δ input distribution.Throughout this paper, we will note how different combinations of input bands used for data affect the results.In particular, we compare different combinations of the Dark Energy Camera (DECam, Flaugher et al. 2015a) ugriz bands, Visible and Infrared Survey Telescope for Astronomy (VISTA; Emerson et al. 2004) JHKs bands, and Wide-Field Infrared Survey Explorer (WISE, Wright et al. ( 2010)) bands.
We only focus on photometric data, although assume that the redshift is known from spectroscopy (similarly to the analysis in Meldorf et al. 2022).This paper is organized as follows.In §2 we describe our SED-fitting code and the model galaxies used in our analyses.In §3 we present our results and discuss, and conclude in §4.All error bars are 1σ unless otherwise stated.

METHODS
Our model SEDs are both created and fitted using the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (BAGPIPES; Carnall et al. 2018) software.BAGPIPES is a fully Bayesian spectral-fitting code used to estimate galaxy properties from photometric and spectroscopic data.BAGPIPES also allows us to reverse this process: estimating the SED for a galaxy with user-given parameters and the measured fluxes in specific photometric bands of the theoretical SED.BAGPIPES allows us to analyze the effect sophisticated dust attenuation curves as well as several parametric SFHs have on measured SEDs.BAGPIPES utilizes the stellar population models derived in Bruzual & Charlot (2003) and relies on the MultiNest nested sampling algorithm (Feroz & Hobson (2008);Feroz et al. (2009Feroz et al. ( ), 2013)), specifically implemented through the PyMultiNest interface (Buchner et al. 2014), to obtain the posterior distributions for desired parameters based on SFH and dust models, prior distributions and observational data provided by the user.

Attenuation Curves
Parametric models of attenuation laws come in a variety of types; this work considers the Noll et al. (2009) modification of the Calzetti law as formulated by Salim et al. (2018).This flexible attenuation law can take a variety of shapes and is governed by three parameters: AV , δ and B. AV is the dust attenuation in the V band (at λ ≃ 5500 Å).It can be considered as a normalization constant for the attenuation curve, and is proportional to dust column density along the line of sight from object to observer.A higher AV indicates a dustier environment, and will result in a higher attenuation curve at all wavelengths.δ is a parameter that governs the difference between the Calzetti and Salim attenuation curves, at δ = 0 they are equivalent.In addition, δ is directly related to the parameter RV via the relationship: where R V,Cal = 4.05.The final parameter, B, represents the strength of a "bump" in total attenuation which peaks at a wavelength of 2175 Å.In previous works (Meldorf et al. 2022) and in tests described below, we find that letting this parameter vary freely has minimal effect on the recovery of physical parameters, and we therefore choose to set it to 0 for most of this analysis.

Simulations
Each model galaxy is given a unique combination of input parameters that govern its dust attenuation curve, SFH and redshift.We simulate galaxies in two different ways, which are used in different parts of the analysis.The first simulation set (hereinafter Sim1) is made in a grid of parameter values, while the second set (hereinafter Sim2) follows a continuous, realistic distribution of parameters.For Sim1, each parameter and its values are given in Table 1.The dust parameters considered are those given in Section §2.1, and truth values are selected to be evenly spaced between the minimum and maximum value we consider.This lets us explore the effect of the SED fitting as a function of the input parameters more broadly than we would with a more realistic distribution of galaxies where the bulk of the objects then to cluster around specific parameters values.Hence, AV has input values of 0.100, 0.525, 0.950, 1.375, and 1.800, and δ is given −1.400, −0.975, −0.550, −0.125, 0.300.The SFH of each galaxy is parameterized as a lognormal function, which was first utilised in Gladders et al. (2013) and revised in Simha et al. (2014): where t is the age of the universe, SFR is the star formation rate, and τ l and T 0,l are free parameters.We follow Diemer et al. (2017) and Carnall et al. (2019) in redefining these parameters in terms of the more intuitive tmax (time at which SFR is maximal) and σSFH (width of the SFH peak at half maximum height) as follows: and We again pick our input values for these parameters to be uniformly spaced throughout the range of possible values, though for these SFH parameters we additionally choose to avoid including values that would be significantly longer than the age of the universe.Thus, tmax has input values 2, 4.66, 7.33, and 10 Gyrs, while σSFR is given 4, 7.33, 10.66, 14 Gyrs as input values.The redshift values chosen are 0.2, 0.4, 0.6, 0.8, so as to match most of the galaxies with spectroscopic redshifts in e.g. the Dark Energy Survey Supernova (DES SN) program.For the total mass formed, we again pick evenly spaced values within the prior range, though for this parameter they are spaced logarithmically, giving log10(M form /M⊙) = 8.5, 9.5, 10.5, and 11.5.The initial mass function used by BAGPIPES is given in Kroupa & Boily (2002).We similarly choose values for the metallicity of the galaxy, Z, to be geometrically uniformly spaced within the prior, specifically 0.5, 1.0, and 2.0 times the solar metallicity (Z⊙).
In order to create our models, we consider every possible combination of all parameters, leading to a total number of models that is the product of the number of values considered for each parameter.For the values in Table 1, this results in 19,200 unique models in total.
Entering these parameters into BAGPIPES yields values for photometric measurements of each model galaxy, specifically we produce photometric data in ugriz bands of the Dark Energy Camera (DECam; Flaugher et al. 2015b), and the Visible and Infrared Survey Telescope for Astronomy (VISTA; Emerson et al. 2004) JHKs bands.We choose to include the same bands used in the DES deep fields, as in the DES SN host studies of Meldorf et al. (2022), and also include the WISE 1,2,3 and 4 bands, since they are available over the entire sky, and specifically in combination with grz optical bands in the Dark Energy Spectroscopic Instrument (DESI) Legacy Survey (Dey et al. 2019).
In order to simulate measurement noise, we assume a signal-to-noise ratio of ten and assign to each data point an error of its flux value divided by ten.This is a conservative choice if we want to use our findings for the DES deep fields or Rubin LSST, where the optical bands reach a greater signal to noise.We then pick a random value drawn from a normal distribution centered at the true flux with a standard deviation of the error to each photometric data point.A similar process is applied to the redshift measurements, where we take the error to be 0.001.
Whereas in Sim1 we allow for every possible combination of every parameter, in Sim2 we specifically select for realistic combinations that could be seen in galaxies.We still use Sim1, although it is not as realistic as Sim2, in order to better explore a wider dynamical range of parameters.To generate Sim2, we first follow the same method as outlined above with the following new additional steps.First, we develop our treatment of SNR to be more realistic.Rather than setting one value of SNR for all galaxies in all bands, we use measured relationships between SNR and band magnitude in real galaxies to determine SNR for each data point.Using the data used in Meldorf et al. (2022) as well as for spatially matched galaxies in the Legacy Survey data, a linear fit between magnitude and SNR was calculated for every band in this work.Then, in the model creation process, each data point for each galaxy is assigned a SNR determined using its magnitude and the derived linear relationship.An error is then drawn randomly from a Gaussian distribution centered at 0 with standard deviation of flux divided by this SNR.This gives this sample of models a distribution of SNR in each band akin to a real galaxy dataset.
In addition, we replace our discrete AV distribution with a continuous one.Instead of looping through five values of AV as we create our models, we select five random values of AV from a uniform distribution ranging from 0.3 to 1.7 for each parameter combination.To clarify, whereas in our Sim1 sample there were five copies of every galaxy which only differed by their value of AV , there are still five copies of every galaxy which only vary in their AV value but now these values are randomly selected rather than predetermined.δ is now also a continuous parameter, randomly drawn between -1.6 and 0.4.We distinguish between Sim2, where both δ and AV take a range of independent values, Sim2a, where δ is fixed, and Sim2b, where δ is calculated from AV based on the Salim & Narayanan (2020) empirical relation described in more detail below.After models are generated we apply a series of cuts to mimic a realistic dataset.We first select only galaxies with a sSFR (the specific star formation rate, which is SFR divided by the stellar mass) in the range 10 −13 < sSFR < 2 × 10 −9 as these are typically expected values (see, e.g., Ilbert et al. 2015).We then also enforce that the magnitude in the i and r bands be consistent with observed values.We cut any galaxies in our sample that do not fall between the 1st and 99th percentiles for the i and r bands in the Meldorf et al. (2022) data.Due to the now random nature of the model creation process, a slightly different number of models passes the selection criteria for each sample, and our Sim2 sample size ranges from 2,883 to 2,945 galaxies for the different cases outlined in §3.2.Moreover, we make simulations that include the 2175 Å bump.We use the Sim2 model setup but draw triple the number of galaxies with each galaxy given a randomly drawn bump strength value in the range (0.15, 0.85); this ensures that for each combination of all other parameters we have three unique bump strength values.
BAGPIPES allows for treatment of dust's re-emission of absorbed light as parameterized by Draine & Li (2007).Their model uses three parameters, qPAH, Umin, and γ.PAH stands for polycyclic aromatic hydrocarbon, and qPAH is the mass fraction of the dust in the form of these molecules.Umin is the minimum of the distribution of starlight intensity that the dust is exposed to, and γ is the mass fraction of the dust that is heated by starlight with an intensity above the minimum Umin (Draine & Li 2007).Umin is defined as a scaling factor of the intensity of the interstellar radiation field as estimated by Mathis et al. (1983) for the solar neighborhood, meaning it is unitless like γ and qPAH.
To analyze the residual distributions of these parameters, we again use the Sim2 model setup using all photometric bands but draw triple the number of galaxies, where each galaxy triplet has the same combination of all other parameters, but a unique dust emission value.Note that if we were to take this approach for each dust emission parameter, we would end up with a factor of 27 more models than in the default Sim2 case, which is computationally expensive to the point of being prohibitive.Thus, we instead draw three random sets of the emission parameters, rather than drawing three of each parameter and considering every possible combination.We assume input values ranging over the values taken from Draine & Li (2007) (qPAH from 0.1 to 4.6, Umin from 0.1 to 24, γ from 0.0005 to 1).Our final set of Sim2 simulation includes 2 different subsets of simulations: one with no dust emission, one with dust emission.
To conclude, our landscape of Sim2 simulations includes: • five different band combinations • three different δ − AV relations (unconstrained, slope 0, and empirical slope) • two bump schemes (no bump, bump) • two dust emission schemes (no emission, variable Draine & Li parameters).

SED fitting
The simulated photometric data of these galaxies are then fed into BAGPIPES so that we can extract galaxy parameters as though they were real data.We assume the same Salim dust law and lognormal SFH as was used to generate the models in order to directly determine how accurately specific parameters are recovered.Prior distributions are chosen such that the prior encapsulates the entire range of possible generated values.Prior distributions on dust parameters are all uniform, ranging from 0 to 2 and -1.6 to 0.4 for AV and δ respectively.The SFH priors are also uniform, where the prior on tmax extends from 0.1 to 15 Gyr, and the full-width-half-maximum prior covers 0.1 to 20 Gyr.The prior on metallicity allows for values between 0 and 2.5, and is again uniform.The prior on the mass formed ranges from 10 1 to 10 15 M⊙ and is uniform in logarithmic space.Finally, the prior for the redshift is a gaussian with a standard deviation of 0.001 centered at the "measured" value of redshift, i.e. the truth redshift with a simulated scatter added on.All of these priors are also given in Table 1.For the dust emission, we assume the prior range that follows the input for the simulation (qPAH from 0.1 to 4.6, Umin from 0.1 to 24, γ from 0.0005 to 1), and give each parameter a logarithmic prior in that range following Chastenet et al. ( 2019).

RESULTS
An example of a corner plot of resultant posteriors for one galaxy is shown in Figure 1.This galaxy was selected due to its residuals (truth value minus recovered value) being close to the average residual for each parameter of interest.The entered truth values for this galaxy were AV = 0.525, δ = −0.975,σSFR = 14.0, log 10 (M⋆/M⊙) = 10.27,Z/Z⊙ = 1.0, tmax = 10.0, z = 0.601.Observing the posteriors for each parameter will give the reader an idea of the degeneracies and uncertainties present in the SED fitting process.For instance, the broad joint 2D posteriors between metallicity and tmax or σSFR reflects the difficulty in fitting galaxy age and metallicity (Worthey 1994).The stellar mass posterior is strongly correlated with the SFH parameters (tmax or σSFR) as the stellar mass is calculated from the SFH.Finally, one can see that there does exist a strong degeneracy between AV and δ as is to be expected.It remains to seem whether this fitting degeneracy is causing false correlations to emerge in the fitted data.A somewhat milder degeneracy is also present between stellar mass and the dust parameter.
Note that unless otherwise specified, in the following dust emission is set to 0 and not fit for.This is because only the longest wavelength band considered here, W4, is signif-icantly affected by dust emission and does not significantly contribute to the dust absorption measurements.In other words, including or excluding W4 in the following where we have turned off dust emission does not have an impact on our results, as it is confirmed by our findings below.For similar reasons, we initially use the simulations with no bump.We separately assess the impact of variable bump strength and dust emission parameters on our results in Sec.3.3 and Appendix A. We now focus on the residual distributions of the parameters defined as the difference between the truth and the recovered values.The distributions of the residuals for each of the seven parameters is shown in Figure 2 for galaxies from the Sim1 data set, using the ugrizJHKs bands.All of the residual distributions follow roughly Gaussian distributions with the exception of the SFH full width half maximum (σSFH), which has four discrete peaks.This seems to be due to the fact that recovered σSFH is uncorrelated with the true σSFH.
For the majority of models, regardless of their σSFH truth value, BAGPIPES tends to yield a σSFH near a value of 13.5 Gyr, corresponding to a very flat SFH, with large uncertainties, averaging to 4.8 Gyr.This is a consequence of our inability to reconstruct the entire SFH from photometric data, which is not surprising (Carnall et al. 2019;Leja et al. 2019).While only 14% of galaxy models were given input σSFH values within 1 Gyr of 13.5, 75% were calculated to have σSFH values in this range.This tendency for all galaxies to recover the same σSFH results in the multi-peaked distribution we recover; since the recovered value is almost constant and the truth values are discrete, the residual distribution becomes almost discrete.Because estimating this parameter is meaningless in this analysis, we ignore it in the following.We will also not focus on any SFH parameters, except for the SFR, which is the only star formation parameter that can be more reasonably estimated (with a 0.18 dex scatter in Sim1) from photometric SED fitting.
For all the other parameters, it is remarkable that the bias distributions have a median which is close to 0, and in all cases within far less than the 1σ of the distribution.Despite the degeneracies, the stellar mass bias only has a 0.1 dex scatter (here considered as half of the central 68 percentile of the bias distribution), close to a typical stellar mass uncertainty (0.11 dex on average in this sample) from this type of measurement (e.g.Palmese et al. 2016Palmese et al. , 2020;;Conroy 2013), while AV and δ have a scatter of 0.25 and 0.28, respectively.The SFR scatter is below 0.2 dex, though the SFH parameters have higher scatters: tmax having a scatter of 2.8 and σSFR having a scatter of 4.6.
Considering the 2D distributions, some parameters show correlations that require further investigation.However, certain correlations, such as the correlation between the tmax and SFR residuals, are expected due to the fact that SFR is calculated directly from tmax and σSFR.Also the correlation between the AV and SFR residuals can be understood since the dust attenuation in UV corrects the observed UV luminosity to get SFR.Given the direction of the degeneracy, the correlation is due to a confusion between "red and dead", low SFR galaxies with low dust content, and young, high SFR, high dust content galaxies.The AV -SFR correlation will be explored in the following subsection.
We note that the SNR used in Sim1 is worse than that from the more realistic distributions in Sim2, therefore also the recovered parameters will have more scatter.By measuring the parameters residuals with Sim2 we recover scatters that are half or less than those we present for Sim1.A similar problem persists for the SFH FWHM.

Dust parameters and star formation rate
In Figure 3, we show the recovered distribution of AV and δ residuals for our model sample Sim1.The contour lines represent the percentiles of the joint distribution from 10 to 90%, and reveal a correlation between ∆AV and ∆δ, meaning that an overestimation in AV is likely to be paired with an overestimation in δ, and the same with an underestimation.A degeneracy is expected due to the definition of these parameters.Because the attenuation is always better estimated when comparing shorter wavelengths to the rest frame NIR, where the effect of dust is minimal, is it reasonable to expect that an overestimation of AV will also lead to an overestimation of δ, while the NIR data points provide a "zero-point' to the estimated attenuation law.One can also understand this in the following way: the total dust luminosity (the difference between the dust-free SED and the observed SED) is better constrained than either AV or δ individually.In order to keep the total dust luminosity constant, the slope of the attenuation law must decrease, which amounts to making δ larger.Another relevant aspect to note to understand why we cannot recover the dust parameters very precisely is that at the redshifts considered here, with u being the shortest band, we are not sampling the rest-frame far UV.The lack of FUV will degrade dust estimation (especially since IR covering the typical dust emission wavelengths is not available either).
The fact that the contour lines elongate along the x-axis close to ∆AV = 0 are due to regions where the input value of AV ≃ 0, as the colors in Fig. 3 show.A low AV corresponds to cases with little dust (or where the effect of dust on the SED is less pronounced), so that δ is difficult to constrain, and it appears unconstrained from our fits.Because low AV values (truth AV = 0.1) lie at one of the edges of the prior (which only allows positive values down to 0), it is expected that most low AV galaxies will lie on the ∆AV < 0 side rather than above 0.
We next seek to explore how the dust residual degeneracies change with the amount of data used.Now moving to our Sim2 sample of galaxies to represent a more realistic galaxy population and considering different combinations of bands used, we generate Figure 4.The effect of more data is immediately evident, the inclusion of more bands (to clarify, we are comparing less versus more bands both using Sim2, these changes are not due to changing from Sim1 to Sim2) both reduces the size of the scatter in ∆AV and ∆δ (from 0.25 and 0.26 in the ugriz case to 0.12 and 0.17 respectively in the All Bands case) while also lifting the degeneracy between the two parameters, as can be determined by comparing the shapes of the contour lines.We find that the inclusion of W1 and W2 provides a slightly better result for the scatter (0.12 for AV and 0.15 for δ) than the addition of W3 and W3 (0.14 for AV and 0.17 for δ) to ugrizJHKs for both AV and δ, while the bias is similar and consistent with 0 for the two cases.In addition, the correlation between the true AV value and the AV residual is still present in these plots, though it is not as extreme in the Figure 3 case, likely due to fewer low AV galaxies being present in the Sim2 case (as the distribution of values peaks above 0.1), amongst other effects.The correlation between the true AV value and the AV residual also lessens as more bands are added.This is likely due to the fact that increasing bands allows us to measure both δ and AV more precisely, reducing the error compared to the range of values they can take and hence ensuring against an artificial trend.Other degeneracies in residuals exist, notably a degeneracy in the residual of AV and of SFR.In addition, both of these residuals seem to correlate with stellar mass as well.
To investigate this relationship, we compare the degeneracies with four different photometric band sets: the standard ugriz, ugriz with the near-infrared JHKs, ugrizJHKs with the WISE bands, and ugriz with WISE bands.Plotting these three parameters together in the four different cases yields Figure 5.In each panel, we show the SFR residual versus the AV residual for all the galaxies in Sim1, while the stellar mass residual is given by the color.
The most immediately evident conclusion from Figure 5 is that a degeneracy between AV and SFR residuals exists, but it can be significantly reduced by introducing additional photometric bands.It is clear that the driving factor in reducing the degeneracy is the WISE bands.Comparing ugriz versus ugrizJHKs and ugriz+WISE versus All Bands subplots we can see there is hardly any difference between these 2 sets of plots, but the difference driven by the introduction of the WISE bands is drastic.The ugriz case results in a ∆AV scatter of 0.35 and using ugrizJHKs gives a scatter of 0.29, while including WISE bands reduces scatter by 80%, ugriz+WISE gives a scatter of 0.062, the inclusion of all bands reduces it slightly further to 0.059.This means an increase in SFR compared to the truth value can produce a galaxy SED which is similar to the true one if coupled with a larger effect of dust attenuation, namely a larger AV , and vice versa with decrements in SFR and AV .The addition of the two 3-5 µm WISE channels significantly helps pinning down the attenuation law further in the infrared, where dust absorption is minimal, while the inclusion of the longer wavelengths at > 10µm W3 and W4 helps to constrain the dust content via its emission.
For what concerns the stellar mass, a degeneracy with SFR was already clear in Figure 2, and obvious since the stellar mass is an integral of SFR over time, while a correlation with AV is less obvious.A slight degeneracy is present as photometry from a more massive, older, redder, and dust-free galaxy can be confused with a less massive, younger, dusty galaxy.Once the dust properties are more precisely constrained (via the introduction of additional bands), this degeneracy also appears to be broken, because the dust correction needed to convert the UV luminosity to SFR becomes more precise.Simultaneously, the stellar mass scatter is significantly reduced.In the ugriz bands case, the distribution has a median and 68th percentiles of 0.03 +0.20 −0.19 , which is reduced to 0.001 +0.08 −0.08 when all bands are included.Scatter and bias are similar whether JHKs or WISE bands are added to ugriz, so it does not seem that either of these 2 sets of bands is really driving the improvement in stellar mass.
We next consider the effect that addition of bands has on the general recovery of both AV and δ.In Figure 6, we plot The degeneracy between A V and SFR residuals is evident in the upper two plots, but the degeneracy is reduced by the use of additional bands.The stellar mass residual is also plotted via the color of the points, demonstrating that any A V -stellar mass degeneracy is also lifted by the use of additional bands.SFR and stellar mass are correlated in all plots, which is not surprising as stellar mass is dependent on the SFR by definition.
the difference between truth and recovered values for AV and δ against their truth values.We bin the AV and δ range into five bins, and plot the median ∆parameter value in each bin as lines.The shaded regions represent the 16th to 84th percentile in each bin.From this figure, we can immediately discern several positive effects that the inclusion of extra bands has on recovered values.In both the AV and δ plots the ugriz relationship exhibits some degree of correlation between the truth value and the residuals.We can see that in both cases the addition of more bands reduces this correlation, nearly eliminating it entirely in the AV case.We find that the Pearson coefficient goes from 0.25 in the ugriz case to 0.08 when all bands are used.For δ, the Pearson coefficient is close to 0.18 in all cases.For both cases we note that the change of the median or 68th percentile of the distribution in the residuals over the entire AV or δ range consiered is typically of the order of 0.1 or less, and significantly smaller than both the distribution scatter and the dynamical range of the model values.With optical bands alone the dust parameters posteriors are significantly prior dominated, and the effect of the prior edges show up as correlations in this parameter space.These results show that analyses attempting to recover dust attenuation parameters from these optical bands alone (e.g.Duarte et al. 2022) should be taken with caution.The inclusion of more bands tightens the distribution in residuals around the zero value.Each successive addition of more bands reduces the width of the 16th-84th percentile regions.Thus, with the addition of more bands not only is the median residual value closer to zero, the scatter, as expected, also shrinks.The per-sistence of the δ -∆δ correlation even when using all of the bands available, although to a lesser extent, (see Figure 6) reflects the difficulty of fitting δ, resulting broad posteriors for each galaxy.Thus, the trend seen in the δ plot is more likely a demonstration of the difficulty of fitting δ in general, rather than an indication of an inherent bias introduced during fitting.Nevertheless, scatter and bias are significantly improved by the addition of bands.
Finally, we can consider the recovery of all parameters of interest and how this is affected by choice of bands.In Figure 7 we plot the overall distribution of residuals for each choice of bands.Again, we notice that in general the distance from zero and the overall spread of the distribution decreases monotonically with the inclusion of additional bands.Considering that ugriz + JHKs and ugriz+WISE have roughly the same number of bands, it is interesting that the latter outperforms the former in 3 out of the 5 parameters considered, while results are comparable for the other two parameters.This implies that the infrared region that WISE covers is particularly important for recovering AV and SFR, while it is less relevant for δ and stellar mass.This connection is unsurprising considering the correlation noted in Figure 5.

Recovering relations between AV and δ
Recognizing and constraining relations between dust parameters can help us understand the evolution of galaxies and the physics behind dust formation.However, since dust parameters are typically extracted using SED fitting codes rather Figure 6.A V and δ residuals versus their truth value for the four different band combinations considered in this work.The central dotted lines give the median value, while the shaded region is the 1σ confidence interval.Inclusion of additional bands narrows the distribution of residuals around zero, indicating more accurate estimates.In addition, the prior tends to cause the distribution of residuals to be non-symmetric about zero and shift with increasing truth value; the inclusion of more bands significantly alleviates this effect.Note that for bins on the edge of the plot, we extend the value taken at the center of the bin to the edges of the plot to represent the entire range of dust values considered; this means that the data does not spontaneously flatten out as the plot seems to imply.than being a direct observable from large galaxy surveys, it is a concern that any correlation measured is due to some intrinsic bias in the methodology rather than a real physical phenomenon.Qin et al. ( 2022) claim that AV − δ correlations are driven by such biases, and test their hypothesis by inputting a flat AV − δ and demonstrating that they recover a non-flat correlation.Here, we attempt to recreate their results.
As we are now specifically checking a result that could affect measurements from real data, rather than exploring how biases in one parameter arise in conjunction with other parameters, we decide to make our galaxy sample a closer representation of realistic galaxies, and use the Sim2 dataset.
For our first experiment, we enter a flat AV − δ relation (what we call Sim2a) with delta fixed at −0.125, similar to the value of δ used by Qin et al. ( 2022): −0.2.The models are calculated in the same way as the Sim2 dataset but we fix delta to be −0.125 rather than selecting multiple random values.We then consider the recovered values of AV and δ after running BAGPIPES using all bands.Using a Markov Chain Monte-Carlo (MCMC) algorithm, we fit a line to the recovered data.Since the number of data points proved too large for the algorithm to compute in a reasonable amount of time, we took a random subset of two hundred galaxies.Our results are shown in the top panel of Figure 8.
Here we recover a slope of m = 0.016 +0.007 −0.006 .This is consistent with zero within 2.5σ.We note that the scatter in δ is larger at low AV , likely due to the fact that lower AV means that δ is harder to constrain.Since the prior on δ is set from -1.6 to 0.4, and the true δ is set to -0.125, a prior-dominated delta measurement is much more likely to have a median below its truth value due to the edge of the prior.To test this, we repeat the analysis with a lower δ value, -0.4.The results of this attempt are shown in the bottom panel of Figure 8.
In this case, the measured slope is m = 0.010 +0.006 −0.006 .The slope is now 37.5% lower, and consistent with zero within 1.7σ.Hence, while both slopes are consistent with a flat relation, it seems that prior effects are driving only minor difference in the recovered slope.In Qin et al. ( 2022) the recovered slope for the experiment similar to this recovers a larger slope than our results.It is unclear what is driving this difference in results.Possible explanations are different priors, the different SED-fitting codes (CIGALE versus BAGPIPES) or the different SFH models (their exponential declining versus our lognormal) used.However, the drastically improved results in our method speaks well to its accuracy in recovering correct parameter values.
Beyond simply not introducing spurious correlations, we can demonstrate that we can recover other dust parameter relations using our BAGPIPES runs.Salim & Narayanan (2020) derives the following relation from galaxies in their sample: which is the combination of equations ( 14) and ( 15) in that work.To determine if we can recover this distribution, we enter this relationship as truth values for another sample of galaxies (what we call Sim2b) and ran BAGPIPES using all of our considered bands.Our recovered fit and the expected fit is shown in Figure 9. Defining m as the coefficient in front of log(AV ) and b as the constant at the end of the equation, we get that m = 1.26 +0.02 −0.02 , which varies from the Salim & Narayanan (2020) by only 1.2σ.Simultaneously we measure that b = -0.054+0.003 −0.003 .This is a statistically significant discrepancy, being slightly higher than 3σ, however it is one that is again likely driven by prior effects.As the majority of the points in this distribution are close to the upper δ prior cutoff or the lower AV cut off, this has the effect of pushing the distribution down and to the left as described above.This explains why the MCMC algorithm returns a curve shifted down from the expected curve.However, this is only introduces an error of 0.01 which is negligible when considered against the typical galaxis uncertainties in AV , which are on average an order of magnitude larger.
One can also consider the effect that choice of bands has on recovery of AV -δ distributions.Repeating the above flat  AV -δ experiment for the four combinations of bands we have been using thus far, we generate Figure 10.Note that the axis has been broken, as the recovered slope for ugriz is about 30 times larger than any other band combination.With the exception of ugriz, none of these deviations from zero are statistically significant as they are all below 3σ.

The 2175 Å bump
As mentioned above, the additional absorption in the attenuation curve known as the 2175 Å bump is not considered for the majority of our analysis, as we have found previously that it had little effect on recovered parameters.We seek to reproduce and quantify that claim with model galaxies in this section.Immediately evident is the difficulty in recovering the value of the bump.Even when using all of the bands considered in this work, we are unable to reliably recover the value of B, instead returning a distribution peaked at the center of the prior (as expected since the posterior will be the same as prior and we are plotting the median); see Figure 11.This behaviour is expected, since to better constrain the bump we would need more bands or spectroscopy at the location of the bump, rather than at longer wavelengths as we are testing here.
Though the bump is poorly constrained itself, it does not introduce any additional error into the other parameters considered.Considering Figure 12, one can determine that the residual distributions for parameters of interest are statistically consistent when fit with and without a bump.Hence, while the bump seems to be poorly fit, including it in a SED fitting process should not harm the other parameters considered.Although the effect is small, we still suggest including the bump modeling (even if UV is not covered) because allowing for that extra attenuation does change the total dust luminosity that needs to be matched to IR luminosity for energy balance.

CONCLUSIONS
In this work, we have presented an analysis of the spectral energy distributions of a range of model galaxies using the BAGPIPES software.Our analysis of galaxy parameter residuals after simulating model galaxies and fitting their SEDs revealed a degeneracy between AV and SFR (as is to be expected, AV is a dust correction to SFR), which can be decreased with the inclusion of additional photometric bands.Specifically, one can reduce ∆AV from -0.10 +0.38 −0.40 to 0.01 +0.11   −0.08 and ∆SFR from -0.12 +0.26 −0.43 to -0.01 +0.06 −0.06 by adding dust IR constraints to a case with neither rest-frame UV nor dust IR.This trend generally extends to all parameters of interest.Comparing recovered residuals for AV , δ, stellar mass, While the medians of all of these distributions are statistically consistent, the inclusion of emission parameters serves to increase scatter in the residuals.For example, the scatter in SFR is 0.038 for both the default and bump cases, while it increases to 0.055 when emission is considered.
SFR, and sSFR revealed the reduction in both mean residual and scatter for these parameters.The inclusion of more bands reduces the correlation between entered AV or δ and their residual.This means that the addition of more bands reduces prior effects that can introduce spurious relations in one's results.We next demonstrated that BAGPIPES does not introduce systematic biases when fitting for dust parameters.In our tests using a flat δ − AV relation and using all bands, we always measured any deviation from expected values as being less than 3σ.Even in the scenario where only optical bands are used, any minimal correlations that we find are never as steep as the observed ones between AV and δ (Salim & Narayanan 2020), making it hard to reproduce the observed correlations with spurious fitting problems alone.Similarly, BAGPIPES is able to recover physical dust parameters with distributions following those found in previous works.Entering a measured AV −δ relationship from Salim & Narayanan (2020) as our truth values, we measured a recovered curve statistically consistent with the expected distribution, though the intercept is mildly biased.This discrepancy is interpreted as due to prior effects and is much smaller than typical statistical uncertainties on δ.
Finally, we demonstrated that including an additional degree of freedom, in the form of a dust bump, should not drastically affect one's determination of relationships between various galaxy parameters, as its inclusion does not introduce any additional degeneracies into the fits.
Our findings are indicative of the fact that correlations such as those found in Meldorf et al. (2022) between the SN Hubble residuals and dust parameters are unlikely to arise from BAGPIPES fitting problems or spurious relations between parameters due to degeneracies, as we do not find evidence for large biases or spurious relations.Similar results have been found in works such as Boquien et al. (2022).
It is important to note that some simplifying assumptions have been made in this work.In particular, we have produced

Figure 1 .
Figure 1.Corner plot of each of the 2D and 1D posteriors for all 7 parameters considered for an example galaxy, using ugrizJHKs bands.The truth value for each parameter is plotted as the blue line, while the 16th, 50th, and 84th percentile of the recovered posterior are shown in the diagonal 1D plots as dashed vertical lines.The truth values for this galaxy were A V = 0.525, δ = −0.975,σ SFR = 14.0 Gyr, log 10 (M⋆/M ⊙ ) = 10.27,Z/Z ⊙ = 1.0, tmax = 10.0 Gyr, z = 0.601.

Figure 2 .
Figure2.Corner plot of the residuals (the truth value minus the recovered value) for each of the seven parameters considered, using ugrizJHKs bands.Each solid line represents zero (i.e. an accurate recovery of the input parameter), while the dotted lines in the 1dimensional histograms represent the median and the 32nd and 68th percentile of the distribution.In order to retain clarity, a ∆log(SFR) < 1 dex cut has been applied, removing 66 of 19,200 galaxies.

Figure 3 .
Figure 3. Recovered A V and δ residuals (the truth value minus the recovered value) using ugrizJHKs bands.The truth value for A V is given by the color bar, while the black dashed lines indicate zero residual, i.e. a perfect recovery of the truth value.The black solid lines indicate the percentiles of the joint distribution from 10 to 90% in steps of 10%.

Figure 4 .
Figure 4. Recovered A V and δ residuals (the truth value minus the recovered value) using Sim2 galaxies in six combinations of bands.The truth value for A V is given by the color bar and the black solid lines indicate the percentiles of the joint distribution from 10 to 90% in steps on 10%.

Figure 5 .
Figure 5. Residuals of A V versus residuals of SFR, where each subplot represents a different choice of photometric bands.The 20, 40, 60, and 80 percent contour lines are plotted in black.The degeneracy between A V and SFR residuals is evident in the upper two plots, but the degeneracy is reduced by the use of additional bands.The stellar mass residual is also plotted via the color of the points, demonstrating that any A V -stellar mass degeneracy is also lifted by the use of additional bands.SFR and stellar mass are correlated in all plots, which is not surprising as stellar mass is dependent on the SFR by definition.

Figure 7 .
Figure 7.The median (point) and 16th and 84th percentiles (error bars) of the distribution of five galaxy parameter residuals (the truth value minus the recovered value), for four different choices of bands.Note how the inclusion of more bands both drives the median of the distribution towards zero and narrows the entire distribution.

Figure 8 .
Figure 8. Results of inputting a flat A V − δ relationship into BAGPIPES using all of the bands in this work and fitting a line to the output.In the top panel, we consider an input of constant δ = −0.125 and in the bottom panel δ = −0.4,represented by red lines.In both cases, the recovered slope (black line) is statistically consistent with zero.

Figure 9 .
Figure9.A similar analysis to that of Figure8, though the input relationship (red line) is that derived inSalim & Narayanan (2020), again using all bands this work considers.The fitted curve (black line) is statistically consistent with their results in terms of the slope of the logarithmic curve, but the intercept is shifted by a mildly significant amount.

Figure 10 .
Figure 10.The recovered slope when inputting a flat A V -δ relationship at δ = −0.4 for four different band combinations.With the exception of ugriz, all of these results are statistically consistent with zero and thus no spurious correlation is being introduced.

Figure 11 .
Figure 11.Distribution of entered versus recovered bump values.Though a uniform distribution of bump values is entered, the recovered bump values tend to occur in the center of the prior, indicating that the bump is extremely poorly constrained.

Figure 12 .
Figure12.Residuals for the parameters of interest for different attenuation curve parameterizations, both without and with bump, as well as including IR dust emission modeling without the bump.

Figure A1 .
FigureA1.A corner plot of the residuals of parameters of interest, fit using all bands with Sim2 model galaxies, with the inclusion of dust emission parameters.This plot demonstrates the difficulty in recovering U min .

Table 1 .
Host galaxy parameters values entered into BAGPIPES to simulate galaxies in the Sim1 set (third column), along with the priors assumed in the SED fitting (last column).