Strategies for Obtaining Robust Spectral Energy Distribution Fitting Parameters for Galaxies at z ∼ 1 and z ∼ 2 in the Absence of Infrared Data

Robust estimation of star formation rates (SFRs) at higher redshifts (z ≳ 1) using UV–optical–near-infrared (NIR) photometry is contingent on the ability of spectral energy distribution (SED) fitting to constrain the dust attenuation, stellar metallicity, and star formation history (SFH) simultaneously. IR-derived dust luminosities can help break the degeneracy between these parameters, but IR data are often not available. Here, we explore strategies for SED fitting at z ≳ 1 in the absence of IR data using a sample of log M * > 10.2 star-forming galaxies from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) for which 24 μm data are available. We adopt the total IR luminosity (L TIR) obtained from 24 μm as the “ground truth,” which allows us to assess how well it can be recovered (as L dust) from UV–optical–NIR SED fitting. We test a variety of dust attenuation models, stellar population synthesis models, metallicity assumptions, and SFHs separately to identify which assumptions maximize the agreement (correlation and linearity) between L TIR and L dust. We find that a flexible dust attenuation law performs best. For stellar populations, we find that Bruzual & Charlot models are favored over those of Eldridge et al. Fixing the stellar metallicity at solar value is preferred to other fixed values or leaving it as a free parameter. For SFHs, we find that minimizing the variability in the recent (<100 Myr) SFH improves the agreement with L TIR. Finally, we provide a catalog of galaxy parameters (including M * and SFR) for CANDELS galaxies with logM*>8 and 0.7 < z < 1.3, obtained using the models we found to be the most robust.


Introduction
The current star formation rate (SFR) and stellar mass are two of the most fundamental properties of galaxies.For an accurate understanding of galaxy formation and evolution, robust constraints of both properties for statistically significant samples are needed.Galaxy properties are often estimated by fitting libraries of model spectra to observed spectral energy distributions (SEDs), a process referred to as SED fitting (for reviews, see Walcher et al. 2011;Conroy 2013).In recent years, the practice of SED fitting has been greatly refined.State-of-the-art SED fitting codes such as CIGALE (Boquien et al. 2019), PROSPECTOR (Leja et al. 2017;Johnson et al. 2021), MIRKWOOD (Gilda et al. 2021), BAGPIPES (Carnall et al. 2018), BEAGLE (Chevallard & Charlot 2016), and SEDfit (Sawicki 2012) are highly flexible, efficient, and accessible, and offer a diverse set of models from which to choose.
In modeling galaxy SEDs, there are three crucial components: a stellar population synthesis (SPS) model, a star formation history (SFH), and a dust attenuation model.The adopted SPS model is convolved with an SFH, which specifies the SFR as a function of time, to produce a model spectrum of a composite stellar population with stars of various ages.The spectrum is then attenuated according to some dust law(s).
Model SEDs with stellar populations of various ages and metallicities and attenuated by varying levels of dust are then fit to observed galaxy SEDs or spectra.
Age, metallicity, and dust each act to redden the spectrum (Bell & de Jong 2001;Carter et al. 2009).As the signal-tonoise ratio (S/N) of the photometry degrades, the age-dustmetallicity degeneracy is exacerbated and the accuracy and/or precision of the SED fitting may be compromised.Progress has been made in identifying sets of models which are physically motivated and help to minimize the effects of degeneracy in the fitting (e.g., Pacifici et al. 2012).It is well established that galaxies host a wide variety of dust attenuation laws with different UV bump strengths and steepness, and allowing for this flexibility in the SED fitting produces more accurate results, at least at low redshifts (Salim et al. 2018).SPS models are now available that include the effects of binary stars, which may be needed to infer accurate physical properties especially for young stellar populations (Eldridge et al. 2017;Stanway & Eldridge 2018).Well-motivated SFH models have proven effective in correcting for bias in stellar mass estimates due to outshining (Buat et al. 2014;Michałowski et al. 2014;Simha et al. 2014;Sorba & Sawicki 2015;Salim et al. 2016;Lower et al. 2020), and can produce self-consistent measurements of the stellar mass growth and SFR density over time (Leja et al. 2019b(Leja et al. , 2020)).Recently, state-of-the-art SFHs have been utilized to reconcile the long-standing discrepancy between the normalization of the star-forming main sequence from simulations with that inferred from observations (Leja et al. 2022).The SED fitting community continues to make great strides in improving the accuracy and reliability of galaxy spectral modeling.
Improvements in accuracy are not necessarily accompanied by improvements in precision, however, and model degeneracies will tend to inflate the uncertainties of physical parameters even when systematic offsets are well accounted for.Complementing UV and optical photometry with mid-and far-infrared photometry, which is sensitive to the emission of dust heated by young stars, can be used to break the age-dustmetallicity degeneracy and greatly improve the constraining power of the SED fitting (e.g., da Cunha et al. 2008;Buat et al. 2014).Unfortunately, however, IR surveys with Spitzer and Herschel were limited in their completeness and sensitivity even for massive galaxies, and for most galaxies one is limited to fitting only their UV-optical-near-infrared (NIR) photometry where age, dust, and metallicity must be simultaneously constrained.It is therefore crucial to understand which model choices ensure that the SED fitting is as precise as possible.The need for good general SED fitting practices is especially important with the imminent influx of survey data from JWST as well as the Vera C. Rubin Observatory.
In this work, we identify some good practices for obtaining robust physical properties from SED fitting of UV-optical-NIR photometry without constraints from the mid-or far-infrared for galaxies at z ∼ 1 and z ∼ 2. Our sample consists of starforming galaxies with well-sampled UV-optical-NIR SEDs and solid mid-infrared photometry.We explore the effects of different models for dust, stellar metallicity, SPS model, and SFHs on the reliability of physical parameters inferred from the SED fitting.In particular, we determine a set of models for which the total energy absorbed by dust in the UV-optical-NIR (L dust ) best agrees (in terms of correlation and linearity) with the total IR luminosity (L TIR ) inferred from mid-infrared photometry.
This paper is organized as follows.We describe our data and sample selection in Section 2. We detail our analysis methods, including choice of SED fitting code and choice of models, in Section 3. We then present our results in Section 4, discuss these results in Section 6, and conclude in Section 7. We assume a flat WMAP7 cosmology (H 0 = 70 km s −1 Mpc −1 and Ω m = 0.27) throughout the paper.
To account for photometric zero-point errors, for each band except Spitzer/MIPS 24 μm (which we do not include in the SED fitting) we add to the reported error an amount equal to 10% of the reported flux.To be included in the SED fitting, we require a galaxy to have a detection in at least one band covering rest-frame UV and be detected in at least five total bands in the rest-frame UV or rest-frame optical-NIR ranges (again not counting Spitzer/MIPS 24 μm).A total of 140,773 galaxies across all four fields satisfy these photometric requirements.

Redshift Selection
CANDELS provides photometric as well as spectroscopic redshifts.We adopt the spectroscopic redshifts where available and the photometric redshifts otherwise; ∼10% of galaxies have spectroscopic redshifts.We impose a redshift cut of z > 0.7 to ensure that the photometry covers the rest-frame UV.We furthermore exclude galaxies with z > 2.3 as the sample size diminishes significantly above this point.In order to discern a potential redshift dependence in our results, we will in some instances divide our sample into two redshift windows, one at 0.7 < z < 1.3 (z ∼ 1) and the other at 1.7 < z < 2.3 (z ∼ 2); we otherwise include all galaxies in the range 0.7 < z < 2.3.We also require galaxies to have a CANDELS redshift that differs from the 3D-HST (Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016) redshift by no more than 0.4 (i.e., |Δz| < 0.4) to avoid including galaxies with poor photometric redshifts.Imposing the upper and lower redshift limits (0.7 < z < 2.3) reduces our sample to 76,593 galaxies, and the |Δz| < 0.4 requirement reduces it further to 60,266 galaxies.

Final Sample Selection
We further require that each galaxy have a 24 μm detection with an S/N of at least five, which reduces our sample size to 3652 galaxies.We also impose a specific star formation rate (sSFR) cut (log sSFR > −10) so that we include only nonquenched/quenching star-forming galaxies in our assessments.We note that there are very few 24 μm detections at log sSFR < −10.We additionally impose a stellar mass cut (log M * > 10.2) to ensure completeness at all the redshifts spanned by our sample (Figure 1).We find that our mass completeness limit is consistent with other estimates from the literature (e.g., Chartab et al. 2020).Note that both the sSFRs and the stellar masses used in these cuts are taken from the fiducial SED fits described in Section 3.2.We finally exclude galaxies with very poor fiducial SED fits (reduced χ 2 > 10).Applying the sSFR, stellar mass, and reduced χ 2 cuts reduce our sample to 2622 galaxies.
Figure 2 shows the sSFR versus stellar mass diagrams for the two redshift windows.Shown in gray is our original CANDELS sample with photometric, redshift, Δz, and r 2 c cuts applied, while the 24 μm-detected sample with the additional sSFR and stellar mass cuts (i.e., our final sample) applied is shown as red points.The sample sizes are 1221 at z ∼ 1 and 757 at z ∼ 2.
Figure 3 shows the completeness of the 24 μm photometry as a function of redshift for the high-mass star-forming galaxies in our sample.We define the completeness as the fraction of galaxies with a 24 μm S/N > 5. Notably, the 24 μm completeness decreases with redshift from ∼60% at z ∼ 0.8 to ∼40% at z ∼ 1.5, before rising back up to ∼60% at z ∼ 2. This is due to a shifting of the rest-frame wavelength to ∼8 μm at z ∼ 2 where strong polycyclic aromatic hydrocarbon (PAH) features are present in the IR spectra.Overall, however, the completeness is relatively uniform across the entire redshift range.

Methods
In this section we describe our choice of SED fitting code, our fiducial set of models which are used as a baseline, the various test models, our method for estimating the total IR luminosity (L TIR ) from the IR photometry, and our treatment of redshift-dependent systematic offsets between L TIR and L dust .

SED Fitting
We use the Python-based SED modeling and fitting code CIGALE (version 2022.0;Boquien et al. 2019).CIGALE allows specification of the SFH, SPS model, nebular emission, and dust attenuation curve.CIGALE uses a Bayesian methodology wherein a probability distribution function (PDF) is constructed for each parameter (e.g., stellar mass) based on the model likelihoods.The nominal value of the parameter is taken to be the mean of this PDF, and the standard deviation of the PDF gives the error (see, e.g., Smith & Hayward 2015).
CIGALE comes with a variety of built-in SED fitting "modules" which encompass the various models used in the SED fitting.For example, two of the dust attenuation modules used in this work are "dustatt_calzleit" and "dustatt_powerlaw," which parameterize the dust attenuation curve in different ways.We also use two modules that were created specifically for this work by modifying existing ones.The first of these is the "dustatt_calzleit_2comp" module, which was created by modifying the built-in "dustatt_calzleit" to allow the slope of the young and old attenuation curves to vary independently.The second is the "bpass" module, which was created by modifying the built-in "bc03" module to use the Binary Population and Spectral Synthesis (BPASS; version 2.2) SPS model library (Eldridge et al. 2017) instead of the Bruzual & Charlot (2003; BC03 hereafter) SPS model library.

SED Fitting Models
In this section we describe the fiducial models used in the SED fitting, which we define as the models that produce the tightest agreement between L TIR and L dust .We also describe each variation on our fiducial model.We independently vary the treatment for dust, metallicity, choice of SPS model, and SFH to isolate their effects on the recovered L dust from one another.
All models include nebular emission (continuum and lines) for which we fix the gas-phase metallicity at Z * = 0.014 and the ionization parameter at log U = −3.4.While inclusion of some nebular emission is necessary to achieve unbiased results (Pacifici et al. 2015;Salim et al. 2016), we find that varying the nebular parameters (e.g., using log U = −2.8)has a negligible impact on L dust and thus the SFR.
For all models explored in this work, we include all available photometry with the exception of IRAC Channel 3 for all galaxies at z < 1.1 and IRAC Channel 4 for all galaxies at z < 1.5 to exclude passbands that include hot dust emission.

Dust Attenuation
A summary of the different dust models tested is shown in Table 1.In our fiducial model the dust attenuation curve has a range of slopes and UV bump strengths, and is labeled as "flexible."This flexible dust law follows the Noll et al. (2009) parameterization as formulated by Salim et al. (2018), and is based on the functional form of the Calzetti et al. (2000) law, but has three important changes that provide it with greater flexibility: (1) it is a two-population law (young and old population; see Charlot & Fall 2000;Wuyts et al. 2009Wuyts et al. , 2011) ) with higher attenuation around young stars (10 Myr and (2) the law can have a range of steepnesses (δ), i.e., the relative attenuation in the UV compared to the optical, which we allow to vary from significantly steeper than the Calzetti curve (δ = −1.5) to shallower than the Calzetti curve (δ = 0.4; note that the Calzetti curve has δ = 0); and (3) a 2175 Å UV bump can be added, for which we allow the strength (i.e., the amplitude) to vary from zero (no bump) to 4.5 (the bump in the Milky Way extinction curve has a strength of three).We allow the normalization E (B − V ) of the curve for the young population to vary from zero to one, with finer intervals at lower values.
We test three dust models in addition to the fiducial model.We first explore the simple assumption of a fixed Calzetti law (Calzetti et al. 2000), which we refer to simply as the "Calzetti" model.This model differs from the fiducial one by fixing the slope δ at zero and bump strength at zero (no bump).The Calzetti model assumes equal attenuation for young and old stars Figure 1.Stellar mass vs. redshift distribution for star-forming galaxies (log sSFR > −10) in our sample with robust IR detections (a 24 μm S/N > 5).To ensure consistent completeness at all redshifts, we include only galaxies with log M * > 10.2 in our final sample.
The next model we consider is a modification of the fiducial dust law where we allow the slope δ of the young and old components to vary independently, which we refer to as the "extraflexible" model.For this model the δ for both the young and old curves have the same range of values as the fiducial model.Like the fiducial model, the reddening of the old curve is normalized to 0.44 times the young curve.
We also consider a flexible attenuation law, which uses as its basis a power-law functional form, rather than the Calzetti one, which we label as "flexible (power law)."We allow the slope to vary from shallow (n = 0.2) to very steep (n = 1.6), again assuming the same slope for the old and young components.We also include a UV bump with the same range of amplitudes as the other models.In CIGALE the power-law curves are normalized in terms of A V instead of E(B − V ) and we allow this normalization to vary from 0 to 2.5 mag, covering a similar range of values in terms of E(B − V ) as the fiducial model.Like the fiducial model, the reddening of the old curve is normalized to 0.44 times the young curve.

Stellar Population Synthesis Models and Metallicity
A summary of our stellar population model variations are shown in Table 2. Our SPS models are from BC03, with the stellar metallicity fixed at solar (Z * = 0.02).We label the Figure 2. sSFR vs. stellar mass for the galaxies in the two redshift windows considered in this study.We separate quiescent (plus quenching) and star-forming galaxies using a log sSFR = −10 cut.The entire sample of CANDELS galaxies on which we performed the SED fitting are shown as a gray density plot.Our final sample, for which we require a robust IR detection (a 24 μm S/N > 5), active star formation (log sSFR > −10), and a relatively high mass (log M * > 10.2), is shown as red points.Current IR observations are not able to detect most of the lower-mass galaxies in HST images.

Normalization CIGALE Module
Flexible (Fiducial) 0 to 4.5 −1.2 to 0.4 Same as young 0.44 We test another three models in addition to the fiducial model.As an alternative to BC03, we consider the BPASS (version 2.2) stellar population models (Eldridge et al. 2017;Stanway & Eldridge 2018).BPASS is notable for including the effects of binary evolution, which are expected to be significant for young, massive stars, and has been used extensively to study the ionizing continua and stellar populations of galaxies at the dawn of reionization.We compare the fiducial model to a BPASS model where the metallicity is also fixed at Z * = 0.02, which we label as "solar Z * BPASS."Using BPASS, we can assess the precision of L dust under different treatments of stellar evolution.
We also test a variation of the fiducial model for which the metallicity is allowed to vary between Z * = 0.004, Z * = 0.008, Z * = 0.02, and Z * = 0.05, and label this model as "free Z * ."We finally test a variation of the fiducial model where the metallicity is fixed at subsolar (Z * = 0.008), which we refer to as "subsolar Z * ."

Star Formation Histories
A summary of the SFH model variations used in this work are shown in Table 3.For the fiducial SFH we employ a parametric double exponential model, which assumes an old population of fixed age with a second (burst) component superimposed.The age of the old population is fixed and set to be ∼500 Myr less than the age of the Universe at the time of observation.The e-folding time of the old population is allowed to vary from 200 Myr to 20 Gyr with a finer grid among shorter e-folding times in order to sample the log sSFR space evenly.We fix the e-folding time of the burst at 20 Gyr (functionally constant); thus our "burst" (in the fiducial model) is not transient but is simply another episode of star formation.We vary the fraction of mass formed in the burst from 0 to 0.5 and vary the age of the burst from 100 Myr to 1 Gyr, again with a finer grid at lower values.Fixing the age of the old population and limiting the burst fraction and burst age prevents outshining of old stars and underestimated stellar masses (Shapley et al. 2005;Pforr et al. 2012;Buat et al. 2014;Michałowski et al. 2014;Salim et al. 2016), and also ensures a more realistic range of colors compared to single population models (Pacifici et al. 2015).The inclusion of the second episode of star formation is needed to allow for high sSFRs not otherwise attainable by a single exponential with a fixed old age (see Ciesla et al. 2017), and effectively produces rising SFHs, which are known to be more realistic than strictly declining models for many high-redshift galaxies (e.g., Reddy et al. 2012).
We test three different SFH models in addition to the fiducial model.For the model labeled as "flexible burst" we allow τ burst to take on values of 10, 200, and 1000 Myr in addition to 20 Gyr.This allows for a wider range of the current SFR at a given mass-weighted age; in other words, for a given burst age, mass fraction, and old e-folding time (which together specify the typical stellar age of the model) different burst e-folding times will lead to different levels of current star formation but with similar stellar population ages.
We include another flexible SFH model for which τ burst is again fixed at 20 Gyr but where we allow for burst ages of 10 and 30 Myr, which allow for younger stellar populations compared to the fiducial model; we label this as the "younger flat bursts" model.Using the flexible burst and younger flat burst models, we will evaluate how the precision of L dust is related to the stochasticity of the recent SFH.
Details of the old (>1 Gyr) SFH are poorly constrained by broadband photometry alone.Nonetheless, we test to see if changing the prescription for the old SFH has a noticeable impact on the agreement between L dust and L TIR .We use the same grid of parameters as for the fiducial run, but with a delayed exponential SFH instead of a declining exponential; we label this as the "delayed old SFH" model.The key difference is that for the delayed SFH, the star formation of the old component rises smoothly, peaks, then declines, whereas for the exponential, the SFH of the old component starts at a maximum value and declines from there.

Total IR Luminosities
To estimate the total IR luminosities from the IR photometry, we make use of the IR templates and associated software of Boquien & Salim (2021), referred to as BOSA.BOSA offers templates based on a sample of 2584 low-redshift star-forming galaxies (galaxies with an active galactic nucleus (AGN) have been excluded) and which have extensive coverage in the IR from the Wide-field Infrared Survey Explorer and Herschel (12-500 μm).
The BOSA templates are notable for including a dependence on the specific SFR as well as IR luminosity.The sSFR dependence helps to ensure the accuracy of L TIR for highredshift galaxies, for which L TIR may be biased when using strictly luminosity-dependent templates (see Elbaz et al. 2011).We note that BOSA is capable of choosing templates based on IR color when multiple IR fluxes are provided, but in the case of a single flux point the templates must also be luminosity dependent.The BOSA code requires at least one IR flux point and some estimate of the sSFR as input.We supply the sSFR from the fiducial SED fitting to obtain L TIR .We also reproduced all results using a fixed log sSFR of −9 (representative of high-redshift star-forming galaxies) and find no change in our conclusions.
We note that 24 μm observed at z = 1 corresponds to a restframe wavelength of 12 μm while at z = 2 the rest-frame wavelength is 8 μm.It has been established that AGN can contribute significantly to galaxy SEDs in the mid-infrared (e.g., Kirkpatrick et al. 2015;Leja et al. 2018), though unobscured (type 1) AGN can have an effect on the UVoptical-NIR photometry as well (e.g., Wuyts et al. 2009).However, in this work we focus mostly on the scatter in L dust versus L TIR , which will be minimized for the best set of models irrespective of any potential moderate AGN contribution.

An Assessment of Redshift-dependent Systematics
Figure 4 shows the difference in L dust estimated from the fiducial SED fitting run and L TIR versus redshift.A polynomial fit (degree = 6), which is shown in Figure 4 as a red line, helps visualize the typical offset at different redshifts.At z  1.4, L dust is typically greater than L TIR by ∼0.2 dex.There is a "break" at z ∼ 1.4 where the systematic offset reverses, beyond which L TIR is greater than L dust by ∼0.25 dex.The "break" at z ∼ 1.4 corresponds roughly to the redshift where strong PAH features are redshifted into the Spitzer/MIPS 24 μm bandpass.The redshift-dependent offsets between L dust and L TIR may be partially due to inflexibility on the part of the sSFR-dependent templates to account for the variation in PAH features as well.Unidentified systematics in the SED fitting may also contribute to the offset (e.g., Leja et al. 2019b).
To try and discern whether the redshift-dependent offsets of Figure 4 are due primarily to the SED fitting or to the IR templates, we perform a mock SED fitting test, which aims to determine how well can SED fitting (still using only UV-NIR constraints) recover "known" IR luminosities.To do so, we first repeat the fiducial SED fitting (see Section 3.2), but include L TIR as a strong constraint.Choosing very small errors (1% of L TIR ) guarantees that the resulting L dust is equal to L TIR .The best-fit fluxes from this IR-constrained fitting now have to conform with this L TIR .For the mock fitting, the mock fluxes are drawn from a Gaussian centered on the best-fit flux assuming the same S/N as the original observations (so, mock error = (observed error/observed flux) × best-fit flux).The mock fitting is performed on the mock fluxes and mock errors, using fiducial models and with no IR constraints.Any differences in the mock L dust , which we refer to as L Mock , and L TIR will then be entirely attributable to the SED fitting.
Figure 5 shows L Mock − L TIR versus redshift.As with Figure 4, we show a polynomial fit (degree = 6) as a red line.Comparing to Figure 4, we see that the differences obtained from the mocks are smaller than the real differences at all redshifts, and lack most of the inflections present in the real data.Thus we conclude that the redshift-dependent systematic offsets between the actual L dust and L TIR are partially due to limitations in the IR templates and partially arising from the SED fitting.The latter are on the order of 0.1 dex, which should not preclude us from addressing the goals of this study.
The goal of this work is to identify best practices for UVoptical-NIR SED fitting using the scatter in L dust versus L TIR , rather than the bulk offsets.We therefore correct for the redshift- Note.For all models, the age of the old population is fixed at ∼500 Myr less than the age of the Universe at the time of observation, and the burst fractions are limited to the range [0, 0.5].
Figure 4. IR luminosity residuals (L dust − L TIR ) vs. redshift for the final sample of galaxies.L dust is the predicted dust luminosity from the fiducial SED fitting, whereas L dust is the same quantity determined from 24 μm observations using IR templates.The red line represents a polynomial fit (degree = 6), which we use to correct redshift-dependent systematics so we can focus on the strength of the correlation between L dust and L dust .These offsets are smaller than the real offsets (Figure 4) at all redshifts, suggesting that the real offsets may be primarily driven by the systematics in deriving L TIR from 24 μm observations using IR templates.
dependent systematics using a polynomial fit (degree = 6), which is shown in Figure 4 as a red line, to place L dust and L TIR on a more equal footing.For each galaxy, we subtract from L dust the value of the polynomial at that galaxy's redshift.The same correction (based on the fiducial L dust ) is applied to all estimates of L dust .

Results
In this section, we test variations on a fiducial SED model in order to evaluate the level of agreement between the dust luminosity (L dust ) predicted from the UV-optical-NIR fitting and the L TIR inferred from the Spitzer 24 μm photometry.We aim to determine which treatments for dust, metallicity, choice of SPS models, and SFH maximize the degree of correlation between L TIR and L dust .We also explore systematics in the estimated stellar mass (log M * ) across different models.
While the goal of the study is to arrive at more precise SFRs, comparing L dust and L TIR instead of SFRs is advantageous for two reasons.First, L TIR includes emission from dust heated from stellar populations of all ages, so we do not need to isolate only the young populations when deriving an SFR estimate from IR emission.Second, on the SED fitting side, the SFR needs to be averaged over some timescale, and it is not obvious which timescale is most appropriate when comparing to SFRs obtained from L TIR (or L TIR in combination with L UV in the popular hybrid method).On the other hand, if L dust determined from the stellar continuum SED fitting agrees with L TIR , then the SED fitting SFRs should be robust too.
To evaluate agreement we primarily look at the standard deviation (σ) of the difference between L dust and L TIR and not the offsets (which are set to zero for the fiducial run by definition, see Section 3.4).For stellar mass comparisons we discuss primarily the mean relative offset (Δ) between the various SED fitting schemes and the fiducial one, since we have no stellar masses that can serve as a ground truth.

Dust Attenuation
We first explore how variations in the assumed dust attenuation law affect the correlation between L dust and L TIR .Figure 6 shows the L TIR and L dust comparison for the fiducial model (top panel) alongside the comparisons for the Calzetti (second panel), extraflexible (third panel), and flexible powerlaw (bottom panel) models.
The Calzetti model has a single slope for the attenuation curve, allows no UV bump, and assumes the same attenuation for young and old stars, unlike the fiducial model, which varies the slope and bump strength and applies extra attenuation to young stars.The Calzetti model results in strongly underestimated L dust among galaxies with relatively low L TIR , which is largely attributable to the absence of steeper slopes.Omitting this low-L dust tail, the rest of the galaxies agree with L TIR relatively well, though the scatter is still higher at z ∼ 2 than for the fiducial dust model.
The extraflexible model allows the attenuation curves applied to old and young stars to vary independently from one another, unlike for the fiducial model where the young and old slopes are kept the same.For extraflexible attenuation, there is a slight decrease in scatter at z ∼ 2. Overall, the extraflexible model shows only a marginal improvement over the flexible fiducial model while being much more computationally expensive due to the expanded grid, thus we do not recommend its use.Assuming a flexible power-law parameterization instead of the fiducial model significantly increases the scatter with respect to the fiducial run in both redshift windows.The fiducial model outperforms both the simpler and more complex models, as well as the alternative power-law model.In this and similar subsequent panels, the scatter σ is the standard deviation in log L dustlog L TIR .

Stellar Population Synthesis Models
We now consider how varying assumptions about the stellar populations impact the precision of L dust .The top panel of Figure 7 shows the results for the fiducial solar Z * BC03 model, while the bottom panel shows the results for the solar Z * BPASS model.We see that using the BPASS SPS models instead of the fiducial BC03 SPS models results in increased scatter between L dust and L TIR at both redshifts.BPASS models appear to predict systematically lower L dust compared to BC03 at both redshifts, despite the stronger ionizing flux of the BPASS models compared to the BC03 models for a given stellar age.The lower values of L dust may be due to the lower stellar masses predicted by BPASS (see Figure 10), resulting in fewer stars and less overall heating.There is also a much more prominent nonlinearity in L dust versus L TIR for the BPASS models at z ∼ 2.

Stellar Metallicity
We now consider how varying the assumptions about the stellar metallicity impacts the precision of L dust .We show comparisons of the different models in Figure 8.The fiducial model, which fixes the metallicity at solar (Z * = 0.02), is shown in the top panel.Allowing the stellar metallicity to vary (free Z * , middle panel) results in increased scatter compared to the fiducial model at z ∼ 1, which is likely attributable to dustmetallicity degeneracy, but marginally decreased scatter at z ∼ 2. The slightly lower scatter at z ∼ 2 for the free metallicity models may suggest that the uncertainties in L dust at higher redshifts are dominated by the photometric S/N rather than modeling uncertainties.The net scatter is lower for the fiducial model compared to the free Z * model, so we find the fiducial model to be favored overall.
While a fixed solar metallicity is often assumed in SED fitting for simplicity, it is not obviously the best choice for galaxies at high redshift, in particular at z ∼ 2 where the gasphase metallicity is lower than at z ∼ 0 by as much as a factor of two (see, e.g., Maiolino et al. 2008;Wuyts et al. 2016).However, we find that for the subsolar Z * model (bottom panel), which fixes the metallicity at Z * = 0.008 (or 40% solar), the scatter is increased compared to the fiducial model and there is a systematic shift to higher L dust ; the systematic shift seems to occur because, with less metallicity, the dust attenuation must increase to match the observed colors.Interestingly, the opposite occurs when the metallicity is fixed  For fixed metallicity, the choice of metallicity is a major systematic.The fiducial model, which fixes the metallicity at solar, is favored overall.Free metallicity models may be slightly disfavored at z ∼ 1 due to dust-metallicity degeneracy, while at z ∼ 2 the S/N of the observations has a greater impact on the uncertainty of L dust than this degeneracy.
at supersolar (Z * = 0.05), which causes a systematic shift to lower L dust .We also do not find any improvement over the fiducial model at either redshift when the metallicity is fixed at even lower values (e.g., Z * = 0.004).

Star Formation History
Finally, we explore how changing the type of SFH model affects the precision of L dust .Figure 9 shows the fiducial model (top panel), the flexible burst model (second panel), the younger flat bursts model (third panel), and the delayed old SFH model (bottom panel).The fiducial model is a twocomponent parameterization featuring a declining exponential old component with a fixed age superimposed with a second episode of recent star formation with an essentially constant SFR and varying mass fraction.We find that the flexible burst model, which differs from fiducial one in that the burst is allowed to decline (variable τ burst ) instead of remaining flat, increases the scatter compared to the fiducial model at both redshifts.The increase in scatter for the flexible burst model is most significant for relatively low-L TIR galaxies at z ∼ 1.
The younger flat bursts model, which differs from the fiducial model in that it allows for burst ages younger than 100 Myr (10 Myr and 30 Myr), significantly increases the scatter among the highest-L TIR galaxies (which tend to be fit with these younger bursts) compared to the fiducial model.Finally, we consider the model that treats the old component as a delayed exponential.We find only marginally increased scatter at z ∼ 1 and marginally decreased scatter at z ∼ 2 compared to the fiducial model; therefore, there is little practical difference between the two, i.e., both represent good choices.

Stellar Masses
So far we have focused on dust luminosity and therefore the SFRs.We now consider how the estimated stellar mass depends, in the relative sense, on the choice of model.The top row of Figure 10 compares the stellar masses from the Calzetti (left), extraflexible (middle), and flexible power-law (right) models to stellar mass from the fiducial model.In terms of the mean offset, the masses from both the Calzetti and extraflexible dust curves are consistent with the fiducial masses of all individual galaxies.The flexible power-law model predicts higher stellar masses by 0.21 dex on average compared to the fiducial model, which is significant considering that the typical error in the stellar mass (fiducial) for galaxies in our sample is ∼0.05 dex.A similar result was found by Lo Faro et al. (2017), who report that masses based on a power-law dust model are 0.15 dex higher than the ones from the Calzetti law.
The middle row of Figure 10 compares stellar masses from the free Z * BC03 (left), subsolar Z * BC03 (middle), and solar Z * BPASS (right) models to stellar mass from the solar Z * BC03 fiducial model.Changing the metallicity treatment does not lead to significant systematic differences in the masses.However, BPASS predicts systematically lower masses compared to BC03, by 0.17 dex on average; this is significant compared to the typical error in the stellar mass (fiducial) for galaxies in our sample, which is ∼0.05 dex.
The bottom row of Figure 10 compares stellar masses from the flexible burst (left), younger flat bursts (middle), and delayed old SFH (right) models to stellar masses from the fiducial model.The masses are quite similar between each model, with minimal scatter and offset.While the flexible burst model has a noticeable one-sided scatter, the bulk of the sample has essentially identical masses.Notably, the delayed old SFH model also produces masses that are highly consistent with the fiducial model.

Catalog of SED Fitting Galaxy Parameters
To facilitate the use of our results in future studies we present a catalog of SED fitting galaxy parameters obtained using our fiducial model.Specifically, we publish our estimates of the stellar mass (log solar masses), SFR (log solar masses per year), far-ultraviolet (FUV) and V-band attenuation (magnitudes), and dust luminosity (log solar luminosities) together with their uncertainties.The catalog can be found online in both the journal (see Appendix and Table 4) and at the website https://salims.pages.iu.edu/candels.Included in the catalog are all galaxies for which we performed SED fitting with redshifts in the range 0.7 < z < 2.3, with reduced 10 r 2 c < , and with differences in redshift between CANDELS and 3D-HST of less than 0.4 (|Δz| < 0.4), amounting to 63,266 galaxies total.The CANDELS identifier (integer index used in all CANDELS catalogs) and CANDELS subfields (e.g., COSMOS and GOODS-S) are also included alongside the CANDELS redshifts used in the fitting (i.e., z best ).

Discussion
In this section we discuss the most significant results outlined in the previous section and place them into the context Figure 10.Comparison of the log stellar mass (log M * ) for various models to log M * from our fiducial model.On average, the Calzetti and extraflexible models predict similar masses as the fiducial model while the flexible power-law model predicts systematically higher masses compared to the fiducial model (top row).The choice of metallicity does not result in systematic differences in mass; BPASS predicts lower masses compared to BC03 (middle row).Masses are highly consistent between each SFH model (bottom row).The scatter σ is the standard deviation in log M *log M *,fiducial while the offset Δ is the mean of log M *log M *,fiducial . of previous studies.We also suggest some additional good practices for consideration in future studies.

Dust Attenuation
Our results are generally consistent with previous findings that galaxy attenuation curves are diverse and not always well described by a fixed law (Kriek & Conroy 2013;Salmon et al. 2016;Salim et al. 2018;Theios et al. 2019;Salim & Narayanan 2020).While a flexible dust treatment is essential for accurate SED fitting results, it is worth considering whether some constraint on the attenuation law could improve the agreement between L dust and L TIR .For instance, it has been shown that there is a relationship between V-band attenuation and the steepness of the attenuation curve (Salim & Narayanan 2020).It has also been shown that a correlation exists between the power-law slope and the strength of the UV bump (Kriek & Conroy 2013).These correlations can be added as a constraint in the SED fitting, still allowing for flexibility but with an informed constraint on the parameter space; this saves computation time and could potentially reduce degeneracy in the SED fitting.We tested both of these constraints separately but found that they offer no improvement in the agreement between L dust and L TIR for our sample.However, such constraints could still potentially be useful if rest-frame UV photometry (essential for constraining the shape of the attenuation curve) is not available.The δ-A V and δ-bump relations could also serve as physically motivated model distributions in Markov Chain Monte Carlo (MCMC) fitting frameworks such as Prospector, which uses a gridless MCMClike approach to fit galaxy SEDs and includes the Kriek & Conroy (2013) slope-bump relation as a built-in distribution (see Leja et al. 2017).
We note that the shortest observed wavelengths included in the photometry are those spanned by the U or u bands (hereafter referred to collectively as U bands).The peak wavelength of the U bands corresponds roughly to 3500 Å, which, at our lowest redshift (z = 0.7), corresponds to a restframe wavelength of ∼2000 Å which is in the near-ultraviolet (NUV) range.One would ideally include both FUV and NUV coverage in the SED fitting since the UV color is sensitive to dust.It is only above roughly z = 1 that the U bands cover the rest-frame FUV; despite this, however, we find that the agreement between L dust and L TIR is actually better at z ∼ 1 than at z ∼ 2 regardless, possibly because of the higher S/N of the z ∼ 1 photometry.

Stellar Population Synthesis Models and Metallicity
We find that fixing the stellar metallicity is a preferred choice given the quality of the data we have in CANDELS and this type of sample; it may not be the best option in some other circumstances, e.g., for galaxies with old stellar populations where the dust has a lesser impact on the galaxy SED and the age and metallicity can be better constrained simultaneously (see Dorman et al. 2003).Our choice of solar metallicity may not be appropriate for galaxies outside the range of our sample, i.e., at lower masses (log M e < 10.2) or higher redshifts (z > 2.3), for which the galaxies have yet to undergo significant enrichment.It has also been found that the assumption of fixed metallicity may lead to systematic biases in the stellar mass, albeit for galaxies below our mass range of log M * < 10.2 (Mitchell et al. 2013).
A more sophisticated approach than assuming a single fixed metallicity would be to impose a constraint on the metallicity according to some mass-dependent and perhaps redshiftdependent relation.Such a technique was employed by Leja et al. (2019b), who, using the Prospector code, sampled the mass-metallicity space according to a normal distribution with mean and σ from the Gallazzi et al. (2005) relation for local galaxies, clipping the distribution to the range 0.00021 < Z * < 0.031 (−1.98 < log (Z * /Z e ) < 0.19); this method effectively weights the models toward the local relation yet allows the fits to take on significantly lower values to account for possible evolution with redshift.
Another approach would be to fix the stellar metallicity in the fitting for a given mass and redshift based on an assumed mass-metallicity-redshift relation, effectively removing the metallicity as a free parameter (since redshift is fixed and the mass is tied to the SFH), which has been proposed by recent studies of the mass-metallicity relation in galaxies (see Thorne et al. 2022).Using VANDELS, Cullen et al. (2019) found that the mass-metallicity relation for z > 2.5 star-forming galaxies is consistent in slope with the local relation, but is shifted in zero-point to lower values by ∼0.6 dex.A redshift-dependent mass-metallicity relation can therefore be introduced into the SED fitting by simply adopting the local relation and scaling it according to the galaxy's redshift.However, direct constraints on the stellar metallicity at intermediate redshifts (i.e., 0.3 < z < 2.5) are severely limited due to the need for highresolution rest-frame UV spectroscopy, which renders it difficult to determine the appropriate scaling with redshift.One possible solution is to use the redshift evolution in the gasphase metallicity-mass relation as a proxy for the scaling of the stellar mass-metallicity relation (see, e.g., Maiolino et al. 2008;Wuyts et al. 2016;Thorne et al. 2022).We note, however, that the resolution of the SPS model metallicity grid for some libraries may be too restrictive to implement a detailed redshift dependence, and our results may suggest that no redshift dependence is even needed, at least for z  2.3, since we find, somewhat surprisingly, that solar metallicity models produce more reliable L dust compared to subsolar models at both z ∼ 1 and z ∼ 2.
Whereas BC03 SPS models are known to compare favorably to other SPS models at low redshifts (e.g., Conroy & Gunn 2010;Hansson et al. 2012;Zibetti et al. 2013), one might expect binary effects to be more relevant for galaxies at high redshift which have younger stellar populations and thus greater fractions of high-mass main-sequence stars.It is thus interesting that BPASS should be disfavored by our results.The effects of binary stars on the stellar emission may be obfuscated if the uncertainties in the UV photometry (in which regime the effects are most pronounced) are sufficiently large, while the treatment of other stellar evolution phases in the SPS models, such as thermally pulsating AGB stars or extreme horizontal branch stars, may have a stronger impact on the SED fitting overall (Maraston 2005;van der Wel et al. 2006;Conroy et al. 2009;Walcher et al. 2011;Conroy 2013;Zibetti et al. 2013;McGaugh & Schombert 2014;Villaume et al. 2015).It may be that BC03 models provide a more realistic treatment of certain aspects of nonbinary stellar evolution compared to BPASS, which could explain why our results favor the BC03 models.However, we note that while BPASS is disfavored in terms of L dust , we have no independent estimate of the stellar mass and so cannot rule out the possibility that BPASS produces more accurate masses.We also cannot rule out the possibility that other properties determined from the fits, such as the strength of the ionizing continuum or specific spectral indices, may also be more realistic when using BPASS.

Star Formation History
The results of our tests on SFH suggest that minimizing the variability of the very recent (100 Myr) SFH produces the most reliable L dust .We find that this is attributable to age-dust degeneracy: when allowing for sharp drops in the recent SFH, models with less dust and a recent drop in the SFR will produce similar colors to a model with constant or rising recent SFR and more dust.The fiducial flat burst model resolves this issue by forcing the recent (<1 Gyr) SFH to either increase or remain flat, effectively avoiding the age-dust degeneracy.Our findings echo those of Wuyts et al. (2011), who find, using single τ SFH models with variable age, that allowing for very small e-folding times (300 Myr) results in less accurate SFRs.However, Curtis-Lake et al. (2021) find that for relatively low-mass galaxies (8  log M *  10.5) at very high redshifts (z ∼ 5), broadband photometry becomes sensitive to fluctuations in the SFR on shorter timescales (∼10 Myr), so in such a regime it may be necessary to allow for more stochasticity in the recent SFH to obtain accurate physical parameters.
We also note that the so-called "poststarburst" or "E+A" galaxies feature recent (<1 Gyr) quenching in their SFHs (e.g., French 2021).A recent quenching is not allowed by our fiducial SFH.It may be necessary to allow for more flexibility in the recent SFH to obtain accurate physical parameters for poststarbursts specifically (see Suess et al. 2022).However, poststarbursts are relatively rare (see French 2021) and they should not significantly affect the results from SED fitting of statistical samples.Furthermore, poststarbursts are sometimes selected based on the ratio of Hα-to-UV SFR, so it is unclear whether L dust would even be significantly affected as the SED fitting may not be sensitive to changes in the SFR on such short timescales, i.e., on the order of 10 Myr (Flores Velázquez et al. 2021;French 2021).
Our findings regarding the use of delayed versus exponential parameterizations for the old SFH are consistent with past studies, which find that the details of the old SFH are effectively unconstrained in broadband SED fitting, though parametric forms where the age of the old population is allowed to vary freely are susceptible to underestimating stellar masses due to outshining (Salim et al. 2016;Carnall et al. 2019;Leja et al. 2019a;Lower et al. 2020).So long as the outshining bias is accounted for (which we accomplish by fixing the age of the old population to near the maximum possible time), the systematic uncertainties in M * and L dust are dominated by the dust, metallicity, stellar population models, and/or recent SFH, rather than the ancient SFH.

Published Catalogs of Physical Properties from SED Fitting
We note that value-added catalogs of physical properties determined via UV-NIR SED fitting are already available for thousands of galaxies at a wide range of redshifts in certain CANDELS fields.One example is the Santini et al. (2015) catalog, which includes physical properties for thousands of GOODS-S and UDS galaxies compiled from various groups who used different SED fitting codes and model assumptions.
Notably, some groups allow the stellar metallicity to vary while others keep it fixed; some groups use SPS models other than that of BC03; some groups use single-component SFHs; and a fixed Calzetti et al. (2000) dust law is assumed in all but one instance (in which case the group also allowed for the SMC dust law).Another public catalog is the Barro et al. (2019) catalog for GOODS-N, which includes stellar masses and SFRs based on both UV-NIR SED fitting as well as combined UV + IR SFRs.Their SED fitting method is based on that of Wuyts et al. (2011).They find the different estimates of the SFR to be generally consistent for galaxies with both UV and IR photometry.While they hold the stellar metallicity fixed and adopt BC03 SPS models, they also assume a fixed Calzetti et al. (2000) dust law and use single-component parametric SFHs with variable ages.For users of these catalogs we caution that the assumption of single-component SFHs may lead to underestimated stellar masses (e.g., Lower et al. 2020), and the assumption of a fixed dust law or variable stellar metallicity may result in biased SFRs for some galaxies as we have demonstrated (see also Salmon et al. 2016;Salim et al. 2018).

Applicability to Other Data Sets
It is worth considering whether the results obtained in this study, which features a sample of high-mass galaxies with abundant UV-optical photometry from CANDELS, is applicable in the hypothetical case of shallower observational data.To test this, we perform a set of mock SED fitting where for each galaxy we artificially increase the photometric errors by a factor of 6.3 (thus mimicking 2 mag brighter limiting magnitudes), then resample the fluxes using the larger errors.We then fit a selection of our models to determine whether their relative performance is affected by the decrease in S/N. Figure 11 shows the results of this mock fitting for the fiducial, free Z * , Calzetti dust, and flexible burst models.We find the relative performance of each model is consistent with our main results: the fiducial model performs best overall, though free Z * outperforms it for z ∼ 2 galaxies.The consistency in results despite the larger errors suggests that our results are applicable even to fields with significantly shallower photometry than that typical of CANDELS.

Choice of SED Fitting Code
We do not explore the impact of different SED fitting codes in this work.However, different codes come prepackaged with different model libraries and the underlying fitting "machinery" (i.e., the inference method used to estimate properties from fitted models) may also differ.These differences may result in systematic uncertainties and, especially in the case of SFRs and dust attenuation (e.g., A V ), systematic offsets which depend on the choice of SED fitting code (Pacifici et al. 2023).We refer to Pacifici et al. (2023) for a discussion of methodological best practices which can be used to mitigate the impact of codedependent systematics.

Conclusions
In this work, we are able to identify some good practices for SED fitting at z  1 to improve constraints on the SFR and reduce the effects of age-dust-metallicity degeneracy.Using a sample of star-forming galaxies with abundant UV-optical-NIR photometry from CANDELS and coverage in the IR with Spitzer, we evaluate which SED fitting models produce the best agreement (in terms of scatter) between L dust from the UVoptical-NIR fitting and L TIR inferred from the 24 μm photometry.Our use of dust (total IR) luminosity as a point of comparison, rather than SFR directly, avoids systematics due to the heating of dust by old stars and the choice of an SFR timescale.
We find generally that for SED fitting there exists a "sweet spot" between assuming a fixed model on the one hand, and having too much freedom in the models on the other hand.Going too far in either direction can result in the introduction of noise or systematics.Whereas exploring different flavors of SED fitting allowed us to identify better choices, we have yet to arrive at a precise match between L dust and L TIR in terms of zero-point differences and the skew (i.e., linearity).Our main conclusions are as follows.
1. Allowing for dust attenuation curve slopes to be flexible, and including the slopes that would be considered rather steep (SMC-like), is essential for getting unbiased L dust (and thus presumably SFRs).On the other hand, assuming a fixed shallow curve can lead to a number of galaxies with severely underestimated dust luminosities.2. The flexible dust attenuation model of Noll et al. (2009), which takes the Calzetti curve and allows its steepness to vary, produces better results than when the form of the underlying curve is assumed to be a pure power law.3. Assuming a fixed stellar metallicity, in particular solar metallicity, produces better results than allowing the metallicity to vary. 4. BC03 SPS models produce L dust estimates that are in tighter agreement with observations than the ones produced by BPASS models.5.An SFH model in which the recent burst of star formation has constant SFR over at least the past 100 Myr produces better results than allowing for a more recent burst, or allowing the burst to decline.6.Whether one assumes an exponentially declining SFH for the first (older) component or the so-called delayed SFH makes essentially no difference, presumably because the shape of the old SFH is poorly constrained by the broadband photometry.7. Stellar masses between different models generally agree to within few hundredths of a dex.The exception is when the attenuation curve slopes are based on power laws, which results in 0.21 dex higher masses than our fiducial dust model; when using BPASS models, we find 0.17 dex lower masses than with BC03 models.8.The relative performances of the different models remain similar when shallower photometry is available.
We also make publicly available our estimates of the stellar mass and SFR, among other parameters.

Figure 5 .
Figure5.IR luminosity residuals (L dust − L TIR ) vs. redshift for the final sample of galaxies.L Mock represents L dust from the SED fitting in which real fluxes have been replaced with the mock fluxes that must by construction produce L TIR .The red line represents a polynomial fit (degree = 6), which helps visualize the redshift-dependent systematics arising from SED fitting alone.These offsets are smaller than the real offsets (Figure4) at all redshifts, suggesting that the real offsets may be primarily driven by the systematics in deriving L TIR from 24 μm observations using IR templates.

Figure 6 .
Figure6.Comparison of L dust from UV-optical-NIR SED fitting to L TIR for different dust attenuation models.The fiducial model outperforms both the simpler and more complex models, as well as the alternative power-law model.In this and similar subsequent panels, the scatter σ is the standard deviation in log L dustlog L TIR .

Figure 7 .
Figure 7.Comparison of L dust from UV-optical-NIR SED fitting L TIR for different stellar population models.The BC03 models are favored compared to the BPASS models.The BPASS models predict systematically lower L dust compared to BC03 at both redshifts.

Figure 8 .
Figure8.Comparison of L dust from UV-optical-NIR SED fitting L TIR for different stellar metallicity models.For fixed metallicity, the choice of metallicity is a major systematic.The fiducial model, which fixes the metallicity at solar, is favored overall.Free metallicity models may be slightly disfavored at z ∼ 1 due to dust-metallicity degeneracy, while at z ∼ 2 the S/N of the observations has a greater impact on the uncertainty of L dust than this degeneracy.

Figure 9 .
Figure 9.Comparison of L dust from UV-optical-NIR SED fitting to L TIR for different types of SFH models.The fiducial model outperforms the other models that allow for more variability in the recent SFH.Adopting a delayed exponential parameterization for the old component, instead of an exponential form, has little effect.

Table 1
Dust Attenuation Law Models Note.Module "dustatt calzleit 2comp" is not part of the standard CIGALE distribution.
fiducial model as "solar Z * ."A Chabrier (2003) initial mass function (IMF) is assumed for all models.
For all models we assume a Chabrier (2003) IMF with an upper mass limit of 100 M e .Module "bpass" is not part of the standard CIGALE distribution.