Panchromatic Photometry of Low-redshift, Massive Galaxies Selected from SDSS Stripe 82

We select galaxies in the 79 deg² SDSS Stripe 82 region that is covered by HerS (Viero et al. 2014) with the Herschel/SPIRE instrument. The SDSS DR16 database provides the spectroscopic redshift, spectral classification (Bolton et al. 2012), and morphological classification (Lupton et al. 2001) of each target. Among the three spectral types (GALAXY, QSO, and STAR), we only select sources classified as GALAXY. We also require the morphological type of the targets to be GALAXY ⁵ when selecting targets in CasJobs. ⁶

Figure 1 shows for the Stripe 82 galaxies the distribution of redshift versus stellar mass (M_*) derived from GSWLC-2 (Salim et al. 2018). A lower redshift cut of z > 0.01 is required in GSWLC-2 to exclude very nearby galaxies whose distances are poorly indicated by their redshift (Tully et al. 2016). Two notable large-scale structures arise from overdensities at z ≈ 0.07 and 0.13 (Figure 1(a)). As one of our primary objectives is to incorporate FIR measurements into the SED analysis in order to obtain more robust SFRs, we require that a significant fraction of the sources be detected at least in the Herschel/SPIRE 250 μm band. Figure 1(b) color codes the data points by the detection rate (${f}_{\det }$) of SPIRE 250 μm as given in the Herschel Extra-galactic Legacy Project (HELP) catalog (Shirley et al. 2021), which considers flux densities with a signal-to-noise ratio (S/N) > 2 as detections. Our sample is defined by the 2681 galaxies with M_* > 10¹⁰ M_⊙ and z < 0.11, which have an average 250 μm detection rate of ${f}_{\det }\gtrsim 0.55$. Beyond z = 0.12, the average detection rate in redshift bins drops to 0.5 or lower. We exclude 13 galaxies with very low surface brightnesses that lack SDSS isophotal measurements down to a surface brightness level of μ_r = 25 $\mathrm{mag}\,{\mathrm{arcsec}}^{-2}$, which are needed for our aperture definition (Section 3), and an additional two galaxies that lie on the edge of the GALEX field, whose distortion effects strongly impact the UV data.

Zoom In Zoom Out Reset image size

To facilitate scientific analysis using our database, we classify the sample into isolated (interaction flag = 0) and interacting (interaction flag = 1) galaxies. We consider a galaxy as interacting or merging if it contains a physically associated companion, which we define as another galaxy with a spectroscopic redshift difference Δz_spec < 0.002 (<600 km s⁻¹) within a projected separation of d ≤ 50 kpc (e.g., Ellison et al. 2010; Calderón-Castillo et al. 2019). Lacking a reliable spectroscopic redshift, we regard it as interacting or merging if it exhibits obvious tidal features connecting to the main galaxy, and its photometric redshift (z_phot), considering its uncertainty, is consistent with the z_spec of the main galaxy. We additionally require that the two targets differ by less than 3 mag in the r band. This ensures that we do not miss galaxy pairs with mass ratios as low as 10:1, a conventional threshold for minor mergers, while at the same time excluding the multitude of confusing, fainter nearby sources that may be associated with clumpy substructures from the main galaxy. The sample contains 96 galaxies with an interacting companion that has M_* ≤ 10¹⁰ M_⊙, and hence formally does not satisfy the stellar mass cut of the main sample. We retain these pairs, so long as the companion is not more than 3 mag fainter than the primary galaxy. Since many of these fainter companions only have photometric redshifts, which have relatively large uncertainties, in these instances, if they do not have spectroscopic redshifts, we simply assume that their redshifts are identical to the spectroscopic redshift of the primary companion.

More than half (∼53%) of the galaxies have strong emission lines in their optical spectra. We diagnose their photoionization source based on the [O iii] λ5007/Hβ versus [N ii] λ6584/Hα line-intensity ratio diagram (Baldwin et al. 1981), adopting the classification criteria of Kewley et al. (2001) and Kauffmann et al. (2003). As we present in Figure 2, among the 1472 galaxies for which all four emission lines are detected with S/N > 3 according to the MPA–JHU catalog, 850 are classified as star-forming galaxies, 247 as AGNs, and 375 as composites (star-forming component mixed with an AGN component). The rest (∼47%) of the sample has little to no detectable line emission, consisting mainly of truly quiescent galaxies with minimal star formation and low-ionization nuclear emission-line regions (LINERs; Heckman 1980), which are predominantly highly sub-Eddington AGNs (Ho 2008, 2009). Since such weak emission-line sources will contribute negligibly to the broadband SEDs, for the purposes of this study will simply consider them inactive galaxies. The above statistics do not account for 17 sources within the redshift cut that have spectral type QSO in the SDSS database, a generic label for sources that exhibit broad emission lines. Although these objects were excluded from GSWLC-2 and therefore do not have determined stellar masses, their stellar masses are expected to exceed 10¹⁰ M_⊙ and therefore satisfy the selection criteria of our sample. From the catalog of Liu et al. (2019), these 17 broad-line AGNs have black hole masses M_BH ≈ 10^6.1–10^8.2 M_⊙, which, according to the M_BH−M_* scaling relation (for all galaxy types) of Greene et al. (2020), correspond to M_* ≈ 10^9.7−10¹¹ M_⊙. In summary, the final sample has 2781 galaxies (Table 1), consisting of 2668 objects that satisfy the original stellar mass cut (M_* > 10¹⁰ M_⊙), 96 supplementary objects of lower stellar mass that are companions of members from the original sample, and 17 broad-line (type 1) AGNs included for completeness.

Our analysis involves panchromatic photometry of low-redshift, massive galaxies selected from the HerS region of Stripe 82. To this end, we assemble the following sets of multiwavelength images.

• 1.  
  Optical. We extract the urgiz images from Stripe 82 ⁷ , whose coadded images are about 2 mag deeper than typical SDSS single-scan images, reaching a limit surface brightness of ${\mu }_{r}\approx 28.5\,\mathrm{mag}\,{\mathrm{arcsec}}^{-2}$ (Fliri & Trujillo 2016). The median seeing ranges from FWHM ∼131 in u to FWHM ≈ 104 in z.
• 2.  
  UV. GALEX conducted an all-sky imaging survey in the FUV (central wavelength ∼1539 Å) and near-UV (NUV; central wavelength ∼2316 Å) bands with a PSF resolution of FWHM = 42 and 53, respectively (Martin et al. 2005). Although the General Release 6 and 7 (GR6+7; Bianchi et al. 2014, 2017) provide Kron (1980) elliptical aperture magnitudes, the measurements are affected significantly by distortions on the edge of the detector (Morrissey et al. 2007) and source blending (Salim et al. 2016). Instead of applying a statistical correction to the GR6+7 data (e.g., Salim et al. 2016), we perform our own aperture photometry after source deblending to measure the total galaxy fluxes using a consistent, common aperture across the SED (Section 3). We download the GALEX image tiles of the Deep, Medium, or All-sky Imaging Surveys (DIS, MIS, and AIS, respectively; Morrissey et al. 2007) from CasJob. ⁸ If a target is covered by different image tiles, we choose the one that has the longest exposure time. To avoid the distortion effect, we do not accept the NUV image if the target is >053 from the image center.
• 3.  
  NIR and MIR. The NIR data come from images taken in the J (1.24 μm), H (1.66 μm), and K_s (2.16 μm) bands from 2MASS (Skrutskie et al. 2006). Although the spatial resolution of the images is moderate (FWHM ≈ 3''), most of the galaxies are still sufficiently resolved that their fluxes are underestimated significantly in the 2MASS Point Source Catalog (PSC; Cutri et al. 2003; Appendix A). Meanwhile, the 2MASS Extended Source Catalog (XSC; Jarrett et al. 2000) only  contains the measurements of ∼60% of our targets. The AllWISE Atlas of the WISE (Wright et al. 2010) all-sky survey provides match-filtered, coadded images in the W1 (3.4 μm), W2 (4.6 μm), W3 (12 μm), and W4 (22 μm) bands with FWHM resolutions of 61, 68, 74, and 12'', respectively. The AllWISE catalog provides fluxes from PSF fitting as well as aperture photometry, but we find that these measurements underestimate the total flux even for marginally resolved galaxies (Appendix A). The fluxes from the "unWISE" catalog (Lang et al. 2016b) derived from forced photometry based on the SDSS r-band profile of the galaxies closely follow our measurements (Appendix A). However, only 76% of our targets have W4 measurements in the unWISE catalog.The above considerations compel us to reanalyze all the NIR and MIR data in a uniform, self-consistent manner, following the same precepts applied to the optical and UV data. We download the 2MASS and AllWISE images, including the AllWISE uncertainty maps, from the NASA/IPAC Infrared Science Archive. ⁹
• 4.  
  FIR. HerS (Viero et al. 2014) observed an area of 79 deg² in the SDSS Stripe 82 region, covering (α, δ)_J2000 = (0^h54^m, −2°) to (2^h24^m, +2°). The SPIRE imaging photometer (Griffin et al. 2010) was used to survey in the 250, 350, and 500 μm bands to an average depth of 13.0, 12.9, and 14.8 mJy beam⁻¹, respectively. Although the rest of the Stripe 82 region is covered by HerMES (Oliver et al. 2012), the achieved depth is significantly shallower, and in this study we only focus on the HerS area. Nearly all of our galaxies are unresolved in SPIRE, which has a PSF FWHM of 182, 252, and 363 at 250, 350, and 500 μm, respectively. We take advantage of the comprehensive photometry offered by HELP (Shirley et al. 2021), which uses a Bayesian approach to deblend galaxies using prior information from higher-resolution measurements at shorter wavelengths. HELP provides a catalog of forced photometry based on the prior information, as well as a "blind catalog" that does not use the prior information. We crossmatch our galaxies to the forced-photometry catalog using a search radius of 1''. Most of the targets that are not found in this catalog are too faint to be detected by Herschel, although an additional ∼70 measurements were located in the blind catalog after enlarging the search radius to 74, which encompasses >99.7% (3σ) of the observed coordinate separations of the crossmatched Herschel counterparts of the SDSS galaxies. We visually verified that most of them are robust measurements and not contaminated by companion sources (Section 3.3). In total, 96% of the galaxies have reliable Herschel/SPIRE measurements from the HELP catalog, with ∼70% of the targets detected at 250 μ m. For the HELP-unmeasured sources, we estimate flux density upper limits at 250 μm based on the number of observation scans (Viero et al. 2014). We adopt a 2σ upper limit of 20.7 mJy for regions scanned twice and 18.7 mJy for those scanned three times.

We take the aperture of a galaxy (with semimajor axis a₂₅ and semiminor axis b₂₅) twice as large as its isophote at 25 $\mathrm{mag}\,{\mathrm{arcsec}}^{-2}$ in r band to ensure a large enough aperture size. We set the dimensions of the image cutouts to 400'' × 400'' if a₂₅ < 50'', and to 1000'' × 1000'' if a₂₅ ≥ 50''. Setting the image size to ≳10 times the galaxy size ensures ample area to compute the background.

The images are preprocessed using the Python package photutils (Bradley et al. 2020). We obtain a preliminary estimate of the background level and its standard deviation using the detect_threshold function to sigma clip the background data with a threshold of 2σ. After smoothing the image to suppress the noise with a 2D circular Gaussian kernel with a FWHM equal to the image resolution, we run detect_sources and generate a segmentation map. A source is considered detected if at least five connecting pixels have a flux more than two times the standard deviation above the background. We then use source_properties to convert the segmentation map into a source list, which includes information such as the centroid, semimajor axis, semiminor axis, and position angle (PA) of the galaxy. With this geometric information in hand, we grow the semimajor and semiminor axes three times and derive the detection mask. An example is shown in Figure 3(a). To determine and remove the final background, we apply the mask and fit the background pixels with a third-order 2D polynomial model.

Zoom In Zoom Out Reset image size

The image segmentation of the target may be blended with nearby sources. We use the function deblend_sources to separate the potentially blended contaminants, setting nlevels ¹⁰ = 32 and contrast ¹¹ = 0.0001 to break up the original segments into smaller segments aggressively based on the saddle point between the flux peaks (Sazonova et al. 2021; Zhao et al. 2022). If the companion sources are deblended, we generate the deblend mask for them. The final mask of the image is a combination of the detection mask without masking the target and the deblend mask (e.g., Figure 3(b)).

In order to measure the total flux of the target, we need to interpolate nearby masked pixels. With the mask applied, we first derive the surface brightness profile of the target using the task ellipse in Pyraf ¹² (Tody 1993). We then use the 2D surface brightness model generated for the target to interpolate the masked pixels that contain significant emissions from the target. We run ellipse twice. In the first run, we fix the center of the ellipse to the centroid of the target according to its SDSS coordinates and fit the ellipticity and PA of the isophotes in logarithmic steps of 0.1 in radius. We choose the isophote with an average surface brightness two times the standard deviation above the background to represent the shape of the galaxy. In the second run, we fix the ellipticity and PA to the values determined in the first run and fit the isophotes in linear steps of 1 pixel to sample the surface brightness profile of the galaxy densely. The 2D surface brightness model is generated from this output, which is used to replace the masked pixels, after adding Gaussian noise equivalent to the standard deviation of the background (Figure 3(c)). In order to define the region to interpolate the masked pixels, we need to determine the outer boundary where the target emission drops to the level of the background noise fluctuations, which is estimated following the methodology of Li et al. (2011). We define the boundary of the galaxy by the semimajor (a_profile) and semiminor (b_profile) axis of the isophote whose surface brightness reaches the background noise fluctuation.

Some galaxies suffer from such severe star contamination that the above standard deblending procedure fails to be effective. Under these circumstances, we use a new method developed to remove the contaminating emission from unsaturated stars (Appendix B). As a brief summary, when the standard deblending in the UV or IR image fails for a star that contaminates the target galaxy, we predict the flux of the star within the aperture in that band and subtract it from the galaxy photometry. We predict the total magnitude of the star in the UV or IR using a newly developed random forest (RF) regression method that uses the star's optical PSF magnitudes as input. The fraction of the star's flux within the aperture is calculated based on its position relative to that of the target galaxy and the PSF light distribution in the respective band. Isolated stars serve as the training sample, using SDSS PSF magnitudes as input to take advantage of the optical data, which have the highest resolution among all the bands and accurate fluxes from the PSF models. We estimate the magnitude of the star in the GALEX, 2MASS, and WISE (W1−W3) bands, which can be predicted based on an R² score greater than 0.94. A small number (30) of galaxies corrupted by nearby, heavily saturated stars are unsalvageable. We removed them from further consideration, noting that this introduces no bias to the final sample.

With the decontaminated, deblended images in hand, we perform aperture photometry across all 14 bands (FUV, NUV, u, g, r, i, z, J, H, K_s , W1, W2, W3, and W4) using a common aperture defined by the isophote at the surface brightness level of μ_r = 25 $\mathrm{mag}\,{\mathrm{arcsec}}^{-2}$. This surface brightness level is trivially satisfied by Stripe 82, which reaches a depth of ${\mu }_{r}\approx 28.5\,\mathrm{mag}\,{\mathrm{arcsec}}^{-2}$ (Fliri & Trujillo 2016). We set the semimajor axis of the aperture to a₂₅ and its semiminor axis to b₂₅ if both axes are larger than three times the FWHM of the PSF; otherwise, the minimum aperture radius is set to 3 times the FWHM of the PSF, sufficient to enclose >90% of the total flux of a point source. This guarantees that most of the flux of compact or edge-on galaxies can be included in the aperture, while obviating the need for aperture correction. Figure 4 shows an example of a galaxy resolved in all the bands. The WISE images are dominated by the confusion noise of the background; source emission that cannot be detected during the preprocessing steps will not bias the aperture photometry.

Zoom In Zoom Out Reset image size

The W3 and W4 bands require special consideration. The coarse PSFs of these two WISE bands (FWHM = 84 and 12'', respectively) imply that many of the galaxies in our sample are expected to be unresolved (Cluver et al. 2014; Rosario et al. 2016). Because the W3 and W4 observations are substantially shallower than those in W1 and W2, our approach of aperture-matched photometry will not be able to detect a significant fraction of the targets. Under these circumstances, PSF fitting should perform better than aperture photometry. We carried out a series of mock tests to determine conservative empirical criteria to recognize unresolved galaxies under conditions appropriate for our sample selection. As detailed in Appendix C, a galaxy can be safely regarded as unresolved in W3 and W4 if, in the r band, R_e is less than 2'' and 4'', respectively, and n ≥ 3.5. By these criteria, 1328 sources are unresolved and unblended in W3; the corresponding number in W4 is 1844. We directly adopt the PSF-fitting measurements from the AllWISE catalog for these sources. AllWISE provides the uncertainty of the PSF fitting even if the target is not significantly detected, which enables us to estimate upper limits for nondetections.

To ascertain the degree to which our choice of aperture captures the total flux of the galaxy, we use GALFIT (Peng et al. 2002, 2010) to generate mock multiband images of a subset of ∼300 relatively isolated galaxies drawn from our parent sample, using as input the best-fit single-component Sérsic profile parameters from Bottrell et al. (2019). After convolving the models with the respective PSF, we add the simulated image to a real background image constructed for each specific band to mimic as faithfully as possible the actual observations. Application of our aperture photometry procedure reveals that the fraction of the input flux recovered is >90% in the UV bands, >90% in g, r, and i, >80% in u, >85% in the z and 2MASS bands, >90% in W1 and W2, and >80% in W3 and W4. For bright objects, such as those with W4 <6.5 mag, the input flux recovery levels are close to unity. However, as objects become fainter, the scatter of the flux recovery increases to a standard deviation of 10%–20% around a median value of ∼1. The underestimation of flux depends on luminosity and is relatively minor for more than 84% of the galaxies when compared to the typical uncertainties (Section 4) associated with stellar mass (0.05 dex) and SFR (0.18 dex).

A minority (∼340) of the galaxies in our sample cannot be deblended using our standard procedure (Section 3.1). These largely consist of advanced mergers and close, projected pairs, as well as a handful of galaxies in very dense environments, such as the centers of rich clusters where multiple galaxies may be crowded together. Under these conditions, aperture photometry (Section 3.2) cannot yield meaningful results, and generally, none of the photometry given in conventional survey catalogs can be trusted. We perform 2D image decomposition to deblend the galaxies using GALFITM (Häußler et al. 2013; Vika et al. 2013), which is a multiwavelength extension of GALFIT (Peng et al. 2002, 2010). The code fits several bands simultaneously, leveraging the structural parameters of the well-resolved images in some bands against those that are more poorly resolved in others. Our fits incorporate 10 of the 14 filters, spanning from the u band to the W2 band. We use the Python package reproject ¹³ to register all the bands to the world coordinate system of the W1 image and rebin them to the lowest common pixel scale of 1375, which is dictated by WISE. Each galaxy is modeled with a single-component Sérsic function. Model parameters across the different bands can be mutually constrained by a Chebyshev polynomial. Consistent with the experience of Häußler et al. (2013), we find that a second-order function describes well the wavelength dependence of R_e and n, with initial guesses taken from the catalog of Bottrell et al. (2019). We keep the axis ratio and PA constant with wavelength, using as initial input a₂₅, b₂₅, and PA₂₅. The magnitude in each band is allowed to vary freely. An example fit is given in Figure 5.

Zoom In Zoom Out Reset image size

The FUV and NUV bands are excluded from the simultaneous multiband fit because they consistently give unstable results, possibly a consequence of the extreme variations in galaxy substructure due to their young stellar populations and the effects of internal dust extinction. We opt, instead, to fit the two UV bands separately, holding constant the structure parameters (R_e , n, axis ratio, and PA) to the values derived for the g band from the simultaneous 10-band fit. Our 2D decomposition also cannot be applied to the W3 and W4 bands, or to any of the three Herschel bands, because the galaxy pairs are too severely blended in these long-wavelength filters. Moreover, we cannot guarantee that the distribution of the dust emission is identical to that of the stellar emission. However, as mentioned in Section 2.2, the HELP forced-photometry catalog provides deblending results based on the optical/NIR prior coordinates, using the full Bayesian posterior probability distribution. We use the HELP results for the blended pairs if the coordinates match. For pairs that can be matched only in the HELP blind-photometry catalog, we visually check their SPIRE 250 μm images and SDSS composite-color images. If the SPIRE source is contaminated by any neighboring objects visible in the SDSS image, we regard the SPIRE photometry as an upper limit for the target. We measure the total flux in W3 and W4 for the blended system within a master aperture, as described below. The total MIR flux is used as a flux upper limit for each component object when fitting its SED (Section 4).

To quantify possible systematic differences between the photometry derived from the model-fitting approach compared to that measured from the aperture-matched technique used for the majority of the sample, we apply both methods to a calibration sample of 50 randomly selected, isolated galaxies that roughly spans the range of r-band magnitudes in the parent sample. The two sets of measurements correlate tightly, although small systematic offsets between them exist. The model-fitting photometry is on average brighter than the aperture-matched photometry by ≤0.1 mag in g, r, i, z, W1, and W2. The offset is larger in the u band (∼0.3 mag), a possible manifestation of Runges phenomenon (Häußler et al. 2022). The J, H, and K_s bands have larger discrepancies of 0.1–0.3 mag owing to their lower sensitivity, with a dispersion comparable to the systematical offsets. We apply these average offsets in the optical and IR bands to the GALFITM-derived multiband photometry for the heavily blended galaxies to mitigate systematic biases with respect to the rest of the sample.

Treating each of the blended galaxies as a single Sérsic component obviously vastly oversimplifies the intricate, complex substructure often displayed by mergers and strongly interacting galaxies. Our approach, albeit crude, suffices to provide an effective first-order, relatively accurate partitioning of the individual fluxes of the constituent members in the blended system. We find that the sum of the fluxes of the individual model components matches the total flux of the blended system to better than ∼0.05 mag in all bands studied. The total, integrated flux of the blended system is computed within a common, master aperture based on the surface brightness sensitivity of the image. We choose the largest a_profile and b_profile (Section 3.1) among all the bands as the final master aperture (see also Shangguan et al. 2019).

The uncertainty of the aperture photometry is

$\begin{eqnarray}&&{\sigma }_{\mathrm{total}}^{2}={\sigma }_{\mathrm{source}}^{2}+{\sigma }_{\mathrm{bkg}}^{2}+{\sigma }_{\mathrm{conf}}^{2}+{\sigma }_{\mathrm{star}}^{2}+{\sigma }_{\mathrm{sys}}^{2},\end{eqnarray} \tag{ 1 }$

where σ_source is the Poisson noise of the photons, σ_bkg is the background noise, σ_conf is the confusion noise, σ_star is introduced by blending stars (Equation (B2)), and σ_sys is the systematic uncertainty of the aperture photometry. We explain each term below.

Assuming that the total flux from the aperture photometry is I, ${\sigma }_{\mathrm{source}}^{2}=I/G$, with G specifying the gain of the image. The background noise can be estimated from an annulus around the target aperture, ${\sigma }_{\mathrm{bkg}}^{2}=N\,{\sigma }_{\mathrm{ann}}^{2}$, where σ_ann is the standard deviation of the background pixels in the annulus and N is the number of pixels in the aperture of the galaxy. As shown in Figure 6(a), the annulus has the same ellipticity and PA as the target's aperture. The inner and outer semimajor axes of the annulus are 1.25 and 1.60 times the semimajor axis of the aperture, respectively, so that the area of the annulus is the same as that of the aperture.

Zoom In Zoom Out Reset image size

As the 2MASS and WISE images are resampled and interpolated, we need to consider their pixel-correlated noise (Jarrett et al. 2000). We adopt

$\begin{eqnarray}&&{\sigma }_{\mathrm{source}}^{2}\,(2\mathrm{MASS})=I/G+{F}_{\mathrm{corrA}}\,{N}^{2}{\sigma }_{\mathrm{ann}}^{2},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{\sigma }_{\mathrm{source}}^{2}\,(\mathrm{WISE})={F}_{\mathrm{corrA}}\,\displaystyle \sum _{i}^{N}{\sigma }_{i}^{2},\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{\sigma }_{\mathrm{bkg}}^{2}={F}_{\mathrm{corrB}}\,N\,{\sigma }_{\mathrm{ann}}^{2},\end{eqnarray} \tag{ 4 }$

where F_corrA = 11.56 and F_corrB = 1 for 2MASS images, ¹⁴ accounting for the image coadding and smoothing. For WISE, F_corrA = F_corrB depending on the aperture size (Cutri et al. 2012). The variance of individual pixels (${\sigma }_{i}^{2}$) for the WISE coadded image can be calculated directly as the quadrature sum of the pixels in the uncertainty map in the same aperture.

We only consider the confusion noise (σ_conf) for WISE images because their coadded images are very deep but poorly resolved (>6''). We randomly sample the background in the image with the same aperture size as that used to measure the target (Figure 6), estimating σ_conf as the standard deviation of the sampled median flux. To reduce the correlation of the sampling, we require that the samples overlap with each other by less than 20% and contain fewer than 20% of the masked pixels. Finally, the mock tests in Section 3.2 show that our aperture photometry can measure at least 80% of the total flux of the galaxy. Based on the results of these mock tests, we adopt σ_sys = 10%–20% as the systematic uncertainty of each band. As discussed in Section 3.2, the 10%–20% flux loss is minor when compared to the typical uncertainties of stellar mass and SFR, particularly for faint galaxies that have even larger uncertainties in their physical properties. It suffices to include a systematic uncertainty in our final flux measurements to account for the scatter and offset in the flux recovery rate.

The full error budget of 2D image decomposition described in Section 3.3 is difficult to estimate. Apart from formal uncertainties, systematic uncertainties can arise from the influence of the S/N and the PSF (Guo et al. 2009; Yoon et al. 2011), especially when bright nuclei are involved (e.g., Kim et al. 2008; Zhuang & Ho 2022), and the background estimation (Huang et al. 2013; Gao & Ho 2017). Because bright AGNs are rare in our sample, we focus on the uncertainties induced by background estimation. We approach this problem empirically, by measuring the parameters for the blended system with different background levels, and calculating the scatter. For each band, we generate an idealized galaxy model using the best-fit parameters given by the 2D multiband decomposition. After adding the Poisson noise of the source and the Gaussian noise of the background to mimic a realistic mock image of a blended, interacting system, we insert the mock image in different locations on a simulated background image specifically designed for each band. The simulated background images for SDSS, 2MASS, and WISE were created by mosaicing relatively empty sky regions. The background for GALEX is more complicated because its depth depends on the diverse exposure times of the different image tiles. For galaxies that are observed in the MIS or DIS, we simulate the background image of GALEX with the Python package make_noise_image, assuming that the large-scale, random noise of the background is dominated by Gaussian noise that can be estimated from the standard deviation of the pixel values in the galaxy-size apertures calculated from the source-masked regions of the FUV and NUV images. The background of the AIS data is primarily influenced by Poisson noise; we combine the real AIS background areas in the same tile to create GALEX background images. We repeat the decomposition multiple times, during each iteration placing the target in a different location in the background. The standard deviation of the multiple trials is used as the final uncertainty.

As discussed throughout the paper, we adopt different methods to measure the flux (F) and uncertainty (σ) of the galaxy in different bands. We follow the same method to estimate the flux upper limits. We consider a source undetected if its flux measurement has S/N < 2. As in Cutri et al. (2012), we estimate the upper limit as 2σ + F, setting F = 0 if the measured F is negative. The only exception comes from the Herschel data. If no measurement is provided by the HELP catalog, we adopt 2σ₂₅₀ as the 250 μm upper limit (Section 2.1; Viero et al. 2014).

We analyze the photometric measurements of our galaxies with the widely used CIGALE code (Boquien et al. 2019). CIGALE models a panchromatic SED from the X-rays to the radio with a flexible collection of model components, including a stellar continuum, nebular emission, dust attenuation and emission, and an AGN continuum (Noll et al. 2009; Ciesla et al. 2015; Boquien et al. 2019; Yang et al. 2020, 2022). The user specifies the model components and a grid of parameters as prior information. The code calculates the likelihood, ${ \mathcal L }=\exp (-{\chi }^{2}/2)$, for each model by comparing the model fluxes with the observed fluxes. CIGALE uses a Bayesian approach to calculate the marginalized probability distribution function using the ${ \mathcal L }$ values of all models. Based on the probability distribution function, CIGALE provides the probability-weighted mean and standard deviation of physical parameters such as stellar mass (M_*), SFR, dust mass (M_d ), and AGN fraction (f_AGN). CIGALE adopts an energy conservation algorithm (Burgarella et al. 2005; Boquien et al. 2019), which requires that the UV and optical emission attenuated by dust reradiates in the MIR and FIR. In accordance with the tutorial provided by CIGALE, ¹⁵ we set the input data and errors as (2σ + F)/2 ± (2σ + F)/2. This approach allows us to include all available models conservatively.

Table 4 summarizes the parameters of the SED model. The stellar emission is represented by simple stellar populations from Bruzual & Charlot (2003) of metallicity 0.004 to 0.02 (Brown et al. 2014; Mountrichas et al. 2021), with the metal-poor end suitable for the low-mass galaxies in the sample. We adopt the Chabrier (2003) stellar initial mass function and a double-decaying exponential function to describe the star formation history. This model reproduces the SEDs of both star-forming and quenched galaxies with a modest number of free parameters (Ciesla et al. 2015, 2016; Salim et al. 2016). We mainly follow Salim et al. (2016) to choose the parameter ranges of the model. The first decaying exponential function reflects long-term star formation of the primary stellar component of the galaxy. It has an e-folding time ranging from τ_main = 600 Myr to 18 Gyr, and the age of the main stellar population is set to t_main = 12 Gyr. The second exponential function depicts the most recent burst of star formation, which can be achieved by setting τ_burst = 20 Gyr. The age of the starburst varies from t_burst = 10 Myr to 5 Gyr, and the fraction of the stellar mass formed in the starburst spans f_burst = 0 to 0.15.

Table 4.  CIGALE Model Parameters for SED Fitting

Parameter	Values
Simple stellar population: Bruzual & Charlot (2003)
Initial mass function	Chabrier
Metallicity	0.004, 0.02
Star formation history: double exponential function
e-folding time of the main stellar population: τmain (Myr)	600, 1000, 1500, 2000, 3000, 4500, 6000, 8000, 12,000, 18,000
e-folding time of the latest starburst: τburst (Gyr)	20
Fraction of stellar mass formed in starburst: fburst	0, 0.001, 0.003, 0.005, 0.01, 0.03, 0.05, 0.1, 0.15
Age of the main stellar population: tmain (Gyr)	12
Age of the latest starburst: tburst (Myr)	10, 100, 500, 1000, 1500, 2000, 3000, 4000, 5000
Nebular emission
Ionization parameter: $\mathrm{log}U$	−3.4
Fraction of Lyman continuum photons absorbed by dust	0.3
Dust attenuation: Calzetti et al. (2000)
Color excess of the nebular lines: E(B − V) (mag)	0.01, 0.15, 0.3, 0.45, 0.6
Reduction factor to apply	0.25, 0.75
Amplitude of the UV bump (Milky Way)	3.0
Slope of the power-law modification: δ	−1.2, −0.8, −0.4, 0
Dust emission: Draine et al. (2014)
Mass fraction of PAHs: qPAH (%)	2.5
Minimum radiation field: ${U}_{\min }$	0.1, 0.5, 2.5, 10, 25
Power-law slope of the dust radiation: α	2
Fraction illuminated from ${U}_{\min }$ to ${U}_{\max }$: γ	0.001, 0.006, 0.05, 0.2
AGN emission: Fritz et al. (2006)
Ratio between outer and inner radius of the torus: ${R}_{\max }/{R}_{\min }$	60
Optical depth at 9.7 μm: τ9.7	1
Slope of the radial coordinate: β	−0.5
Angle between equatorial axis and line of sight: ψ (deg)	0 (type 2), 60 (type 1)
Full opening angle of the torus: θ (deg)	100
Contribution of the AGN to the total LIR: fAGN	0, 0.001, 0.02, 0.1, 0.3, 0.6

Download table as:  ASCIITypeset image

For the nebular emission component, we set the ionization parameter to $\mathrm{log}\,U=-3.4$ and the fraction of Lyman continuum photons absorbed by dust to 0.3, which have been found to match the observed emission-line equivalent widths of SDSS spectra (e.g., Inoue 2001). The color excess of the nebular lines for both the young and old populations ranges from 0.01 to 0.6, with reduction factors of 0.25 or 0.75 compared to the attenuation of the stellar continuum. As in Brown et al. (2014), we use the starburst attenuation curve of Calzetti et al. (2000), modified by a power-law term with exponent δ = −1.2 to 0 to steepen it (Salim et al. 2018). The impact of the different attenuation laws will be discussed later. The amplitude of the UV bump (Fitzpatrick & Massa 1986) is fixed to the Milky Way value of 3.

The MIR and FIR emission of star-forming galaxies is based on the dust emission model of Draine & Li (2007), as updated by Draine et al. (2014). Draine et al. (2014) constrain the dust composition and size distribution based on Spitzer observations. The emission templates are calculated assuming that the dust is exposed to two different environments: (1) a diffuse radiation field described by a constant minimum radiation field intensity U_min; and (2) a photodissociation region with the radiation field intensity ranging from U_min to U_max. The probability distribution of the dust mass in the radiation field is dU/dM ∝ U^α , and the fraction of dust mass in the photodissociation region is γ. Following Draine et al. (2014), we fix ${U}_{\max }={10}^{7}$ and α = 2, which have a minimal effect on the SED fit (Draine et al. 2007; Aniano et al. 2012), and set ${U}_{\min }=0.1\mbox{--}25$ and γ = 0.001–0.2. The dust composition and the mass fraction of the dust in polycyclic aromatic hydrocarbons (PAHs) is given by the parameter q_PAH. We set q_PAH = 2.5% because our broadband SED cannot constrain it (see also Shangguan et al. 2018).

We include the AGN module of Fritz et al. (2006) for the AGNs and composites in the sample. The module consists of a power-law continuum from the accretion disk in the UV/optical and thermal radiation from the dusty torus in the NIR and MIR. Many parameters of the model cannot be fully constrained by the fit. Similar to Ciesla et al. (2015), we fix most of the parameters to commonly adopted values, varying only the amplitude because our main goal is to extract the SFR, M_*, and M_d for the host galaxy. We fix the angle between the equatorial axis and the line of sight to ψ = 0° for the type 2 AGNs and composites so that the torus is edge-on and there is no contribution to the UV/optical bands from the AGN; for the type 1 sources, we choose ψ = 60°, which is consistent with the distribution of torus inclination angles obtained for nearby, optically bright type 1 AGNs (Zhuang et al. 2018). We fix the radial and angular dust distribution with β = −0.5 and γ = 0, as typically found by Fritz et al. (2006). The optical depth at 9.7 μm is set to τ_9.7 = 1, because it cannot be constrained by the broad SED. We only vary the amplitude of the dust torus model, such that the fractional contribution of the AGN emission to the total IR luminosity, ${f}_{\mathrm{AGN}}={L}_{\mathrm{IR}}^{\mathrm{AGN}}/{L}_{\mathrm{IR}}^{\mathrm{total}}=0\mbox{--}0.6$, where ${L}_{\mathrm{IR}}^{\mathrm{AGN}}$ and ${L}_{\mathrm{IR}}^{\mathrm{total}}$ are integrated from 1 to 1000 μm.

The targets can be well fitted by CIGALE, achieving ${\chi }_{\nu }^{2}\lt 2$ for ≳97% of the cases (Table 5). This fraction rises to ∼99% if the blended merger pairs are excluded, owing to the relatively smaller flux uncertainties derived by the GALFITM modeling. The derived physical properties, as well as the ${\chi }_{\nu }^{2}$ for the sample, are presented in Table 5. Figure 7 shows SED fits for a star-forming galaxy, a typical active galaxy with a low level of AGN activity (f_AGN = 7.1%), and an exceptionally strong active galaxy (f_AGN = 26.1%).

Zoom In Zoom Out Reset image size

Table 5. Physical Properties of the Sample

ID	${\chi }_{\nu }^{2}$	$\mathrm{log}\,\mathrm{SFR}$	$\mathrm{log}\,{M}_{* }$	$\mathrm{log}\,{M}_{d}$	fAGN	$\mathrm{log}\,{L}_{{\rm{O}}{\rm\small{III}}}$
		(M⊙ yr−1)	(M⊙)	(M⊙)	(%)	(erg s−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
1	2.35	−1.70 ± 0.28	10.61 ± 0.03	6.76 ± 0.99	0	38.94 ± 0.20
2	3.72	−0.92 ± 0.19	10.84 ± 0.03	6.53 ± 0.47	0	38.97 ± 0.34
3	0.76	−0.10 ± 0.12	11.26 ± 0.03	8.21 ± 0.29	0	39.53 ± 0.07
4	0.44	−2.23 ± 0.78	10.37 ± 0.03	6.41 ± 0.60	0	38.67 ± 0.33
5	0.52	−1.43 ± 0.32	10.27 ± 0.04	6.23 ± 0.53	0	39.05 ± 0.14
6	0.69	−1.90 ± 0.55	10.28 ± 0.03	6.75 ± 0.64	0	39.15 ± 0.11
7	0.45	−1.76 ± 0.54	10.54 ± 0.04	6.42 ± 0.53	0	38.76 ± 0.26
8	0.68	−1.10 ± 0.22	10.63 ± 0.05	6.59 ± 0.54	0	39.35 ± 0.17
9	0.74	−1.42 ± 0.39	10.37 ± 0.05	6.74 ± 1.17	0	39.03 ± 0.20
10	0.39	−1.59 ± 0.62	10.59 ± 0.02	7.56 ± 0.37	0	38.56 ± 1.00

Note. Column (1): index number. Column (2): reduced χ² of SED fit. Column (3): SFR. Column (4): stellar mass. Column (5): dust mass. Column (6): fractional contribution of the AGN emission to the total IR luminosity; f_AGN = 0 if the galaxy is inactive (see Section 4 for details). Column (7): luminosity of [O iii] λ5007, not corrected for extinction.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Figure 8 quantifies the deviation between the data and the best-fitting SED model for each of the bands, for clarity focusing only on sources detected with S/N > 2 in each band. The data and model agree well for most of the bands across the full range of stellar masses (M_* ≈ 10¹⁰–10^11.2 M_⊙). The scatter is relatively larger in H, K_s , and W4 because they have larger measurement uncertainties. Only the 500 μm band shows substantial systematic deviation as a function of M_*. This is likely related to the "500 μm excess" (Galliano et al. 2003; Boselli et al. 2012; Smith et al. 2012; Ciesla et al. 2014; Nersesian et al. 2019). The Draine et al. (2014) model used in our fitting corresponds to an effective graybody emissivity index of β ≈ 2, but a flatter index of β = 1.5 or even less is necessary to reproduce the flatter slopes of the FIR SED, especially in low-metallicity galaxies (Boselli et al. 2012; Smith et al. 2012; Lamperti et al. 2019; Smith et al. 2019). This may explain the tendency for the lower-mass galaxies in our sample to exhibit systematically stronger 500 μm emission.

Zoom In Zoom Out Reset image size

Our sample has a median f_AGN = 10% for type 1 AGNs, f_AGN = 4% for type 2 AGNs, and f_AGN = 2% for composites. To avoid the uncertainty of applying the AGN model to objects with a very low AGN contribution (Ciesla et al. 2015), we ultimately exclude the AGN component from the fits when f_AGN < 10%. The AGN component in such weakly active systems can be overestimated by up to a factor of two from degeneracy with low dust emission due to moderate star formation. For the overall sample, galaxies with f_AGN ≥ 10% comprise 59% of the type 1 AGNs, 22% of the type 2 AGNs, and 10% of the composites; the respective median values of f_AGN for these classes are 15%, 13%, and 13%.

The main goal of our SED analysis is to derive the SFR, M_*, and M_d for the sample galaxies. Emissions from an AGN's dusty torus may contaminate the NIR and MIR portions of the SED. To evaluate the magnitude of this effect, we compare the SED fitting results of the AGN and composite targets, with and without the AGN torus module included (see Section 4 for details). Among the galaxies with f_AGN ≥ 10%, the mean differences of their SFR, M_*, and M_d before and after adding the AGN module are 0.09, 0.01, and −0.01 dex, respectively, with a maximum deviation of 0.5, 0.13, and −0.6 dex. In contrast, for the galaxies with f_AGN < 10%, we find little difference in the resulting SFRs (<0.03 dex), M_* (<0.01 dex), and M_d (<0.02 dex) compared to the uncertainties. This is not surprising given the weakness of the AGNs in our sample. To illustrate this point, we estimate the bolometric luminosities of the AGNs and composites from the observed [O iii] λ5007 luminosity, using the line measurements from the MPA–JHU catalog for the main sample and from the catalog of Liu et al. (2019) for the broad-line AGNs. We adopt a bolometric correction of 3500 as recommended by Heckman et al. (2004). Other [O iii]-based bolometric corrections can be contemplated (see the discussion in Kong & Ho 2018), but they do not alter our basic conclusion: the type 2 AGNs and composite sources in our sample have L_bol ≈ 10⁴⁰–10⁴⁴ erg s⁻¹, which are much lower than the threshold for quasars (L_bol ≥ 10⁴⁵ erg s⁻¹; Reyes et al. 2008). Only the handful of supplementary type 1 AGNs have L_bol ≈ 10⁴⁴–10⁴⁵ erg s⁻¹.

Last, but not the least, the optical emission-line diagnostics may miss heavily obscured AGNs, which, nevertheless, may be identifiable by their MIR colors (e.g., W1 − W2 > 0.8 mag; Stern et al. 2012). This color criterion is met by 23 galaxies in our sample, although their optical emission-line ratios classify them as inactive. We fit the SEDs of these 23 sources with and without the AGN module of CIGALE. Most show negligible differences.

We conclude that black hole accretion is energetically negligible in galaxies with f_AGN < 10% and that it has essentially no impact on the physical parameters of primary interest to this study (SFR, M_*, and M_d ) derived from the SED fitting.

To investigate the effects of nondetections (upper limits) in the different bands on our results, we fit mock SEDs generated from the SEDs of real targets and study the deviation of the derived physical parameters from the true input values as the number of upper limits increases. We select as the test sample galaxies detected in all the bands and take their best-fit parameters as the true input values. We randomly scale down the amplitude of the observed SED so that their r-band flux densities fall between 0.15 and 2 mJy, while preserving the observed uncertainty of each band. The flux density becomes an upper limit if it is below two times the uncertainty. This scaling guarantees that the mock SED does not have upper limits in u, g, r, i, z, and W1, as is the case in most of our observed data. In total, we randomly generate ∼2000 mock SEDs with nondetections in the FUV, NUV, J, H, K_s , W2, W3, W4, 250 μm, 350 μm, and 500 μm bands, closely mimicking the actual observations. We require that there be at least one detection in the other wavelengths to isolate the effect of nondetection on the bands of interest, since the flux densities of an SED scale together. We also create 100 mock SEDs with all bands detected after rescaling the flux, to serve as a control sample. We then use CIGALE to fit the mock SEDs with the same model discussed in Section 4.1. The input parameters are very well recovered if there are no nondetections in the mock SEDs in the control sample.

Figures 9 and 10 summarize the deviation (output − input) of M_*, SFR, and M_d from our tests. These physical parameters are not significantly affected by the 2MASS bands, and for the sake of clarity we do not include their results in the plots. The stellar mass is very well recovered (${\rm{\Delta }}\,\mathrm{log}\,{M}_{* }\approx 0.03\,$dex) even when the GALEX, Herschel, or even WISE bands are upper limits (Figure 9(a)). This is because all galaxies are detected in the five SDSS bands, which play the key role in constraining the stellar mass from the SED fitting. The SFR can be estimated fairly accurately using upper limits in the GALEX or Herschel bands individually (${\rm{\Delta }}\,\mathrm{log}\,\ \mathrm{SFR}\approx 0.03\,$dex), which is small compared to the typical uncertainty of ∼0.18 dex. SFR is moderately underestimated (0.2 dex) if all of the MIR and FIR data are upper limits. These mock galaxies typically have relatively low input SFRs (≲0.3 M_⊙ yr⁻¹), and the uncertainties of their output SFRs are large (∼0.6 dex), such that CIGALE derives SFRs consistent with the inputs, although the overall distribution is slightly underestimated. This situation affects ∼400 objects in the actual sample, and for these objects we apply a systematic correction to their SFRs according to this test.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to understand the impact of the upper limits on M_* and SFR, we use the mock SEDs to simulate the consequences of dropping the upper limits altogether, successively testing the effect of each band. Dropping the upper limits in the UV or IR bands makes little difference to M_* (Figure 9(a), bottom panel). The situation, however, differs for the SFRs, which can be significantly overestimated when we neglect the upper limits (Figure 9(b), bottom panel). All else being equal, the SFRs are better constrained by including upper limits than not. Even though the median bias is moderate, ∼0.1 dex if GALEX or Herschel is omitted and ∼0.3 dex even if all IR (2MASS, WISE, and Herschel) bands are excluded, neglecting upper limits leaves a long tail of positive residuals: ∼0.7 dex for GALEX, ≳1 dex for Herschel, and ≳1.5 dex for all IR bands. Our results agree with those of Noll et al. (2009) and Ciesla et al. (2015), who concluded that SFRs in star-forming and AGN host galaxies can be biased too high by more than 0.8 dex when IR filters are not used, on account of overestimation of the FUV and visual attenuation factors. Many works (e.g., Salim et al. 2007; Hoversten et al. 2011; Jeong et al. 2012; Brown et al. 2014; Salim et al. 2016) derive SFRs from panchromatic (UV to FIR) SEDs, but undetected bands often are not included because upper limits are usually not available. For instance, the SED fitting of Salim et al. (2016) only utilizes data from GALEX and SDSS when W3 and W4 are not detected. The bottom panel of Figure 9(b), which mimics the assumptions of Salim et al. (2016), illustrates the degree to which SFRs can be overestimated, and that the biases can be mitigated by incorporating upper limits into the analysis.

Upper limits in GALEX (Figure 10(a)) or WISE (Figure 10(b)) have little influence on the accuracy of the dust masses from SED fitting, compared to the typical 0.42 dex uncertainty. The dust mass can be reasonably well constrained if there is a detection in at least one Herschel/SPIRE band (Figure 10(c)). However, if all three SPIRE bands are upper limits, the resulting dust masses are ∼0.24 dex lower than the input values. In contrast, removing all the SPIRE data significantly overpredicts M_d by ∼1 dex or more (red dashed line in Figure 10(c)). This is not unexpected. The FIR data are needed to constrain the peak of the cold dust SED (Draine et al. 2007), which not only helps to derive the correct FIR luminosity but also to prevent unphysically large attenuation through the energy-balance constraint. For galaxies that are not detected in any of the Herschel bands, using our mock tests as a guide, we apply a systematic correction of 0.24 dex to their dust masses, although we stress that the correction is small compared to the statistical uncertainty of M_d for individual galaxies (0.42 dex).

We note that our use of galaxies whose SEDs have been detected in all bands to mimic galaxies with nondetections in some bands may not accurately represent all galaxy types. For example, quiescent galaxies may have systematically different SEDs compared to the mock sample. However, the main purpose of this series of mock tests is to demonstrate the importance of including upper limits in SED fitting instead of simply disregarding them. This point remains valid regardless of possible sample selection effects. While the empirical corrections derived from our mock tests may not be accurate if the SED of the measured galaxy differs significantly from that of our SED-detected sample, this is likely a higher-order effect, which is beyond the scope of this work. We advise caution when using dust mass and SFR corrections for individual quiescent galaxies.

As mentioned in Section 4.1, we adopt Bayesian estimates of M_*, SFR, and M_d based on the likelihood of each SED model. This method takes into account the degeneracy of parameters among different models if they provide indistinguishable likelihoods. The uncertainties of the likelihood-weighted physical properties are calculated with consideration of the upper limits in the data (Boquien et al. 2019). Our choice of parameter ranges covers typical values of nearby galaxies whose SEDs have been well fit (e.g., Brown et al. 2014; Ciesla et al. 2015). Moreover, we test the systematics and scatter introduced by the upper limits; manually converting detections to upper limits (Section 4.4) demonstrates that our physical parameters of interest are minimally affected by the upper limits. Our optical photometry from SDSS delivers robust M_* for all of the targets, with typical uncertainties of 0.05 dex. The simulated values of SFR and dust mass are consistent within 2σ, with 92% and 99% of the true values, respectively. We show that including the upper limits in SED fitting reduces the bias of SFR and dust mass measurements. Although the derived quantities are usually consistent with the input values, considering the uncertainties, the sample-averaged SFRs and dust masses can be underestimated at the level of ∼0.2 dex when the target is not detected in the WISE and Herschel bands (Section 4.4). We correct for this source of systematic bias for the faint sources. These galaxies usually have SFR and/or dust mass uncertainties >30% of the best-fit values from the SED fitting. We note that the sample distributions of SFR and M_d may be biased in regions of low values.

We stress, however, that our tests do not account for uncertainties stemming from prior model assumptions. For example, at a fixed stellar initial mass function, the stellar population model and star formation history can introduce a 0.3 dex uncertainty to the stellar mass and SFR (Conroy 2013). It is difficult to push the sensitivity of the specific SFR (sSFR) derived from integrated SEDs down to values ≲10⁻¹² yr⁻¹ because of uncertainties associated with the contribution of UV emission from evolved stars and dust heating (O'Connell 1999; Conroy 2013). Nearly 30% of the galaxies in our sample have an sSFR close to or below 10⁻¹² yr⁻¹. One should regard their SFRs with some caution.

The stellar masses from our work agree well with those of GSWLC-2 with a negligible zero-point offset and scatter (median and standard deviation ${\rm{\Delta }}\mathrm{log}{M}_{* }=0.01\pm 0.11$ dex; Figure 11(a)). The 0.1 dex scatter is consistent with the measurement dispersion of M_* induced by nondetections in the UV or other bands (Section 4.4). In comparison with the values in the MPA–JHU catalog, our stellar masses are on average slightly larger (${\rm{\Delta }}\mathrm{log}{M}_{* }=0.11\pm 0.11$ dex; Figure 11(b)), but we note that a systematic trend is visible in the residuals, such that toward the low-mass end our masses tend to be higher. This discrepancy likely arises from the difference in the star formation history used (Salim et al. 2016). Adopting the delayed exponential star formation history assumed in the MPA–JHU catalog slightly reduces the systematic differences between their stellar masses and ours (${\rm{\Delta }}\mathrm{log}{M}_{* }=0.08\pm 0.12$ dex; Figure 11(c)), but cannot totally remove the trend. The residual offset of 0.08 dex owes to to the fact that our optical fluxes are slightly brighter (e.g., 0.08 mag in u and 0.05 mag in g) than the modelmag values used in the MPA–JHU catalog (Appendix A; Figure A2). This effect influences less-massive galaxies more significantly because they have a higher fraction of young stars (Sextl et al. 2023), resulting in a larger stellar mass offset toward lower mass. The inclusion of UV photometry, as implemented in GSWLC-2, significantly improves the stellar mass measurements, which otherwise may be biased by only fitting the SDSS photometry, particularly for galaxies with a young stellar population.

Zoom In Zoom Out Reset image size

The heavily blended, merging galaxies in our sample deserve closer scrutiny. For these sources, which comprise 7% of the sample, our stellar masses show minor systematic differences compared to the values given in GSWLC-2 (${\rm{\Delta }}\mathrm{log}{M}_{* }=0.03\pm 0.19$ dex), and MPA–JHU (${\rm{\Delta }}\mathrm{log}{M}_{* }=0.11\pm 0.23$ dex). This is not surprising because the stellar mass estimates are mainly driven by the optical photometry based on SDSS images (Section 4.4), and compared to the uncertainties, there are no significant statistical differences between the photometry derived from our 2D deblending (Section 3.3) and the SDSS model magnitudes used in other catalogs.

For the sample as a whole, our SFRs show a relatively large scatter and modest zero-point offset with respect to GSWLC-2 (ΔlogSFR = −0.12 ± 0.56 dex; Figure 12(a)). Closer inspection reveals that most of the discrepancy is confined to the 45% of objects that have low star formation activity: while ΔlogSFR = −0.05 ± 0.23 dex at SFR > 0.3 M_⊙ yr⁻¹, ΔlogSFR = −0.40 ± 0.71 dex at SFR ≲ 0.3 M_⊙ yr⁻¹. Now, distinct from our work, which strives to provide full spectral coverage from the FUV to FIR for the SED fits, GSWLC-2 did not incorporate the J, H, K_s , W1, and W2 bands into their analysis. To evaluate the extent to which this difference in spectral coverage contributes to the discrepancies between the two studies, we repeat our SED fits without these NIR bands. The resulting SFRs are not strongly affected, for the difference between the results for the full and partial SEDs is only ΔlogSFR = 0.002 ± 0.06 dex. We suspect that the systematic deviation between our results and those of GSWLC-2 stems from the manner in which the two studies treat the SED longward of ∼20 μm. On the one hand, we rigorously include in the SED fit all (14) bands from the FUV to 500 μm, using new, customized, more accurate photometric measurements and their associated uncertainties. Upper limits are rigorously included in our sample. Salim et al. (2018) did not have direct access to FIR data and instead extrapolated the total IR emission from MIR flux densities. On the other hand, the unWISE photometry they used for the W4 or W3 (if W4 is unavailable) band is ∼0.3 mag brighter than the profile-fitting photometry we use for the unresolved galaxies that dominate the low-SFR population. For the resolved objects, the unWISE W4 and W3 photometry is ∼0.2 mag brighter than our aperture photometry, as explained in Appendix A (Figure A5). GSWLC-2 predicted the total IR luminosity from unWISE W3 photometry for >40% and W4 photometry for >20% of the low-SFR galaxies, which accounts, at least in part, for the systematic difference in SFR. Further comparisons are not straightforward because GSWLC-2 amalgamated photometry is derived using multiple methods. For instance, the unWISE photometry was measured by model fitting (Lang et al. 2016b), which differs from the method used for treating the UV bands. Among the population of galaxies with SFR ≲ 0.3 M_⊙ yr⁻¹ for which our two studies show the most glaring inconsistency, ∼36% of the sources are not detected in W3 or W4 by unWISE. Without the MIR data, the SFRs of GSWLC-2 default to the UV/optical SED fitting from Salim et al. (2016), which likely overestimated the SFRs. To demonstrate this point, in Section 4.4, we fit the mock SEDs with nondetections in WISE, while intentionally dropping all of the data at wavelengths longer than the SDSS z band. Under these circumstances, which mimic the procedure of Salim et al. (2016), we confirm that the median value of the SFR indeed gets systematically overestimated, even up to ∼1.5 dex. Therefore, we caution that the SFRs from GSWLC-2 may be significantly overestimated in the low-SFR regime (SFR ≲ 0.3 M_⊙ yr⁻¹). Section 4.4 emphasizes the importance of properly including upper limits for the WISE bands and bands at longer wavelengths to avoid overestimating the SFR.

Zoom In Zoom Out Reset image size

An even more complex situation emerges when we contrast our measurements with those from the MPA–JHU database. First, recall that the total SFRs from the MPA–JHU catalog come from the combination of the contribution measured from the fiber spectrum plus that obtained from SED fitting of broadband UV and optical photometry outside the fiber. This method gives SFRs of considerably large uncertainty. As shown in Figure 12(b), sources with SFR ≳ 0.5 M_⊙ yr⁻¹ track our measurements fairly well, with ΔlogSFR = −0.00 ± 0.29 dex. The scatter is comparable to the uncertainties of the SFRs in the MPA–JHU catalog in this SFR range. At the highest SFRs, the MPA–JHU values have a mild tendency to be higher than ours. This trend can be explained by our choice of dust attenuation. Although we can achieve better agreement between the two sets of measurements by using the original attenuation law of Calzetti et al. (2000), a modified attenuation law with a steeper slope has been found to be preferable by many investigators (e.g., Buat et al. 2011; Salim et al. 2018; Qin et al. 2022). The departures between our measurements and those of the MPA–JHU catalog become alarmingly severe (ΔlogSFR = −0.47 ± 0.71 dex) once SFR ≲ 0.2 M_⊙ yr⁻¹. In the range SFR ≈ 0.001–0.1 M_⊙ yr⁻¹, the SFRs in the MPA–JHU catalog remain more or less fixed at SFR ≈ 0.1 M_⊙ yr⁻¹. The low-SFR sources also have vastly larger uncertainties in the MPA–JHU catalog (median 0.70 dex) than in ours (median 0.45 dex). Within the MPA–JHU catalog, the fiber SFRs for AGNs, composites, and unclassifiable objects (with low-S/N spectra) derive from a relation between the 4000 Å break (D4000) and the sSFR calibrated from emission-line galaxies for which both quantities could be measured. The objects in our sample with low SFRs and large offsets belong to this category. As shown in Figure 11 of Brinchmann et al. (2004), when D4000 is large (D4000 ≳ 1.6), it starts to become less sensitive to sSFR and exhibits a large spread. In particular, the mean relation exhibits an unphysical rise in sSFR at D4000 > 2. The large uncertainty and potential bias of the D4000–sSFR relation may boost the SFR in quiescent galaxies and account for the systematical underestimation of SFRs observed at SFR < 0.2 M_⊙ yr⁻¹.

Lastly, we again examine the subset of heavily blended galaxies that we analyzed by 2D decomposition (Section 3.3). The deviation of SFRs for the blended components is slightly larger than that for isolated objects. Our measurements of the SFRs of the individual galaxies in the blended pairs are 0.14 ± 0.69 dex lower than those in the GSWLC-2 catalog and 0.19 ± 0.60 dex lower than those in the MPA–JHU catalog.