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

The broadband spectral energy distribution of a galaxy encodes valuable information on its stellar mass, star formation rate (SFR), dust content, and possible fractional energy contribution from nonstellar sources. We present a comprehensive catalog of panchromatic photometry, covering 17 bands from the far-ultraviolet to 500 $\mu$m, for 2685 low-redshift (z=0.01-0.11), massive ($M_*>10^{10}\,M_\odot$) galaxies selected from the Stripe 82 region of the Sloan Digital Sky Survey, one of the largest areas with relatively deep, uniform observations over a wide range of wavelengths. Taking advantage of the deep optical coadded images, we develop a hybrid approach for matched-aperture photometry of the multi-band data. We derive robust uncertainties and upper limits for undetected galaxies, deblend interacting/merging galaxies and sources in crowded regions, and treat contamination by foreground stars. We perform spectral energy distribution fitting to derive the stellar mass, SFR, and dust mass, critically assessing the influence of flux upper limits for undetected photometric bands and applying corrections for systematic uncertainties based on extensive mock tests. Comparison of our measurements with those of commonly used published catalogs reveals good agreement for the stellar masses. While the SFRs of galaxies on the star-forming main sequence show reasonable consistency, galaxies in and below the green valley show considerable disagreement between different sets of measurements. Our analysis suggests that one should incorporate the most accurate and inclusive photometry into the spectral energy distribution analysis, and that care should be exercised in interpreting the SFRs of galaxies with moderate to weak star formation activity.


INTRODUCTION
The galaxy population displays a bewildering variety of morphology, internal structure, kinematics, gas content, star formation rate and efficiency, and level of central black hole accretion. At the same time, the diversity of galaxy properties is interwoven by a number of empirical scaling relations that link them to each other, and many physical parameters are strongly coupled to the galaxy mass (e.g., Kennicutt & Evans 2012;Kormendy & Ho 2013;Cappellari 2016;Saintonge & Catinella 2022). Environments also play an important role (e.g., , as does, of course, cosmic epoch (e.g., shangguan@mpe.mpg.de Shapley 2011;Madau & Dickinson 2014). Within this backdrop, low-redshift galaxies serve as a valuable laboratory to explore the myriad outstanding complexities of galaxy evolution. Apart from providing witness to the latest stage in the cosmic lifecycle of galaxies, the local Universe offers the obvious advantage of proximity. Nearby galaxies are bright and well-resolved. In the last two decades, the advent of large-area sky surveys has furnished a rich inventory of multiwavelength data suitable for probing the statistical physical properties of nearby galaxies, whose panchromatic coverage includes, among others, the ultraviolet (UV) by the Galaxy Evolution Explorer (GALEX; Martin et al. 2005), the optical by the Sloan Digital Sky Survey (SDSS; York et al. 2000), the near-infrared (NIR) by the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), the mid-arXiv:2307.13461v1 [astro-ph.GA] 25 Jul 2023 infrared (MIR) by the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010, and, over much more limited but still sizable areas, the far-infrared (FIR) by Herschel (Pilbratt et al. 2010). Despite their relatively shallow depth and moderate resolution, these large-area databases constitute an invaluable resource for probing, at the very minimum, several fundamental properties of nearby galaxies on global scales, through analysis of their integrated, broadband spectral energy distribution (SED). Broadband SED fitting furnishes key physical parameters such as stellar mass (M * ), star formation rate (SFR), dust mass (M d ), and even an assessment of the contribution by active galactic nuclei (AGNs) to the total luminosity (e.g., da Cunha et al. 2008;Conroy 2013;Boquien et al. 2019). Coupled with optical imaging and spectroscopy (e.g., from SDSS), one can explore the rich interplay between star formation, black hole accretion, and the interstellar medium on the one hand, and galaxy structure and environment (e.g., tidal interactions, mergers, environment) on the other hand.
The SDSS Stripe 82 region is one of the largest areas observed with relatively high and uniform sensitivity observed across a wide range of wavelengths. Scanned repeatedly by the SDSS survey for ∼ 80 times, the coadded optical images are about 2 magnitudes deeper than the SDSS legacy survey (Abazajian et al. 2009;Jiang et al. 2014), rendering them especially sensitive to lowsurface brightness features in galaxy outskirts (e.g., Fliri & Trujillo 2016). Apart from the standard all-sky coverage provided by GALEX, 2MASS, and WISE, numerous dedicated multiwavelength surveys have been conducted in the Stripe 82 region, among them the VISTA-CFHT survey in the NIR (Geach et al. 2017), the Spitzer-IRAC Equatorial Survey in the MIR (Timlin et al. 2016), high-resolution imaging with the Very Large Array in the radio (Hodge et al. 2011), and Chandra observations in the X-rays (LaMassa et al. 2016). Of most relevance for securing wavelength coverage that impacts the derivation of SFRs and dust masses, a portion of Stripe 82 was observed in the FIR by the Herschel Stripe 82 Survey (HerS; Viero et al. 2014) and the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012), which makes Stripe 82 the largest area scanned by Herschel to a uniform depth (Lutz 2014;Shirley et al. 2021). These Herschel observations have much better sensitivity than the all-sky FIR surveys executed by the Infrared Astronomical Satellite (Neugebauer et al. 1984) or the AKARI Far Infrared Surveyor (Kawada et al. 2007). The above-summarized databases distinguish Stripe 82 as one of the best regions to study the comprehensive SED of a sizable population of representative low-redshift galaxies.
Stripe 82 has been the target of diverse investigations. Most previous studies of low-redshift galaxies in this survey area largely capitalized on the advantages afforded by the deep optical and NIR imaging to probe their low-surface brightness features, morphology, and local environment. Examples include the construction of a large, mass-limited sample of massive galaxies (Bundy et al. 2015), searching for tidal features in early-type systems (Kaviraj 2010;Yoon & Lim 2020), surveying for ultra-diffuse galaxies (Zaritsky et al. 2021), investigating the host galaxies of nearby AGNs (Matsuoka et al. 2014) and their incidence of mergers (Karhunen et al. 2014;Hong et al. 2015), and measuring internal galactic substructures, such as bulges and disks (Bottrell et al. 2019;Sachdeva et al. 2020), disk truncations and extended halos (Peters et al. 2017), and asymmetry (Yesuf et al. 2021). Ellison et al. (2016) used Herschel observations of a subset of Stripe 82 galaxies as a training set for artificial neural network prediction of total IR luminosity for SDSS galaxies. In a similar spirit, Rosario et al. (2016) assessed the systematics of optical SFRs with the aid of IR-based SFRs computed from Herschel and WISE observations of low-redshift Stripe 82 galaxies. Nevertheless, the valuable repository of multiwavelength data for Stripe 82 has yet to be fully exploited to characterize systematically the comprehensive (UV to FIR) SEDs of its constituent galaxies to derive their basic physical properties. Value-added databases of SDSS galaxies, such as that maintained by MPA-JHU (Kauffmann et al. 2003;Brinchmann et al. 2004;Salim et al. 2007) and the GALEX-SDSS-WISE Legacy Catalog 1 and 2 (Salim et al. 2016(Salim et al. , 2018, certainly encompass the Stripe 82 field. However, the SED analysis that underpins these catalogs did not incorporate the broadest or most inclusive wavelength coverage in the photometry used, in large part because most SDSS galaxies lack Herschel observations. A more subtle, often overlooked limitation of previous efforts is that they usually amalgamate measurements from extant primary source catalogs, which employed different methods of photometry and uncertainty estimation. This can introduce hard-to-quantify systematic uncertainties and bias the SED fitting results (e.g., Hill et al. 2011;Clark et al. 2018). Such systematic uncertainties are especially pernicious for low-redshift galaxies, as many catalogs performed photometric measurements assuming that the sources are point-like [using point-spread function (PSF) fitting photometry], while, in actuality, the galaxies are resolved or marginally resolved Calderón-Castillo et al. 2019;Miller et al. 2020). One can measure galaxy fluxes consistently across different bands by modeling its extended surface brightness distribution. For instance, Lang et al. (2016b) measured fluxes in the four WISE MIR bands (3.4 to 22 µm) using the code Tractor (Lang et al. 2016a) by constraining the galaxy with the SDSS r-band "exponential + de Vaucouleurs" model. Alternatively, the GALAPAGOS code  can fit the light distribution of the galaxy with a parametric model to yield its total flux simultaneously in multiple bands Psychogyios et al. 2020). The latest version of the code (Häußler et al. 2022) even permits simultaneous multi-band decomposition of the bulge and disk components. ProFuse (Robotham et al. 2022) offers similar functionalities. While effective and efficient, the necessarily idealized and oversimplified fits usually invoked in parametric methods might introduce systematic uncertainties depending on the degree to which the models reproduce true light distribution of the galaxy.
The complications of parametric models can be obviated by adopting a model-independent strategy for aperture photometry. The aperture size, shape, and orientation can be chosen to optimally capture the desired boundaries of the source, and they can be fixed for consistency across different bands. As an example, Wright et al. (2016) conducted forced-aperture photometry of galaxy images in 21 bands (UV to FIR) for the Galaxy and Mass Assembly sample (GAMA; Baldry et al. 2010;Driver et al. 2011Driver et al. , 2016 by using as prior an aperture defined from a high-resolution optical image and then calculating the working aperture for other bands, after convolving the prior aperture with the PSF of the corresponding bands. Clark et al. (2018) adopted a similar approach for the extensive panchromatic (UV to submillimeter) data for the DustPedia sample of nearby galaxies. Developed in the same spirit, the source detection and photometry code ProFound (Robotham et al. 2018) starts with a segmentation map created from the stacked optical and NIR images and then enlarges it to generate the source aperture according to the surface brightness depth and resolution of the image.
We present a comprehensive catalog of panchromatic photometry for 2781 low-redshift (z = 0.01−0.11) galaxies selected from the HerS region of Stripe 82, covering 17 bands from the far-UV (FUV) to 500 µm. Of these, 2668 sources belong to a primary sample of active and inactive galaxies with stellar masses M * > 10 10 M . We perform SED fitting to derive stellar masses, SFRs, dust masses, and AGN luminosity fractions, which will serve as the foundation for a series of forthcoming works. This paper describes the technical details of our approach to securing robust and physically consistent measurements of total galaxy fluxes and associated uncertainties across many bands, and, where necessary, proper upper limits.
As most of the galaxies in our sample are marginally to fully resolved from the UV to the MIR, but not in the FIR, we develop a hybrid approach of matchedaperture photometry, which can be used effectively for the 14 shorter wavelength bands (FUV to 4.6 µm), along with profile-fitting photometry for the longer wavelength bands. Interacting galaxies and sources in crowded regions, such as the cores of galaxy groups and clusters, are decomposed properly using multi-band, simultaneous two-dimensional (2D) image decomposition. Special care is devoted to estimating the full error budget, which includes uncertainties associated with background variance and contamination by foreground stars, which we treat using a new method. We design a series of mock tests to evaluate the impact of including certain critical bands and their respective upper limits in the SED analysis. Comparing our new measurements with those of previous works reveals that SED fitting-based SFRs of galaxies with moderate to low star formation activity should be used with considerable caution.
The paper is structured as follows. Section 2 introduces the sample definition and the multiwavelength data used in this work. Our methodology is introduced in Section 3, which includes our techniques for image preprocessing and deblending, aperture-matched photometry, model-fitting photometry for heavily blended galaxies, and estimating uncertainties and upper limits. Section 4 presents the method for SED fitting, the results obtained therefrom, the effect of AGNs, the influence of nondetections, and parameter uncertainties. Section 5 compares the stellar masses and SFRs derived in this work with those of other widely used catalogs. The implications of this work are discussed in Section 6, and the main conclusions are summarized in Section 7. Some technical details are relegated to Appendices A-C. We adopt a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 and Ω Λ = 0.7. Stellar masses and SFRs reported in this study are based on the Chabrier (2003) stellar initial mass function.

Sample Definition
We select galaxies in the 79 deg 2 SDSS Stripe 82 region that is covered by the Herschel Stripe 82 Survey (HerS; Viero et al. 2014) with the Herschel/SPIRE instrument. The SDSS DR16 database provides 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.  (5): Flag for the photometric method: 0 = no photometry because of serious contamination from saturated stars or projected background galaxies, optical isophotal size is unavailable, or strongly affected by distortion effects from the edge of the GALEX field; 1 = aperture photometry; 2 = photometry by 2D decomposition; 3 = aperture photometry for supplementary object; 4 = photometry by 2D decomposition for supplementary object. Col. (6): Flag to indicate whether the target is in an interacting system: 0 = isolated; 1 = merger pair system; 2 = cluster central region; 3 = uncertain (a merger system or an isolated galaxy with subtructures); 4 = merger system consisting of > 2 galaxies. Col. (7): Index number of the merger companion if interaction flag = 1; if the merger system includes more than two galaxies, the companion is the closest galaxy. Col. (8) We also require the morphological type of the targets to be GALAXY 1 when selecting targets in CasJobs 2 . Figure 1 shows for the Stripe 82 galaxies the distribution of redshift versus stellar mass (M * ) derives from the GALEX-SDSS-WISE Legacy Catalog 2 (GSWLC-2; Salim et al. 2018). A lower redshift cut of z > 0.01 is required in the GSWLC-2 to exclude the very nearby galaxies whose distance is poorly indicated by their redshift (Tully et al. 2016). The two notable large-scale structures arise from overdensities at z ≈ 0.07 and 0.13 ( Figure 1a). 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 1b color codes the data points by the detection rate (f det ) of 1 https://www.sdss.org/dr16/algorithms/classify/#Star/ GalaxyClassification 2 http://skyserver.sdss.org/CasJobs/ SPIRE 250 µm as given in the HELP catalog (Shirley et al. 2021), which considers flux densities with S/N > 2 as detections. Our sample is defined by the 2681 galaxies with M * > 10 10 M and z < 0.11, which have an average 250 µm detection rate of f det 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 mag 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 effect strongly impacts the UV data.
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 −1 ) within a projected sepa-  (Salim et al. 2018) and redshift for galaxies in Stripe 82. We code the targets by (a) the number density of SDSS galaxies in each stellar mass-redshift bin and (b) the detection rate in Herschel/SPIRE 250 µm (HELP catalog; Shirley et al. 2021). Galaxies with S/N > 2 are considered detections. We bin stellar mass in steps of 0.5 dex and redshift in steps of 0.01. The detection fraction is defined as the ratio between the number of detections over the number of galaxies in each bin. The cuts of redshift and stellar mass for our sample are shown as black dash lines (see text for more details).
ration 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 magnitudes 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 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 magnitudes 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 red- shifts 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 active galactic nuclei (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 (LIN-ERs;Heckman 1980), which are predominantly highly sub-Eddington AGNs (Ho 2008(Ho , 2009). Since such weak emission-line sources will contribute negligibly to the broadband SED, 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. Al-though these objects were excluded from the GSWLC-2 and therefore do not have stellar masses, their stellar masses are expected to exceed 10 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 11 M . In summary: the final sample has 2781 galaxies ( Table 1), consisting of 2668 objects that satisfy the original stellar mass cut (M * > 10 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.

Data
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.
• 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 = 4. 2 and 5. 3, respectively (Martin et al. 2005). Although the General Release 6 and 7 (GR6+7; Bianchi et al. 2014Bianchi et al. , 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 total galaxy fluxes using a consistent, common aperture across the SED (Section 3 The AllWISE catalog provides fluxes from PSFfitting 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 5 .
• FIR: The Herschel Stripe 82 Survey (HerS; Viero et al. 2014) observed an area of 79 deg 2 in the SDSS Stripe 82 region, covering (α, δ) J2000 = (0h 54m, −2 • ) to (2h 24m, +2 • ). The SPIRE imaging photometer (Griffin et al. 2010) was used to survey the 250, 350, and 500 µm bands to an average depth of 13.0, 12.9, and 14.8 mJy beam −1 , respectively. Although the rest of the Stripe 82 region is covered by the Herschel Multi-tiered Extragalactic Survey (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 18. 2, 25. 2, and 36. 3 at 250, 350, and 500 µm, respectively. We take advantage of the comprehensive photometry offered by the Herschel Extra-galactic Legacy Project (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 7. 4, 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 HELPunmeasured 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 3 times.

METHOD
Apart from the three Herschel FIR bands, whose fluxes, owing to limitations of spatial resolution, we simply adopt from the HELP catalog (Section 2.2), for all the other 14 bands from FUV to W4 we apply a uniform method to measure the total flux using a common aperture defined by the optical profile of the galaxy. This section describes the steps necessary to preprocess the images, including source masking and detection, background subtraction, deblending, and removal of contaminants. We derive aperture-matched photometry for galaxies that are isolated or that can be deblended adequately. A minority of merging galaxies or galaxies in exceptionally crowded environments require more sophisticated 2D parametric decomposition. The photometric measurements are summarized in Tables 2 and 3. All flux and magnitude measurements in our catalog have been corrected for Galactic extinction, based on the extinction curve outlined by Fitzpatrick & Massa (2007) and an assumed value of R V = 3.1.

Preprocessing and Deblending
We take the aperture of a galaxy (with semi-major axis a 25 and semi-minor axis b 25 ) twice as large as its isophote at 25 mag 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 25 < 50 , and to 1000 × 1000 if a 25 ≥ 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 5 connecting pixels have a flux more than 2 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, semi-major axis, semi-minor axis, and position angle (PA) of the galaxy. With this geometric information in hand, we grow the semi-major and semi-minor axes 3 times and derive the detection mask. An example is shown in Figure 3a. To determine and remove the final background, we apply the mask and fit the background pixels with a third-order 2D polynomial model.
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 6 = 32 and contrast 7 = 0.0001 to aggressively break up the original segments into smaller segments 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 3b).
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 pro- gri Figure 4. Illustration of the aperture (dashed ellipse) for a resolved galaxy in images from FUV to W4, color-coded by the flux intensity in each filter. The aperture is twice the size of the isophote of the galaxy at 25 mag arcsec −2 in the r band. Optical images come from the coaddition of SDSS Stripe 82 data produced by Fliri & Trujillo (2016). The first panel shows the SDSS gri color composite image provided by the SDSS SkyServer.
file of the target using the task ellipse in Pyraf 8 (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 average surface brightness 2 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 densely sample the surface brightness profile of the galaxy. 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 3c). 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 fluctuation, which is estimated following the methodology of Li et al. (2011). We define the boundary of the galaxy by the semi-major (a profile ) and semi-minor (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 of the total magnitude of the star in the UV or IR using a newly developed random forest 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 optica 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 2 score greater than 0.94. A small number (30) of 8 https://pypi.org/project/pyraf/ 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.

Aperture-matched Photometry
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 mag arcsec −2 . This surface brightness level is trivially satisfied by Stripe 82, which reaches a depth of µ r ≈ 28.5 mag arcsec −2 (Fliri & Trujillo 2016). We set the semi-major axis of the aperture to a 25 and its semiminor axis to b 25 if both axes are larger than 3 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.
The W3 and W4 bands require special consideration. The coarse PSFs of these two WISE bands (FWHM = 8. 4 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 tar-   (2), (3), (5), and (7): Flux density and associated uncertainty for W1, W2, W3, and W4 bands. The deblended flux densities for interacting systems are listed as separate entries. Cols. (4) and (6): Flag for extended object for W3 and W4, respectively: 0 = unresolved; 1 = resolved according to method of Appendix C; 2 = blended, interacting system. Col. (8): Flag for cross-matching SDSS coordinates to the HELP catalog: "forced" = object located within 1 from an object in the HELP forced-photometry catalog; "blind" = object located within 7. 4 from an object in the HELP blind-photometry catalog; null otherwise. Cols. (9)-(11): Flux density and uncertainty of Herschel/SPIRE 250, 350, and 500 µm bands. The value listed is the median flux density of the posterior distribution from the HELP catalog (Shirley et al. 2021 get 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 to generate mock multi-band 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 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).

Model-fitting Photometry for Heavily Blended Galaxies
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 either. We perform 2D image decomposition to deblend the galaxies using GALFITM Vika et al. 2013), which is a multiwavelength extension of GALFIT (Peng et al. 2002. 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 9 to register all the bands to the world coordinate system of the W1 image and rebin them to the lowest common pixel scale of 1. 375, 9 https://reproject.readthedocs.io/en/stable/celestial.html 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 secondorder 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 25 , b 25 , and PA 25 . The magnitude in each band is allowed to vary freely. An example fit is given in Figure 5.
The FUV and NUV bands are excluded from the simultaneous multi-band fit because they consistently give unstable results, possibly a consequence of the extreme variations in galaxy substructure due to the young stellar population 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, 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 compositecolor 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 aperturematched 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 Runge's 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 multi-band 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).

Estimating Uncertainties and Flux Upper Limits
The uncertainty of the aperture photometry is The mask instead of the image is displayed. We require that the red ellipses contain < 20% of the masked pixels.
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, σ 2 source = I/G, with G specifying the gain of the image. The background noise can be estimated from an annulus around the target aperture, σ 2 bkg = N σ 2 ann , 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 6a, 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 semi-major axis of the aperture, so that the area of the annulus is the same as that of the aperture.
As the 2MASS and WISE images are resampled and interpolated, we need to consider their pixel-correlated noise (Jarrett et al. 2000a). We adopt where F corrA = 11.56 and F corrB = 1 for 2MASS images 10 , accounting for the image coadding and smooth-10 https://wise2.ipac.caltech.edu/staff/fmasci/ApPhotUncert corr. pdf ing. For WISE, F corrA = F corrB depending on the aperture size (Cutri et al. 2012). The variance of individual pixels (σ 2 i ) 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 on 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 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;, and 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 to 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 limit. 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 σ 250 as the 250 µm upper limit (Section 2.1; Viero et al. 2014).

CIGALE Model
We analyze the photometric measurements of our galaxies with the widely used CIGALE code (Boquien et al. 2019). CIGALE models the panchromatic SED from the X-rays to the radio with a flexible collection of model components, including the stellar continuum, nebular emission, dust attenuation and emission, and the AGN continuum (Noll et al. 2009;Ciesla et al. 2015;Boquien et al. 2019;Yang et al. 2020Yang et al. , 2022. The user specifies the model components and a grid of parameters as prior information. The code calculates the likelihood, L = exp(−χ 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 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 11 , we set the input data and errors as (2σ + F )/2 ± (2σ + F )/2. This approach allows us to conservatively include all available models. 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 starforming and quenched galaxies with a modest number of free parameters (Ciesla et al. 2015(Ciesla et al. , 2016Salim et al. 2016). We mainly follow Salim et al. (2016) to choose the parameter range of the model. The first decaying exponential function reflects the long-term star formation of the primary stellar component of the galaxy. It has an e-folding time ranging from τ main = 600 Myr to , 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 Aniano et al. 2012), and set U min = 0.1−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 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, and L total IR are integrated from 1 to 1000 µm.

Results of SED Fitting
The targets can be well fitted by CIGALE, achieving χ 2 ν < 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 GALFITM modeling. The derived physical properties, as well as the χ 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%). 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 − 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 flatter slopes of the FIR SED, especially in low-metallicity galaxies (Boselli et al. 2012;Smith et al. 2012Smith et al. , 2019Lamperti et al. 2019). This may explain the tendency for the lower mass galaxies in our sample to exhibit systematically stronger 500 µm emission.

AGN Contribution
Our sample has a median f AGN = 10% for the type 1 AGNs, f AGN = 4% for the type 2 AGNs, and f AGN = 2% for the composites. To avoid the uncertainty of applying the AGN model to objects with very low AGNcontribution (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 2 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 SFR, M * , and M d for the sample galaxies. Emissions from the AGN 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 difference 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), or M d (< 0.02 dex) compared to the uncertainties. This is not surprising given the weakness of the AGNs in our sample.  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 40 − 10 44 erg s −1 , which is much lower than the threshold for quasars (L bol 10 45 erg s −1 ; Reyes et al. 2008). Only the handful of supplementary type 1 AGNs have L bol ≈ 10 44 − 10 45 erg s −1 .
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 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 pri-mary interest to this study (SFR, M * , and M d ) derived from SED fitting.

Reliability of SED Fitting with Nondetections
To investigate the effects of nondetections (upper limits) in 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 to 2 mJy, while preserving the observed uncertainty of each band. The flux density becomes an upper limit if it is below 2 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 (∆ log M * ≈ 0.03 dex) even when the GALEX, Herschel, or even WISE bands are upper limits (Figure 9a). 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 (∆ log SFR ≈ 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 −1 ), 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.
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 9a, bottom panel). The situation, however, differs for the SFRs, which can be significantly overestimated when we neglect the upper limits (Figure 9b, 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) andCiesla et al. (2015), who concluded that SFRs in star-forming The effect on dust mass is moderate for GALEX and WISE but substantial for Herschel. Combining WISE detections in the W3 and W4 bands, together with upper limits in Herschel/SPIRE, can provide useful dust masses after correcting for systematic bias.
and AGN host galaxies can be biased too high by more than 0.8 dex when the 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 data when W3 and W4 are not detected. The bottom panel of Figure 9b, 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 10a) or WISE (Figure 10b) 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 10c). 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 10c). This is not unexpected. The FIR data are needed to constrain the peak of the cold dust SED , 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 mass, 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 SED has 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 higherorder effect, which is beyond the scope of this work. We advise caution when using dust mass and SFR corrections for individual quiescent galaxies.

Parameter Uncertainties
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 uncertainty of the likelihood-weighted physical properties is calculated with consideration of the upper limit 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 upper limits; manually converting detections to upper limits (Section 4.4) demonstrates that our physical parameters of interest are minimally affected by 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 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 distribution of SFR and M d may be biased in the region of low values.
We stress, however, that our tests do not account for uncertainties stemming from prior model assumptions. For example, at fixed stellar initial mass function, the stellar population model and star formation history can introduce 0.3 dex uncertainty to the stellar mass and SFR (Conroy 2013). It is difficult to push the sensitivity of sSFR derived from integrated SEDs down to values 10 −12 yr −1 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 sSFR close to or below 10 −12 yr −1 . One should regard their SFRs with some caution.

COMPARISON WITH PREVIOUS WORK
This section compares the stellar masses and SFRs derived from our SED fitting with those from two widely used catalogs, the MPA-JHU database (Kauffmann et al. 2003;Brinchmann et al. 2004;Salim et al. 2007) and GSWLC-2 (Salim et al. 2018). The MPA-JHU catalog derives M * from SED fitting of optical photometry, using model magnitudes (ModelMag) of the five bands of SDSS (Salim et al. 2007). The total SFR of the galaxy comes from a combination of the fiber-based SFR of the central region (∼ 1 − 6 kpc for the 3 -diameter fiber in the redshift range z < 0.11 of our sample) of the galaxy ) and the SED-based SFR for the galaxy outskirts (Salim et al. 2007). SFR inside the fiber is estimated from the extinction-corrected Hα emission if the galaxy is purely star-forming according to the Baldwin et al. (1981) [O iii] λ5007/Hβ versus [N ii] λ6584/Hα diagram. If, on the other hand, the target shows nuclear activity or has weak Hα emission, the SFR is estimated from the strength of the 4000Å break inside the fiber, and from SED fitting outside. Salim et al. (2018) collect broadband (FUV to MIR) photometric data from preexisting catalog measurements provided by GALEX, SDSS, and WISE to derive M * and SFR from SED fitting. With the aid of a spectral template from Chary & Elbaz (2001), they convert the MIR flux densities from W3 or W4 to total IR (8 − 1000 µm) luminosity, which is used in the SED fitting with a customized version of CIGALE. A correction for AGN contribution to the IR luminosity is applied based on the strength of the [O iii] λ5007 line (Salim et al. 2016).
We derive the dust mass of low-redshift galaxies in Stripe 82 for all the sources using panchromatic (UV to FIR) SED fitting, primarily relying on data from the HerS survey (Viero et al. 2014). Our tests in Section 4.4 suggest that the upper limits in undetected bands usually can give useful constraints on the dust mass. Bertemes et al. (2018) also estimated dust masses for 78 massive galaxies in the Stripe 82 region with CO(1-0) measurements, all of which had 3σ detections in WISE W4 and Herschel 250 and 350 µm. Our dust mass measurements agree with those of Bertemes et al. (2018) to 0.05 ± 0.19 dex.

Comparison of Stellar Masses
The stellar masses from our work agree well with those of GSWLC-2 with negligible zero point offset and scatter (median and standard deviation ∆ log M * = 0.01 ± 0.11 dex; Figure 11a). The 0.1 dex scatter is consistent with the measurement dispersion of M * induced by nondetecions in 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 (∆ log M * = 0.11 ± 0.11 dex; Figure 11b), 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 (∆ log M * = 0.08 ± 0.12 dex; Figure 11c), 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 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.
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 (∆ log M * = 0.03 ± 0.19 dex), and MPA-JHU (∆ log M * = 0.11 ± 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.

Comparison of Star Formation Rates
For the sample as a whole, our SFRs show a relatively large scatter and modest zero point offset with respect to GSWLC-2 (∆ log SFR = −0.12 ± 0.56 dex; Figure 12a). Closer inspection reveals that most of the discrepancy is confined to the 45% of objects that have low star formation activity: while ∆ log SFR = −0.05 ± 0.23 dex at SFR > 0.3 M yr −1 , ∆ log SFR = −0.40 ± 0.71 dex at SFR 0.3 M yr −1 . Now, distinct from our work, which strives to provide full spectral coverage from 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 ∆ log SFR = 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 are ∼ 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 −1 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 SFR. 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 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 −1 ). 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. 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 12b, sources with SFR 0.5 M yr −1 track our measurements fairly well, with ∆ log SFR = −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. Al- though 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 (∆ log SFR = −0.47±0.71 dex) once SFR 0.2 M yr −1 . In the range SFR ≈ 0.001 − 0.1 M yr −1 , the SFRs in the MPA-JHU catalog remain more or less fixed at SFR ≈ 0.1 M yr −1 . 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 specific SFR (sSFR) calibrated from emission-line galaxies for which both quantities could be measured. The objects in our sample with low SFR 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 un-certainty 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 −1 .
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.

DISCUSSION
We study the panchromatic (UV to FIR) SEDs of massive galaxies in the SDSS Stripe 82 region using a set of newly measured aperture photometry. For galaxies that are spatially resolved, as is the case for the majority of our low-redshift sample, it is important to adopt a common, matched aperture for all bands to measure the total flux of the source. Careful consideration is given to deblending physical or projected neighboring galaxies, as well as to removing contaminating foreground stars. We pay close attention to the calculation of rigorous uncertainties, both statistical and systematic, and to the derivation of upper limits for nondetections. It is critical to include upper limits into the SED analysis in order to obtain unbiased measurements of physical parameters, especially the SFR and dust mass.
We derive stellar masses that are largely consistent with those reported in the GSWLC-2 and MPA-JHU catalogs. We observe a mild deviation between the MPA-JHU stellar masses and ours among low-mass galaxies, which can be attributed to differences in the choice of star formation history. The SFRs measured by our SED fitting are consistent with those provided by the GSWLC-2 and MPA-JHU catalogs for galaxies with relatively strong star formation activity (SFR 0.3 M yr −1 ). Both catalogs deviate from our measurements for galaxies with SFR 0.3 M yr −1 . It is not easy to fully explain the deviations between our SFRs and those from GSWLC-2 because the latter combines a set of heterogeneous photometric measurements, which at wavelengths longer than ∼ 20 µm are based on extrapolation from an average spectral template. By contrast, our SEDs are constructed from optimized measurements of actual observations up to 500 µm. The systematic deviations between our results and those from the MPA-JHU catalog can be traced to several factors. Apart from the increased uncertainty incurred from the hybrid method of deriving total SFRs by combining central fiber measurements with (incomplete) photometric measurements in the galaxy outskirts, the reliance on the relation between D4000 and sSFR to estimate SFRs for galaxies with weak emission lines and AGNs introduces not only significant scatter but also bias.
As a first exploratory application of our database, we examine the distribution of the Stripe 82 galaxy sample in the SFR − M * plane ( Figure 13). We construct the main sequence using the star-forming galaxies that have sSFR ≡ SFR/M * > 10 −11 yr −1 (Wetzel et al. 2012;Kukstas et al. 2020), calculating the median SFR in stellar mass bins of 0.15 dex, following Saintonge et al. (2016). We fit the SFR−M * relation with a second-order polynomial, which suffices to describe the data because of the limited number of galaxies at M * 10 11.3 M . The best-fit relation is where x = log (M * /10 10 M ). Our main sequence is overall close to, but slightly above that reported by Saintonge et al. (2016): at M * = 10 10.5 M , ∆ log SFR ≈ 0.1 dex ( Figure 13). Saintonge et al. (2016) derive SFRs using a combination of GALEX and AllWISE (W3 or W4) photometry. The W3 and W4 photometry, obtained using a Kron (1980) aperture, can miss 10%−50%  Saintonge et al. (2016) and Renzini & Peng (2015). Our main sequence is derived following the method in Saintonge et al. (2016).
of the total flux of highly concentrated galaxies (Graham & Driver 2005). This may explain their lower SFRs relative to ours because galaxies with M * 10 10 M usually have high concentration (Wuyts et al. 2011;Popesso et al. 2019).
We confirm the flattening of the slope of the main sequence relation, which is particularly notable for galaxies more massive than ∼ 10 10.5 M . Our main sequence lies substantially below the relation reported by Renzini & Peng (2015), by ∆ log SFR ≈ −0.3 dex at M * = 10 11 M , and is in better agreement with that of Saintonge et al. (2016). We suspect that part of the discrepancy, apart from differences resulting from the SFR estimator used, stems from the sample definition. The study of Renzini & Peng (2015) excluded galaxies classified as AGNs and composites. Since nearby AGNs preferentially inhabit earlier type, massive galaxies with relatively evolved stellar populations (Ho et al. 1997(Ho et al. , 2003, these sources mainly occupy the high-mass region with SFRs lower than those of star-forming galaxies (Leslie et al. 2016). Additional factors may contribute to the flattening of the main sequence at the high-mass end, including the growth of bulges (Wuyts et al. 2011;Abramson et al. 2014;Feldmann 2017;Belfiore et al. 2018), underestimation of SFR in edge-on disks (Morselli et al. 2016), and starvation of cold gas in massive halos (Popesso et al. 2019). We will explore these issues more fully in a forthcoming work.
To illustrate the importance of estimating proper SFRs for investigating the statistical physical properties of galaxies, in Figure 14 we show the distribution of Stripe 82 galaxies on the SFR − M * plane using the parameters derived from this study, compared to those based on the GSWLC-2 and MPA-JHU catalogs. While the distribution of the star-forming galaxies around the main sequence is largely consistent among the three sets of measurements, in detail there are subtle differences. The ridge line of star-forming galaxies at the massive end of GSWLC-2 is slightly lower than our main sequence. This is a consequence of the fact that GSWLC-2 systematically overestimates SFR over the range ∼ 0.1 − 1 M yr −1 (Figure 12a), whereby galaxies that are located in the so-called green valley according to our measurements become artificially elevated to the star-forming sequence in GSWLC-2. The contours of star-forming and green valley galaxies from our measurements are consistent with the contours derived from the MPA-JHU catalog (Figure 14b). In contrast, the distribution of quenched galaxies (e.g., SFR 0.1 M yr −1 ) show dramatic differences between the GSWLC-2 and MPA-JHU catalogs, and between these two catalogs and ours. The discrepancies stem from the different methods and systematic biases incurred when measuring the SFR, as discussed in Section 5.2. The grey dotted diagonal line marks log (sSFR/yr −1 ) = −12, below which the SFRs may not be trustworthy (see Section 4.5).

SUMMARY
We construct panchromatic SEDs of 2781 low-redshift (z = 0.01 − 0.11) galaxies in the SDSS Stripe 82 region, among them 2668 that have stellar masses M * > 10 10 M , that have been observed with the Herschel/SPIRE instrument at 250, 300, and 500 µm.
In combination with other survey data from GALEX, SDSS, 2MASS, and WISE, we develop a hybrid approach of matched-aperture photometry, which can be used effectively for the 14 shorter wavelength bands (FUV to 4.6 µm), along with profile-fitting photometry for the longer wavelength bands. This approach is optimized for relatively bright, nearby galaxies that are partly or fully resolved. We apply simultaneous, multiband, 2D image fitting to properly decompose interacting/merging galaxies and sources in crowded environments. A new method, based on random forest regression, is introduced to decontaminate foreground stars from images. Apart from obtaining reliable and physically consistent total galaxy fluxes and their associated uncertainties across multiple bands, we estimate robust upper limits for the nondetections and fully incorporate them into the final SED analysis. Mock tests evaluate the influence of certain critical bands and their respective upper limits. We derive stellar masses, SFRs, dust masses, and AGN luminosity fractions, which will be the subject of a series of forthcoming works.
The main results are as follows: 1. Our method of galaxy photometry yields selfconsistent, panchromatic total fluxes with systematic uncertainty less than ∼ 20%.
2. Although 40% of our sample with strong optical emission lines have some level of nonstellar nuclear activity according to their spectral classification, in the vast majority of the cases the AGN contributes negligibly to the broadband SED fits and can be neglected.
3. The stellar mass can be well constrained solely using photometry from the five SDSS bands. By contrast, estimating accurate SFRs requires the detection of or, at the very minimal, upper limits in the MIR or even bands of longer wavelengths. Measuring the dust mass requires detection in ar least one Herschel/SPIRE band, and whenever possible, even upper limits in the FIR bands should be used to avoid excessively large systematic biases.
4. Comparison of the stellar masses from our SED fitting reveals good agreement with those from the GSWLC-2 and MPA-JHU catalogs. The same is true for the SFRs of galaxies that occupy the star-forming main sequence, which, for our relatively massive galaxies correspond to SFR 0.3 M yr −1 . However, for galaxies with weaker star formation activity, including those on and below the green valley, the SFRs from the GSWLC-2 and MPA-JHU catalogs disagree with each other, and both sets of measurements are systematically overestimated relative to ours. SED fitting-based SFRs of galaxies with moderate to low star formation activity should be used with considerable caution.
LCH was supported by National Key R&D Program of China (2022YFF0503401), the National Science Foundation of China (11721303, 11991052, 12011540375, 12233001) and the China Manned Space Project (CMS-CSST-2021-A04, CMS-CSST-2021-A06). We acknowledge Yingjie Peng, Linhua Jiang, Chao Ma, and Kai Wang for their discussion and contributions to various aspects of this work. We thank the referee for a very constructive report that greatly helped to improve the paper. In this appendix, we compare our photometric measurements with those provided in the catalogs of GALEX, 2MASS, and WISE. For clarity, only sources that have flux S/N > 2 as given by the GALEX, 2MASS, and WISE catalogs are presented in the comparison. We use the fit left censoring 12 method in the KaplanMeierFitter class of Python to perform the survival analysis for the comparisons, calculating the median and 1σ of the flux differences after accounting for upper limits.
12 https://lifelines.readthedocs.io/en/latest/Survival%20analysis% 20with%20lifelines.html#left-censored-data-and-non-detection 2. SDSS: Figure A2 compares our measurements with the SDSS modelmag values. The SDSS values are fainter than our measurements by 0.08 mag in u and by 0.05 mag in g, although these systematic offsets are small compared to the scatter. In the r, i, and z bands, our measurements agree well with the modelmag values. The modelmag photometry for the five SDSS bands are based on the better of two model fits in the r band that assume either an exponential or de Vaucouleurs function. However, the bluer bands, which are dominated by younger stars, may have different light profiles. Our non-parametric measurements overcome this problem and recover the total flux in the SDSS bands.
3. 2MASS: We compare our measurements with two sets of 2MASS catalog results. Figure A3 shows good consistency (∆J = −0.04 ± 0.20 mag, ∆H = −0.04 ± 0.25 mag, ∆K s = −0.04 ± 0.24 mag) between our measurements and the "total flux" from the 2MASS XSC (Jarrett et al. 2000a). The XSC fitted the galaxy light profile with a Sérsic modified exponential function in the J band, which was then applied to all bands and integrated out to 4 times the disk scale length to measure the total flux of the galaxy. The good consistency confirms that our measurements enclose the total flux of the galaxy. However, only about 60% of our sample is included in the XSC. While PSF-fitting magnitudes are available for all sources in the 2MASS PSC (Cutri et al. 2003), the large, systematic discrepancy (∆J = −0.95 ± 0.45 mag, ∆H = −0.95 ± 0.48 mag, ∆K s = −0.88 ± 0.41 mag) with our measurements ( Figure A4) indicates that most of the galaxies are resolved significantly in 2MASS images.
4. AllWISE: For resolved galaxies flagged by ext flg > 0 in the AllWISE catalog, our measurements are about 0.5 mag brighter than the profilefitting results ( Figure A5). Even for galaxies classified as unresolved (ext flg = 0), our measurements are still 0.2 mag brighter in W1 and W2. For W3 and W4 bands, we only compare the measurements of resolved galaxies selected according to the criteria outlined in Appendix C. Our aperture photometry is about 0.4 mag brighter than the profile-fitting results. The elliptical aperture photometry of AllWISE, which is based on the 2MASS extended source aperture, still systematically misses some flux compared to our measurements. A similar result was found by Cluver et al. (2014).

B. REMOVING BLENDED EMISSION FROM FOREGROUND STARS
We develop a new method to remove the contamination of a foreground star if its emission cannot be deblended from the target galaxy following the standard deblending procedure described in Section 3.1. While the galaxy and foreground stars are well resolved and can usually be deblended in SDSS images, they become heavily blended in WISE images on account of its much coarser PSF. Blending can even be problematic in GALEX and 2MASS images, too. As robust measurements exist of the positions and magnitudes of the galaxy and the star(s) from SDSS, we can remove the contamination of the stellar emission if we can accurately estimate the magnitude of the star(s) in the other bands based on the SDSS measurements.
We use random forest (RF) regression (Breiman 2001) to predict the magnitude of the star in the UV and IR bands based on their five-band SDSS optical magnitudes. The RF algorithm is an effective machinelearning model for regression problems. It is a supervised learning algorithm that links the input individual subsamples of the data set and outputs predicted values by building multiple decision trees. The decision trees work through minimizing the Gini Impurity (Pedregosa et al. 2012) at each stage of the process. Models are constructed to uncover the highly nonlinear correlation between the targeted output value and the input subsamples that may contribute to it (Bluck et al. 2022). We will show that the optical SED from the five bands of SDSS contains enough information to strongly constrain the stellar type and accurately predict the magnitude of the star in the other bands.
We employ the Python function RandomForest Regressor 13 from the package scikit-learn (Pedregosa et al. 2012). To construct the training samples and test the feasibility of the method, we collect unsaturated (z > 13 and u > 16 mag) field stars from the SDSS DR16 database (class=STAR) that do not have companion sources within 20 , a distance larger than 3 times the FWHM of the PSF. We also exclude stars near the Galactic plane (|b| < 10 • ) because dust extinction strongly affects the optical SED of the star. We crossmatch the SDSS coordinates of the stars to the GALEX GR6+7, the 2MASS PSC, and the AllWISE catalog. Limiting to measurements with S/N > 3, we find a sample of more than 21,000 stars for the UV and 54,000 stars for the IR, which covers the bands J through W3 but not W4 and the FIR, for which stellar emission is typically very faint.
We arbitrarily divide the stars into a training sample and a testing sample, roughly in the number ratio 8:2. For each of the eight bands (FUV, NUV, J, H, K s , W1, W2, W3), we use the training sample to train an RF model to make predictions based on the SDSS bands. Then we use the testing sample to evaluate the effectiveness of the model. We adopt n estimators=300 instead of the default value of 100 when training the FUV and NUV data, because we find that this choice significantly improves the prediction for the UV bands but not for the IR bands. Otherwise, we adopt the default parameters of RandomForestRegressor. Figure B1 compares the  measured and predicted magnitudes of stars in the testing sample. For each band, we calculate the R 2 score, where k is the number of data in the test sample, and y i , y i andŷ are the truth value, prediction, and mean of the truth value. The best possible R 2 score for a model is 1. Our models show very high consistency between data and prediction. We find little systematic uncertainty between the measured and predicted magnitudes for stars with z < 18.5 and u < 18.5 mag. As the contamination is negligible if the star is fainter than these limits, we only consider the contamination from stars brighter than these criteria. The scatter between the measured and predicted magnitudes is typically 0.1 mag, except for FUV (0.4 mag) and NUV (0.3 mag). We take the scatter as the uncertainty of the prediction. For each target, we search for the stars with z < 18.5 and u < 18.5 mag inside twice the aperture size from the center of the galaxy. If the star cannot be deblended from the target in the UV and IR images following the procedure of Section 3.1 (e.g., Figure B2), we use the SDSS magnitudes and RF model to predict the magnitude of the star in these images. We calculate the fraction of a star's flux (f star ) within the aperture of the target in the unmasked region from the PSF models, which are provided by the GALEX and WISE web-  sites 14 . The PSF of 2MASS varies from image to image, and for each image tile we extract stars with S/N > 3 in an uncrowded field to generate the PSF model using the IRAFstarFINDER function from photutils 15 . 14 http://www.galex.caltech.edu/researcher/techdoc-ch5.html and https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/ sec4 4c.html 15 https://photutils.readthedocs.io/en/stable/epsf.html The uncertainty introduced by star contamination is where σ RF is the scatter of the RF prediction in the corresponding band. We exclude targets whose SDSS images contain a saturated star (z < 13 or u < 16 mag) within twice their aperture size.

C. RESOLVED SOURCES IN W3 AND W4
In the AllWISE catalog, sources are flagged as resolved (ext flg > 0) if they are included in the 2MASS   extended source catalog or PSF-fitting provides poor results (χ 2 ν > 2). This subjective flag cannot be used to guide our selection of galaxies that are resolved in W3 and W4. To obtain a more reliable empirical guide, we apply the PSF-fitting method to mock galaxies to determine the critical optical size above which PSF fitting fails to yield a robust measure of galaxy flux. Guided by the measurements of Bottrell et al. (2019), we construct mock galaxy images using a single-component model with Sérsic index n = 1 − 5 and effective radius R e = 1 − 10 . To mimic the actual observations, we set W3 = 8.5 − 11.5 mag and W4 = 5.5 − 7.5 mag, and we convolve the mock images with the corresponding PSFs. Employing the Levenberg-Marquardt algorithm for the PSF fitting, we follow the same approach as All-WISE 16 and use as an initial guess a circular aperture of radius 8. 25 for W3 and 16. 5 for W4. The fit considers confusion uncertainty and pixel value uncertainty. Figure C1 illustrates that PSF fitting can measure > 80% of the galaxy flux in both W3 and W4 if the galaxy has high Sérsic index (n ≥ 3.5). However, if the galaxy is less centrally concentrated, the flux measured by PSF fitting drops quickly below 80% if the galaxy has R e > 2 in W3 or R e > 4 in W4. Therefore, we consider the galaxy possibly resolved in W3 and W4 if it has R e larger than 2 and 4 , respectively, and n < 3.5 in the SDSS r band. For galaxies without R e measured by Bottrell et al. (2019), we establish an empirical correlation between R e and the Petrosian radius R 50 in the r band ( Figure C2), for which a linear regression LinMix (Kelly 2007) yields R 50 = 0.49 R e + 0.89.
Thus, a galaxy is resolved in W3 if R 50 > 1. 9 and in W4 if R 50 > 2. 9. To facilitate other investigators who might wish to correct the PSF-fitting fluxes of galaxies from the AllWISE catalog, we fit linear functions to the flux recovery fractions in W3 and W4 as a function of R e , separately for galaxies with n < 3.5 and n ≥ 3.5. The best-fit functions are shown in Figure C1.  Figure C2. Relation between Sérsic Re and Petrosian R50 for our galaxies. The black dashed line is the 1-to-1 relation. The blue line indicates the best-fit linear relation.