Identifying and Characterizing the Most Heavily Dust-Obscured Galaxies at $1\leq z\leq 4$

There is a well known degeneracy in reddening due to redshift, stellar age, and dust attenuation. Any successful photometric redshift estimation must therefore allow for the full range of these effects. Marchesini et al. (2010) were the first to show using the photometric redshift code EAZY (Brammer et al., 2008) that introducing a template that is both old and dusty moves a significant number of sources that would be incorrectly identified as massive galaxies at $z>3$ to $z\sim 2-3$. The effects introduced by including such a template were further explored in Muzzin et al. (2013b) and Marchesini et al. (2014). With sufficient dust, the ability to detect the Lyman break is reduced, and the Balmer and 4000 Å breaks become broadened, resulting in both less accurate and less precise photometric redshifts. These findings illustrate the importance of including the widest possible range in stellar population parameters/templates when calculating photometric redshifts, or equivalently, that the adopted template set spans all the possible ranges in the color-color spaces as real galaxies.

We choose to re-derive the redshifts for our UltraVISTA sample rather than use the default EAZY redshifts from the catalog in order to ensure that we include photometric solutions for the most extreme cases of dust obscuration. For an initial sample to perform our revised redshift modeling, we select sources from the UltraVISTA DR3 v3.3 catalogs with $K_{S}<25$ and colors $H-K>0$ and $H-[3.6]>0.5$. These color cuts are intentionally relaxed compared to those used to select extremely red objects (EROs) and similar objects in the literature in order to be as comprehensive as possible in our selection of dusty objects. These cuts will only exclude blue and obviously unobscured sources from our sample. Additionally, a signal to noise ratio cut at 5 is applied to the $H-$ and $[3.6]-$ bands in order to ensure a reliable color selection (the magnitude cut ensures a signal to noise ratio greater than five for the $K_{S}$-band). Figure 1 shows the distribution of $H-K$ and $H-[3.6]$ colors as a function of redshift for the full UltraVISTA parent sample ($\sim 60,000$ sources) with the color cuts indicated as red horizontal dashed lines. Redshifts in this figure are the standard v3.3 version calculated with EAZY. Additionally, each panel shows the color evolution with redshift for the EAZY template set used to derive the redshifts. The twelve templates shown are meant to span a wide range of galaxy types such that the SED of an observed source can be fit through linear combination of the templates. The color tracks show that EAZY requires the reddest templates in order to match the observed colors for the most extreme sources, and that even for these objects, the template set is pushed to its limit. There is an inherent trade-off in constructing a template set that can act as a basis for a diverse sample of galaxies, simultaneously covering the full range of observed SEDs and limiting the number of templates to prevent over-fitting. It is expected that edge cases may not be fit well by the templates. For this reason, we choose to re-derive redshifts for the sample meeting the above color cuts allowing for a wider range in stellar population model parameters.

The 30,205 sources which meet both color selections are modeled with FAST++ (Kriek et al., 2009; Schreiber et al., 2018) with a Calzetti dust attenuation law (Calzetti et al., 2000), a Chabrier (2003) IMF, delayed exponentially declining star formation history, and allowing $0\leq$A${}_{V}\leq 6$ in order to redetermine their redshifts. FAST++ utilizes the full list of 49 photometric bands spanning the UV to Spitzer IRAC 8$\micron$ in the UltraVISTA DR3 catalog. In addition, FAST++ also allows the user to provide an IR luminosity to help break the age-dust degeneracy inherent when fitting UV-NIR data. Providing a luminosity requires knowing the redshift beforehand. Our methodology is to fit the UV-8$\micron$ photometry to determine the redshift, then perform a second fit with the redshift fixed while providing an IR luminosity obtained by scaling the observed MIPS 24$\micron$ flux using the calibration of Dale & Helou (2002). This obviously does not provide an optimum estimate of LIR, but does provide a way to help break the age-dust degeneracy.

Through comparison with the original EAZY redshifts we decided upon the following approach. Given the default FAST++ setup with these modeling assumptions, we find that $\sim 2\%$ of sources have failed fits. These failed fits either return models which fail to match the observed photometry or have unrealistic physical parameters. For these sources, we rerun FAST++ without a template error function to obtain a redshift. With this redshift, we run FAST++ with the redshift fixed at this value and re-instituting the template error function to derive the stellar population parameters. For the majority of initially failed fits, this process provides reasonable results. Figure 2 shows the results of our comparison of $z_{\rm FAST++}$ and $z_{\rm EAZY}$ with coloring according to the A_V derived by FAST++. The black line indicates the 1:1 relation, and the dashed lines are at factors of 0.5 and 1.5 below and above. In general there is good agreement. Although there are always systematic differences between photometric redshift codes, the agreement for the bulk of the sample shows our FAST++ redshifts to be reliable. We remove outliers lying outside the dashed lines as we find they are almost universally failed fits. Of particular interest to this work are sources for which FAST++ prefers a solution with high levels of obscuration at lower redshift, which can be prominently seen in the $z=1-3$ range. These sources would potentially be missed with modeling that does not allow sufficient levels of dust obscuration.

Figure 3 shows the relations between $z$, A_V, and M_∗ for the sources we model with FAST++. Contours and small gray points show the distributions for all sources we remodeled, while colored points indicate our final selection. Using the output parameters from FAST++, we select sources with $1\leq z\leq 4$, A${}_{V}\geq 3$, and log(M${}_{*}/\hbox{M${}_{\odot}$})\geq 10.5$. The attenuation cut can be seen in Figure 3 to select only the most highly obscured sources from the parent sample, with the distribution sharply dropping off after this point. We found no reliable candidates beyond $z=4$, and so adopt this value as our limiting redshift. Finally, we chose a mass limit that allows us to be mass-complete in our highest redshift bin. This additional selection is a requirement to properly interpret any trends with stellar mass, including the star-forming main sequence and size-mass relation (see Sections 4.3 and 4.4). This process results in a final sample of 65 sources. Figure 3 illustrates that these sources indeed lie at the extremes of physical parameter space. These sources are matched to the Herschel FIR data as described in Section 2. Table 1 lists the detection rates for the MIPS 24 $\micron$ and Herschel bands for both our final sample and sources from the UltraVISTA catalog with $K_{S}<25$ in the same redshift range. The extraordinarily high detection rates for our sample compared to the parent $K_{S}$-band selected sample indicate that our selection method is efficient in selecting galaxies with strong dust emission based only on the UV-to-8 $\micron$ photometry.

To test the dependence of our selection on the assumed attenuation curve, we redetermine the redshifts and physical properties of our final sample using the Kriek & Conroy (2013) attenuation curve. With the Kriek & Conroy (2013) curve, we find that 15 of the 65 sources have disagreeing redshifts using the criterion $\Delta z/(1+z)>0.1$. When comparing the physical parameters for sources with agreeing redshifts, we find a median offset of 0.2 mag lower A_V and 0.16 dex higher M_∗. Three of the 15 sources with disagreeing redshifts remain within the sample selection criteria. Combined with the MAGPHYS modeling (below) which uses both a different dust attenuation law and star formation history, we conclude that our SED modeling technique efficiently selects ultra-dusty galaxies, but that the precise physical properties derived will obviously depend on the assumed attenuation law. We conclude that since our selection is based on these derived properties, the exact sample selection is in part driven by model choice, but that most galaxies which meet our selection are genuinely very dusty.

At this stage we introduce a comparison sample of sources detected with Herschel . From the parent UltraVISTA catalog, we select sources with identical mass and redshift cuts to the ultra-dusty sample, but in place of the A_V selection, we require a 3$\sigma$ detection in at least one of the five Herschel bands (note this also requires a MIPS $24\micron$ detection). This results in a sample of 1616 sources that are selected based on their dust emission rather than dust attenuation.

As an additional check of the robustness of our selection, we examine the rest-frame UVJ colors of our final sample of sources in Figure 4. Points are colored by their A_V value. We find that all but four of our sources (94%) fall in the dusty star-forming region of the UVJ diagram. The remaining four reside in the quiescent region, although we point out that these sources lie near the boundary. We conclude that our selection based on the FAST++ A_V is generally consistent with a UVJ color selection when the attenuation reaches our chosen level of 3 mag. Comparatively, in Martis et al. (2019) we demonstrate that dusty galaxies with less severe attenuation A${}_{V}\gtrsim 1$mag) occupy the quiescent and non-dusty star-forming regions of the UVJ diagram in higher numbers. We also show the UVJ colors of the Herschel comparison sample in gray. The Herschel sample occupies almost exclusively the star-forming region of the diagram, but with significantly bluer colors than the ultra-dusty sample. This suggest our sample is extreme in color and dust attenuation not only with respect to the overall galaxy population, but compared to FIR-selected galaxies as well.

We remodel the final sample of 65 sources using the high-z extension of MAGPHYS (da Cunha et al., 2008, 2015) including the matched Herschel data in order to obtain intrinsic, unattenuated SEDs to which we compare the observed SEDs. The MAGPHYS modeling utilizes most of the UV-MIR bands as well as fully incorporating the MIPS 24$\micron$ and Herschel measurements for 49 bands total. Sources without matched Herschel counterparts are assigned upper limits for the Herschel bands. MAGPHYS treats the upper limits by setting the flux to zero and the error to the upper limit for any nondetections. We use the photometric redshifts from FAST++ when computing the models. For the values of the physical parameters, we adopt the medians of the output distributions. See Appendix A for further description of the MAGPHYS modeling.

Figure 5 shows the comparison of physical parameters derived by FAST++ and MAGPHYS. We observe that the stellar population parameters are comparable to the FAST++ output values, although we do note a systematic increase in the derived stellar masses of $\sim 0.2$ dex when using MAGPHYS, with the offset being more pronounced for our extreme dusty galaxies. From the above testing, we found that a differing attenuation law can cause discrepancies of this magnitude, but the degree to which the star formation history or inclusion of FIR data contribute is not clear. We note that the star formation histories used in our FAST++ modeling and MAGPHYS are qualitatively similar, allowing a rising segment at early times, and exponential decline at late times. The SFR from FAST++ and MAGPHYS are in generally good agreement, although it appears there is a tail of sources with SFR from MAGPHYS being larger, potentially caused by the inclusion of the IR information and better accounting of dust-obscured star formation. A_V from FAST++ and MAGPHYS are also surprisingly well correlated, albeit with a large scatter of 0.35 mag. In general, A_V and SFR are notoriously the most uncertainly derived properties from SED modeling (e.g., Muzzin et al., 2009).

The distributions of key physical parameters of the MAGPHYS fits are shown in Figure 6. Since MAGPHYS assumes a different treatment for the dust (Charlot & Fall, 2000), the fact that only eight of the 65 sources in the final sample fall below $A_{V}=3$ mag reinforces the robustness of our selection of heavily dust-obscured objects. All 65 objects are classified as (U)LIRGs according to the output L_dust from MAGPHYS, further supporting this conclusion. The SFRs for this sample span a wide range with a median of 95 M_⊙ yr^-1. We also show the physical parameter distributions derived from MAGPHYS for the Herschel comparison sample, scaled to match the number of objects in the ultra-dusty sample. Several interesting points immediately emerge. The Herschel sample tends to slightly lower stellar masses ($\sim 0.2$ dex), in broad agreement with observed correlations between stellar mass and dust attenuation (e.g. Whitaker et al., 2010). Interestingly, the comparison is well matched in star-formation rate and dust luminosity, meaning that the extreme optical attenuation levels do not significantly bias the sample to higher levels of dust emission than typical Herschel samples. Finally, we find that our ultra-dusty sample exhibits typical A_V close to two full magnitudes higher than the Herschel sample. We see that the high-$A_{V}$ tail of the Herschel sample reaches these attenuation levels, but only in small numbers. Such a radical difference suggests that these objects either differ from other samples of dusty star-forming galaxies in terms of their dust properties or represent the high obscuration tail of currently known populations of dusty galaxies. This may be caused by different dust composition, geometry, mass, or any combination of the above. Some possible explanations are explored below. For more discussion of how the Herschel comparison sample relates to our selection, please see Appendix D.

A selection at A${}_{V}\geq 3$ implies that fewer than $\sim 5\%$ of optical photons are able to escape their region of production. Such extreme attenuation levels are usually only observed within the central regions of starbursts or in some SMGs (e.g., Casey et al., 2017; Calabrò et al., 2018; Schreiber et al., 2018), but the presented analysis suggests a sample of objects with similar attenuation averaged over the entire galaxy. To illustrate the extreme effect of dust on these galaxies’ SEDs, we show in Figure 7 the attenuated (red) and unattenuated (blue) SEDs for our sample derived from the MAGPHYS modeling. The solid curves show the median SEDs for all sources in the sample shifted to the rest-frame. The shaded regions indicate the range from the $15^{th}$ to $85^{th}$ percentiles. The SEDs are normalized to the rest-frame $J$-band, effectively normalizing by stellar mass. Individual photometric observations are shown as small gray points. For the Herschel bands, 3$\sigma$ upper limits are indicated by triangles in the case of non-detections. The median photometry as well as its $15^{th}$ to $85^{th}$ percentiles are plotted as large black points with errorbars. The downward arrows for the Herschel points indicate that the 3$\sigma$ limits are used in the calculation of the median. By visualizing the SED in terms of luminosity, one can see that almost all of the energy output for these objects occurs at IR wavelengths. In agreement with the FAST++ A_V selection, the attenuated and unattenuated SEDs’ luminosity differ by a factor of $\sim 100$ in the optical region, and even more in the UV (the normalization in Figure 9 somewhat reduces the appearance of this). We also show for comparison the median best fit FAST++ SED and $15^{th}$ to $85{th}$ percentiles in orange. In the overlapping region, the models agree well, with FAST++ preferring somewhat larger luminosities in the FUV.

The physical properties of our sample do not appear to align with typical Herschel -selected or SMG galaxy samples. Typical SMGs have A${}_{V}\sim 2$ (Hainline et al., 2011; Simpson et al., 2014; da Cunha et al., 2015). Whereas there are examples of high levels of obscuration similar to those we observe in some of these samples (Ma et al., 2015, 2019), they also typically have SFRs on the order of hundreds of solar masses per year, which few of our sources reach. Casey et al. (2017) present NIR spectroscopy of a sample of 450 and 850 $\micron$-selected sources with a median stellar mass of $4.9\times 10^{10}$ M_⊙, SFR of 160 M_⊙ yr^-1 and A_V of 5, which much more closely align with the derived properties of our sample. Importantly, they find that 11/35 of their spectroscopically confirmed sources have incorrect or missing photometric redshifts, stressing the need for special care to be taken in determining photometric redshifts of such heavily obscured sources. They argue that the spatial disconnect of obscured and unobscured regions poses the largest problem in this regard. Addressing this issue for our sample would require high resolution imaging of the thermal dust continuum to identify the physical regions responsible for the observed IR emission.

Simpson et al. (2014) and Elbaz et al. (2018) also argue strongly that the spatial configuration of stars and dust presents a risk to interpreting panchromatic SED modeling. In their analysis of a sample of ALMA-selected SMGs, the former find that the stellar light extends to radii 3-4x further out than the dust emission of the central starburst, meaning that energy balance arguments will no longer apply, and making estimates of physical properties unreliable. In fact, when they estimate dust extinction using a hydrogen column density inferred from their gas mass estimates, they find an average A${}_{V}=540$ mag. Clearly this disagrees strongly with any results obtained through SED modeling. Resolution of this issue at high redshift will require detailed, spatially resolved studies of the stars and dust, which will likely require observations with the James Webb Space Telescope (JWST ).

We wish to note the presence of active galactic nuclei (AGN) for two reasons. First, an AGN can contribute significantly to the SED of a galaxy, potentially leading to biased estimates of stellar and dust properties. Second, we are intrinsically interested in the prevalence of AGN within this population. We identify AGN within our final sample using the same criteria as Martis et al. (2019). Briefly, objects are flagged as AGN based on an IRAC color selection, radio excess, and X-ray emission. In all, 12 out of 65 (18%) meet at least one of the selection criteria. Of these 12, eight meet the IRAC color selection, seven meet the radio selection, and two are selected via their X-ray emission. Given their extreme obscuration levels, it is not surprising that the minority of the identified AGN (2/12) are X-ray AGN. Two of the sources identified as AGN have SFRs in excess of 1000 $\hbox{M${}_{\odot}$}yr^{-1}$ according to both the FAST and MAGPHYS modeling, indicating that at least some sources may have their SFRs overestimated due to an AGN contribution. In all figures, AGN are identified as different symbols/colors.

Previous work has shown that both dust mass and attenuation correlate with SFR (e.g. da Cunha et al., 2010; Nagaraj et al., 2021, 2022, but see van der Giessen et al. (2022)) Given the significant amount of dust attenuation in this sample, it would be reasonable to expect high SFRs. Additionally given we expect some correlation between dust attenuation and emission, highly attenuated sources would also have significant levels of dust emission, which can arise from star formation, AGN, or older stellar populations. If star formation is the primary driver of high levels of dust emission, then the source may lie above the star-forming main sequence.

Figure 8 shows the relation between SFR and M_∗ for our final sample as well as the comparison Herschel sample along with several recent determinations of the star-forming main sequence from the literature. The SFR and M_∗ shown are those derived from MAGPHYS. We find that the majority of our sources are consistent with lying on the main sequence and that they are not separated in parameter space from the Herschel comparison sample. This may be a somewhat surprising result, since starbursts in the relevant redshift range are frequently observed to be heavily obscured and individually detected with Herschel or in the sub-millimeter (e.g., Schreiber et al., 2015; Miettinen et al., 2017; Elbaz et al., 2018). Martis et al. (2019), however, shows that a selection in dust obscuration leads to a sample with diverse levels of star formation. This work suggests that this holds true for stricter selections in dust obscuration as well. These results are also consistent with those of Penner et al. (2012) who show that their sample of DOGs which have $10^{12}\lesssim L_{\rm IR}\lesssim 10^{13}$ exhibit a range in star formation activity. Only 24% of their sample have sSFRs which indicate the dominance of a compact starburst.

Despite the modest SFRs for our sample, we still observe significant levels of IR emission. The two component dust model utilized by MAGPHYS allows slightly older stellar populations to contribute to IR emission by heating dust in the diffuse interstellar medium rather than dust in birth clouds around star-forming regions. Recent work suggests that intermediate age stars can contribute significantly to MIR emission (Utomo et al., 2014; Leja et al., 2019; Roebuck et al., 2019; Martis et al., 2019). The level of dust heating by intermediate age stars will depend on the lifetime of dust grains generated by different sources. The observed correlation between SFR and dust mass (da Cunha et al., 2010) suggests that at least some dust is produced quickly by supernovae during an episode of star formation, but the degree to which older asymptotic giant branch stars can continue to produce dust well after star formation has ended remains unclear (Goldman et al., 2021). Depending on the dust lifetime, it may be possible that some of our sources are post-starburst galaxies which have begun to decline in SFR before the large amount of dust associated with the burst has been destroyed. Given the high levels of obscuration, and hence the relatively featureless SEDs, measuring precise star formation histories for these objects to test this hypothesis will prove extremely challenging.

In order to better understand the physical arrangement giving rise to these conditions, we investigate the HST $F814W$ and $F160W$ images of all sources available from the COSMOS (Koekemoer et al., 2007; Massey et al., 2010) and DASH (Drift And SHift Momcheva et al., 2017; Mowla et al., 2019) programs’ coverage of the UltraVISTA survey area. Figure 9 shows three example sources. From left to right we show the UltraVISTA $K_{S}$, $F814W$ and $F160W$ images, as well as the UltraVISTA catalog number and MAGPHYS physical properties. We observe a range of morphologies, with some sources appearing round and compact, whereas others are more elongated suggesting a possible disk-like morphology observed edge-on. There are indications of possible mergers/interactions in progress for source IDs 89509, 159508, and 198121.

To quantify these observations, we perform profile fitting of the $K_{S}$ and $F160W$ images with the code GALFIT (Peng et al., 2010). Each source is modeled with a single Sersic profile except for the sources identified above as clearly decomposing into multiple sources in the the higher resolution HST imaging. In these cases we model each component with a Sersic model. Given that HST imaging is only available for 23 out of 65 of our sources, we reference the fitting results from the $K_{S}$-band images. Appendix C shows that despite the lower resolution, differences in the fit parameters for the quantities we discuss here do not affect our conclusions.

The top panel of Figure 10 shows the distribution of axis ratios for our sample measured from the GALFIT modeling of the $K_{S}$-band images. We find examples of sources across the full range of galaxy shapes. The peak of the distribution occurs around an axis ratio of 0.4. Elongated sources with these lower axis ratios are likely to be disks observed edge-on. Previous work shows that viewing a star-forming disk edge-on results in a higher observed attenuation level (Wang et al., 2018), so this provides a natural explanation for the extreme attenuation level of our sample.

This leaves open at least two possibilities for our sources with high axis ratios. They may represent the same population of disks, but viewed face-on, or they may be more compact starbursts such as the "blue nuggets" observed in this redshift range (Osborne et al., 2020; Zolotov et al., 2015). To begin to address this question we show in the bottom panel of Figure 10 the effective radius versus axis ratio for our sample. If these "blue nuggets" do form the rounder portion of our sample, we may expect an anti-correlation between size and axis ratio, such that rounder objects are smaller. If we discount the two sources with the highest axis ratios and large errors, then one may tentatively observe this anti-correlation. We further address this question in the next section.

In a study of starburst galaxies at $0.5<z<0.9$, Calabrò et al. (2018) find a wide range of obscuration values depending on location within the galaxy, with the starburst core reaching A${}_{V}=2-30$ mag. Given that the majority of their sample are morphologically classified as mergers, they argue that these extreme obscuration values may serve as a useful tool for identifying mergers at higher redshift. Our available HST images do not allow us to rule out merger classifications for the majority of our sources, but since only three show signs of neighboring sources, mergers do not appear necessary to generate the levels of obscuration we observe.