The statistical properties of 28 IR-bright dust-obscured galaxies and SED modelling using CIGALE

Abstract

The aim of this study is to characterize the physical and statistical properties of a sample of infrared-bright dust-obscured galaxies (DOGs) by fitting their spectral energy distributions (SEDs). We examined 28 DOGs at redshifts 0.47 ≤ z ≤ 1.63 discovered by combining images of the Subaru Hyper Suprime-Cam (HSC) survey, the VISTA Kilo-degree Infrared Galaxy (VIKING) survey, and the Wide-field Infrared Survey Explorer (WISE) all-sky survey, and detected at Herschel Spectral and Photometric Imaging Receiver (SPIRE) bands. We have detected a significant active galactic nucleus (AGN) contribution to the mid-infrared luminosity for |$71\%$| of DOGs. Our DOGs contain several types of AGNs; the majority of AGN emission originates from Type 2 AGNs. Our DOG sample shows very high values of stellar mass [log (M_*/M|${}_\odot $|⁠) = 11.49 ± 1.61] compared with other samples of DOGs selected at infrared wavelengths. Our study is promising to identify a new type of DOGs called “overweight” DOGs (ODOGs). ODOGs may indicate the end of the DOG phase, and then they become visible quasars. Principal component (PC) analysis is applied to reduce the number of dimensions of our sample, removing the dependency on the observed variables. There are two significant PCs describing |$72.7\%$| of the total variance. The first PC strongly correlates with redshift, dust luminosity, dust mass, and stellar mass, while far-ultraviolet (FUV) attenuation strongly correlates with the second PC, which is orthogonal to the first one. The partial correlation between the resulted physical parameters is tested, supporting the reliability of the correlations.

1 Introduction

Following the initial mid-infrared studies of extragalactic sources (Low & Kleinmann 1968), the all-sky survey of the Infrared Astronomical Satellite (IRAS) achieved a great accomplishment by discovering many galaxies producing much more energy in the infrared region than in the optical range (Neugebauer et al. 1984). This was the first time that these kinds of galaxies had been seen. After analyzing the data, researchers found that the impact of active galactic nuclei (AGNs; Lutz et al. 1998; Genzel et al. 1998) on the infrared (IR) light at low redshifts appears to rise with this IR luminosity. Furthermore, at higher redshifts, the high luminosity of ultra-luminous infrared galaxies (ULIRGs; L_IR > 10¹²L|${}_\odot $|⁠; Sanders & Mirabel 1996) is not yet completely understood and shows substantial variability in the AGN-to-starburst ratio (e.g., Sanders 1999; Desai et al. 2007; Menéndez-Delmestre et al. 2009; Pozzi et al. 2012). This turns out to be a significant challenge for us in comprehending these sources.

In the research conducted on the formation and evolution of galaxies, it has been hypothesized that enormous galaxies may have grown alongside the supermassive black holes (SMBHs) that reside at their galactic centers (Sanders et al. 1990; Hopkins et al. 2006, 2008). The gas-rich mergers or violent disc instabilities are what set off dense starbursts (Barnes & Hernquist 1992; Hopkins et al. 2008; Dekel et al. 2009). These events also supply the fuel for the accretion of the central black hole. The host galaxy and the SMBH develop in tandem, passing through phases characterized by starburst dominance, AGNs/QSOs, starburst composites, and AGN dominance until the AGN feedback becomes powerful enough to cause the ejection of gas and dust. This results in the cessation of star formation and the activity of AGNs in a very short timeframe and ultimately leads to the emergence of a passively evolved galaxy (Sanders & Mirabel 1996; Granato et al. 2004; Hopkins et al. 2006, 2008; Alexander & Hickox 2012). During the extreme star formation process, a tremendous quantity of dust is produced, which plays a significant part in sculpting the observable spectral energy distribution (SED) of a massive galaxy going through distinct stages. The vast majority of ultraviolet (UV) and optical photons are absorbed by dust, and these photons are then re-emitted at wavelengths that are far-infrared (FIR) and sub-millimeter (sub-mm). As a consequence of this, starburst-dominated and AGN–starburst composite systems would seem to be IR-luminous, precisely like the populations that have already been detected, which are: ULIRGs, sub-mm galaxies (SMGs; Blain et al. 2002) and dust-obscured galaxies (DOGs; Dey et al. 2008). The analysis of galaxies that are bright in the IR and have a high redshift can help in understanding the difficult situations that occur during the initial stages of the formation of massive galaxies.

The color-based parameters were suggested by Dey et al. (2008) for the optimal evaluation of a statistically meaningful sampling of dusty ULIRGs at z ∼ 1.5–3. They chose a set of visually dim (22 < R_Vega < 27) DOGs and mid-IR bright [F_ν(24 μm) > 0.3 mJy] DOGs, and determined those sources to have F_ν(24 μm)/F_ν(R) > 1000. Star-formation rate (SFR) density and SMBH growth rate peak simultaneously in the number density of DOGs at z ∼ 1–2 (Richards et al. 2006; Dey et al. 2008; Madau & Dickinson 2014). This lends credence to the hypothesis that the most active objects regarding the co-evolution of SMBHs and galaxies are connected to DOGs. In this perspective, DOGs that have a strong IR luminosity might potentially provide refuge for an SMBH that is quickly growing, and as a result, these DOGs are important in order to understand the evolution of galaxies and SMBHs.

In accordance with the SEDs of DOGs, DOGs are divided into two different subclasses: “bump DOGs” and “power-law (PL) DOGs” (Dey et al. 2008). The PL DOGs display a power-law characteristic in their SEDs, while the bump DOGs reveal a rest-frame 1.6 μm spectacular bump in their SEDs. The bump DOGs belong to galaxies that are operating in the star-forming (SF) phase (Desai et al. 2009; Bussmann et al. 2011), whereas the PL DOGs are believed to relate to galaxies that are operating in the AGN phase (Fiore et al. 2008; Bussmann et al. 2009b; Melbourne et al. 2012). When there is a greater amount of mid-infrared (MIR) flux density, the percentage of PL DOGs through all DOGs will rise (e.g., Dey et al. 2008; Toba et al. 2015).

The transformation of space points into a place with fewer dimensions is the fundamental purpose of the approaches for reducing the dimensions of space. This must be done without the loss of any vital information. Unsupervised and supervised methods are the two primary methods that may be used, respectively. When using the unsupervised method, it is not necessary to provide labels for the various data classes. While using the supervised method, the dimensionality reduction algorithms consider the class labels. The unsupervised method’s strategy that has proven to be the most successful so far is called the Principal Component Analysis (PCA; Turk & Pentland 1991; Belhumeur et al. 1997).

PCA has been applied rather often in the field of astronomy to investigate ISM turbulence, galaxy nuclei (Steiner et al. 2009), and protostellar jets (Cerqueira et al. 2015), as well having other applications. PCA is performed on a collection of data that contains several variables. It does this by increasing the variance, which creates additional orthogonal variables called main components (PCs). One way to represent the original variables is to use a linear combination of principal components. PCA is a method that is used to distinguish unnecessary information from information that is important within the original variables.

This paper presents a sample of 28 DOGs published in the Noboriguchi et al. (2019) catalogue. These DOGs are selected by combining the WISE all-sky data and deep optical imaging data obtained with the Subaru Hyper Suprime-Cam (HSC). We have used the Herschel Spectral and Photometric Imaging Receiver (SPIRE) point source catalog (FIR). Optical and NIR counterparts exist for all utilized DOGs. By using the optical to FIR spectra multi-wavelength datum, we can determine the AGN proportion which contributes to their IR brightness, and global physical properties of DOGs by fitting their SEDs, including stellar mass, SFR, dust mass, and luminosity, and others. The statistical properties are studied for our sample via the PCA method by using packages from the R project.¹

The format of this paper is as follows. Section 2 contains a description of the data and the sample selection. The SED fitting approach and the models used in the Code Investigating GALaxy Emission (CIGALE) can be found in section 3. Section 4 includes the results about the main physical properties of the DOGs, concerning the fractional AGN contribution to the infrared emission and the AGN types with regard to DOGs' activity. A comparison of the DOGs to other galaxy samples is presented in section 5. The statistical properties of our sample are shown in section 6. Finally, the summary of this work and the conclusions of our results obtained from our analysis are presented in section 7. In this paper, the cosmological model used is that of a flat Universe with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7.

2 Data

2.1 Sample selection

Our sample is a sub-sample of IR-bright DOGs revealed by Noboriguchi et al. (2019). The selection criteria for DOGs from Noboriguchi et al. (2019) (hereafter, DN19) was used in a similar manner to that of Toba et al. (2015). In brief, they are chosen by combining the Subaru HSC deep optical imaging data and the WISE all-sky data. However, because WISE’s spatial resolution is much lower than that of HSC, it is not easy to distinguish WISE sources from their optical counterparts. As a result, between all HSC-detected objects, the NIR imaging data was used to select carefully any relatively red sources, similar to the method used in Toba et al. (2015).

2.2 Used catalogs

In this study, we used the Herschel point source catalog (HPSC; Marton et al. 2015) beside catalogs which were already used in Noboriguchi et al. (2019): the HSC catalog data (optical), VIKING catalog data (NIR), and WISE catalog data (MIR). In figure 1, we show an example for three DOGs of our sample that have images at Herschel bands.

The HSC is a wide-field optical imaging camera with a field of vision that is 1.°5 in diameter. This camera was put in the main focus of the Subaru Telescope (Furusawa et al. 2017; Komiyama et al. 2017; Miyazaki et al. 2018). The HSC survey team produced the S16A catalog, which is based on data collected between 2014 March and 2016 April. Five different bands are used to monitor the S16A wide2’s total survey area of 456 deg². 178 deg² within the monitored region contains the anticipated full depth information for each of the five bands (g, r, i, z, and y; Aihara et al. 2017).

Utilizing the VISTA telescope (Dalton et al. 2006) and the VISTA Infrared Camera, the VIKING survey obtained image data over a large region in the NIR using five different bands (Z, Y, J, H, and K_s). The research conducted by Noboriguchi et al. (2019) made use of the DR2 catalog.² The increased assessment area in VIKING DR2 results in a larger number of items being found in the DR2 data. This is due to the fact that the survey depths of VIKING DR1 and DR2 are comparable. The magnitudes in the VIKING NIR have also been adjusted to account for galactic extinction (Schlegel et al. 1998).

The MIR bands (3.4, 4.6, 12, and 22 μm) of all-sky images were obtained by the WISE satellite. The ALLWISE catalog (Wright et al. 2010; Cutri et al. 2014) was used in this study.

The HPSC (Marton et al. 2015) was constructed from all 6878 usable SPIRE scan map observations in a standard configuration including all applicable calibration observations. The total of 1693718 objects split via filter band into 950688, 524734, and 218296 objects for 250 μm, 350 μm, and 500 μm, respectively. These totals also contain sizeable percentages of somewhat extended sources; these are (in order of increasing wavelength) |$29\%$|⁠, |$38\%$|⁠, and |$45\%$|⁠. A comparison of the source extractions of SPIRE's prime calibrator Neptune with model data demonstrates that the catalog fluxes are highly compatible with the radiation model used for the calibration. The flux values of four different photometric methods are described to provide good guidance to scientists: two are exclusively applicable to point sources (TML, Sussextractor; Savage & Oliver 2007), one is to classify potentially extended sources and as a sanity check (DAOPHOT; Stetson 1987), and one is to derive a best guess for extended sources, assuming a Gaussian elliptical profile. Such values are largely in agreement for point-like sources, and when they deviate from each other they will lead the astronomer to potential source extension (Schulz et al. 2017). The products explained in this catalog are made available online through the Infrared Science Archive (IRSA)³ and the Herschel Science Archive (HSA).⁴

2.3 The selection criteria of the DN19

In the process of choosing the DN19, clean samples were made, each of which consisted of only correctly identified objects. The number of objects that can be found in the HSC S16A wide2 forced catalog is 293520279. Objects that had photometry that had been incorrectly measured were rejected from participation in the sample because of the way that Toba et al. (2015) chose its participants. In order to get rid of the items in the sample that had inconsistent photometry, the ones that had a signal-to-noise ratio (S/N) <5 in any of the bands were thrown out. Objects that lacked enough exposure in any of the five HSC bands were rejected from inclusion in the sample. As a direct result of this, the clean sample for the HSC consisted of 16680947 items.

The number of objects found in the VIKING DR2 catalog is 46270162. By using the criteria “primary_source = 1,” which was used to produce the VIKING clean sample, it was possible to get rid of objects that are not one-of-a-kind. Using the criteria “kspperrbits = 0,” objects that were influenced by substantial noises were removed from the clean sample of the VIKING. In addition, objects that had a S/N < 5 in the K_s band were eliminated. As a direct result, 13455180 objects served as the clean sample for the VIKING.

The total number of objects in the ALLWISE catalog is 747634026. By using the criteria “w4sat=0” and “w4cc_map=0,” the ALLWISE clean sample does not include any objects with inadequate photometry due to significant noise. Objects in the W4 band with an S/N ratio of less than 3 were likewise omitted from the ALLWISE clean sample. Consequently, 9439990 objects were left as the ALLWISE clean sample.

Cross-matching was done between the clean samples to select IR-bright DOGs. The overlap region is ∼105 deg² in total. The selection criteria for DOGS that were decided upon were as follows: (i − [22])_AB ≥ 7.0. Thus, the DN19 contains 571 IR-bright DOGs (see figure 1 in Noboriguchi et al. 2019).

3 Methodology

Our primary objective is to select and analyze a representative sample of DOGs so that we may better understand the fundamental aspects of their physical and statistical features. To accomplish this objective, it was essential to take into consideration the characteristics of these galaxies throughout a wide range of wavelengths. By fitting their SEDs, taking into account photometric information ranging from the visible to the FIR, we were able to identify the primary underlying physical features of IR-bright galaxies.

In this study, the criteria of the DOGs’ selection is to select the first 28 brightest DOGs at 22 μm which have AB magnitudes between around 10 mag and 14 mag and have Herschel data points at 250, 350, or 500 μm. The mean value of redshift for our DOGs is equal to 1.14 which implies DOGs are not objects from the local Universe.

The SEDs cover the rest-frame wavelength ranges from FUV to FIR (from 0.15 to 106 μm). They have been modeled using version 11.0 of the CIGALE software package (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019).⁵ When constructing the SEDS, CIGALE makes sure to strike a balance between the amount of energy absorbed by dust and the amount that is re-emitted in the IR. Figure 2 shows an example of the SED of DOG17 in our sample.

Considering the nature of the objects we investigate in this paper, the most essential finding is that CIGALE can model both star formation and AGN-related emissions efficiently and reliably. Furthermore, CIGALE includes a very valuable methodology for testing the validity of the parameters and their potential degeneracies by creating a mock catalogue that is analyzed in a similar method to the data that has been observed. For a more comprehensive description of the code, we recommend readers visit the CIGALE web page.

Based on the above-mentioned features of CIGALE, we utilized it to figure out our DOGs’ physical characteristics, such as stellar masses, SFRs, dust attenuation, IR luminosities, dust masses, AGN fraction, and others.

The SED fits of 28 DOGs are shown in figures 23–27 in appendix  1, and the list of used modules and input parameters of CIGALE are shown in table 2 in appendix  2.

The reduced chi-squared (χ²) value of the best model is used to determine the quality of the fitted SEDs, marginalized over all parameters. The probability distribution function (PDF) values of the produced concerned parameters (e.g., L_dust, SFR, M_*, IR power-law slope, the AGN torus angle with respect to the line of sight, and AGN fraction) are calculated in this paper. The values of the reduced χ² for our DOGs are shown in figure 3.

In order to construct the mock catalog, the best fit for each source is considered. After that, CIGALE is able to adjust each variable by including a value that is drawn from a Gaussian distribution with the same standard deviation as the uncertainty in the observation.

The analysis of this mock catalogue is then carried out in precisely the same manner as that of the primary observations. It is possible to make accurate estimates of those physical characteristics whose precise and estimated values are comparable to one another. A similar reliability check has been carried out by many works (e.g., Buat et al. 2011; Małek et al. 2014; Toba et al. 2019, 2020, 2022). The comparison of the findings of the mock with the actual catalogs in figure 4 reveals that CIGALE provides a very clear approximation of the crucial physical parameters that we use for our analysis, since the Pearson product-moment correlation coefficient (r) values are greater than 0.5.

4 The main physical properties of DOGs obtained from the SED fitting

The main physical parameters of our DOG sample are derived based on the CIGALE SED-fitting code shown in appendix  3 (tables 3, 4, and 5) and their analysis is presented in the following subsections.

4.1 Dust luminosity

By using the Draine and Li (2007) model, we were able to determine the dust luminosity of our DOGs. The median value of the dust luminosity is equal to 5.1 × 10¹² L|${}_\odot $|⁠. One peak can be identified in the distribution of dust luminosity, which can be seen in the upper center panel in figure 5. This peak indicates a relatively large sample of galaxies that have high-z and star-forming rate galaxies, with the greatest point centered at approximately log (L_dust/L|${}_\odot $|⁠) ∼ 12.7.

4.2 Dust attenuation

We have computed the amount of attenuation caused by dust in order to investigate how the highly dusty component influences the stellar portion of the spectra. Using CIGALE, we can determine the degree of obscuration of stellar luminosity for the stellar population utilizing V and FUV filter data. The dust attenuation for the star population of our sample was found to be relatively strong compared with normal galaxies (∼0.7 mag), LIRGs (∼1.6 mag), and ULIRGs (∼6.0 mag) (Małek et al. 2017), as shown in the last two plots of figure 5.

4.3 Star formation rate

Historically, IR and sub-mmm studies of chosen galaxies have concentrated on utilizing the re-radiated energy from dust at 8–1000 μm as a proxy for SFR. The foundation of this relationship is clearly explained in the pioneering study of Kennicutt (1998). Dust luminosity is a good tracer of SFR, the reason for that is the vast majority of star formation being obscured, and dust emission is formed by the absorption of photons emitted by massive stars. It is interesting to check the relation of the SFR with dust luminosity produced by CIGALE findings because CIGALE seeks to account for both obscured (radiated in the FIR) and unobscured (radiated in the UV) star formation (see figure 6). Additionally, CIGALE accounts for that portion of dust luminosity which is heated by earlier stellar populations. Our results in equation (1) are marginally consistent with those of Kennicutt (1998) in equation (2) with 2σ confidence.

In figure 8, we compare the SFR range for five DOGs’ samples at different redshifts with our sample. Extreme SFRs of the order of 10⁴  M|${}_\odot $| yr⁻¹ would be needed to power the luminosities seen in hot DOGs (Eisenhardt et al. 2012; Assef et al. 2015).

4.4 Stellar mass

The computed mean value of log (M_*) for our DOGs is equal to 11.49 ± 1.61 M|${}_\odot $|⁠. The stellar masses of our DOG sample obtained by CIGALE are consistent with stellar masses calculated by Małek et al. (2017) (in their table 2) for ULIRGs from an ADF-S sample and located at similar redshifts. However, our DOGs are more massive than other DOG samples by factors of 9, 4, and 8, respectively, with respect to the following studies selected from infrared surveys and published by:

• Bussmann et al. (2009b) (in their table 5); for 12 DOGs in the Boötes Field of the NOAO Deep Wide-Field Survey at z ≈ 2 it is shown that the average of their log (M_*) is 10.53 M|${}_\odot $|⁠.

• Bussmann et al. (2009a) (in their table 6); 31 DOGs (24 μm bright) were studied at z ≈ 2 identified in the Boötes Field of the NOAO Deep Wide-Field Survey. We found that the mean value of their sample of log (M_*) is 10.93 M|${}_\odot $|⁠.

• Toba et al. (2017) (in their table 1); seven DOGs were discovered with the IRAS faint source catalog and the Akari FIR all-sky survey bright source catalog at z < 1. These DOGs have a mean stellar mass of log (M_*) = 10.59 M|${}_\odot $|⁠.

The M_* of our sample as a function of L_IR is shown in figure 9.

4.5 IR power-law slope

If photometric points at rest-frame λ ≤ 70 μm are available, the MIR power-law slope (α) may be restricted, and alpha can be determined. Otherwise, a fixed value of α = 2.0 is in agreement with the vast majority of sources and is analogous to the MIR component of SED templates (Casey 2012). The typical values of alpha vary from 0.5 to 5.5. However, there is an extra cut-off that must be set on α values at short wavelengths to prevent a non-physical explanation due to the divergent luminosity in the NIR at α < 1.

Figure 5 shows the α parameter values for our DOGs. The majority of our resultant α values are between 2.0 and 3.0. A weak linear correlation is seen in figure 10.

4.6 Fractional contribution of AGN IR emission and types of AGNs

As the part that AGNs play in DOGs is not totally understood, one of the primary goals of our study was to assess the AGN contribution to the DOGS ’ IR emission and determine the types of AGNs that are associated with the activity of DOGS. Using the Fritz, Franceschini, and Hatziminaoglou (2006) templates, we’ve calculated the AGN IR emission’s fractional contribution using the central source’s point-like isotropic emission and dust radiation with a toroidal shape proximity the central engine. The toroidal obscurer absorbs the AGN emission, where it is either re-emitted at wavelengths between 1–1000 μm or dispersed by the obscurer itself.

Figure 5 shows the AGN fractional contribution to our sample of DOGs (the middle figure in the second column). We find that the AGN contribution is higher than |$20\%$| for 20 DOGs (about |$71\%$| of the sample). According to Ciesla et al. (2015), the AGN emission fractional contribution to the L_dust bounded by CIGALE is always overestimated for a low proportion of Type 1 AGNs (AGN_frac < 0.1), while it is underestimated for Type 2 and intermediate types of AGNs (AGN_frac > 0.2). Because the majority of our DOGs had an AGN_frac value greater than 0.2, we have concluded that the actual contribution is likely underestimated, and as a result, the values that we have obtained may be regarded as lower limits.

The mean value of the AGN fraction of our DOGs is 0.32, and the brightest galaxy at 22 μm, HSC J140638.20+010254.6 (ID = DOG_1), has 0.69 as a maximum AGN fraction among the whole sample. In the study of Gabányi, Frey, and Perger (2021) investigating the radio properties of a large sample of 661 DOGs, this is among the very few radio-detected sources, suggesting that the radio emission, at least partly, is attributed to AGN activity. Interestingly, for the other two radio-detected DOGs in the Gabányi, Frey, and Perger (2021) sample that we studied, HSC J140738.47+002731.4 (DOG_5) and HSC J141546.89−011451.7 (DOG_7), we also found relatively large AGN_frac values; 0.59 and 0.47, respectively (see also Toba et al. 2017). The model of Fritz, Franceschini, and Hatziminaoglou (2006) describes different types of AGN model templates, including a Type 1 AGN (known as an unobscured AGN), a Type 2 AGN (also known as an obscured AGN), and a template that sits in the middle of them and is known as an intermediate type. According to Fritz, Franceschini, and Hatziminaoglou (2006) (also Mountrichas et al. 2021), the torus angle of an AGN relative to the line of sight is defined as: ψ = 0.°001 for Type 2 AGN and ψ = 89.°9 for Type 1. As shown in figure 5, |$80\%$| of the sample have ψ < 5° and the rest is between 60° and 90°. Following the analysis of the findings acquired from the AGN component, we have concluded that there are many distinct types of AGN components for DOGs. AGNs of Type 1 and Type 2 and intermediate types are included in our sample; nevertheless, Type 2 AGNs are responsible for the majority of the fractional contribution to the total amount of AGN emission. This shows that Type 1 AGNS produce only a minor contribution to the MIR emission for galaxies with high dust luminosities, and the high brightness of the dust component of the spectrum is connected to their star-forming activity.

Figure 11 illustrates the link between the torus angle with regard to the line of sight and the AGN fractions that contribute to the MIR emission of the sample. The prevalent pattern indicates that DOGs that include Type 2 AGNs are likely to have a larger AGN contribution.

4.7 Metal abundance

In this study, the SSP library we utilized (Bruzual & Charlot 2003) can be used with a wide range metallicities (Z) in discrete values of (0.0001, 0.0004, 0.004, 0.008, 0.02, and 0.05). Tremonti et al. (2004) examines 50000 galaxies from the Sloan Digital Sky Survey (SDSS) to explore a tight mass–metallicity (MZ) relation in the local Universe. Surveys of distant galaxies have allowed investigating of the MZ relation at intermediate (Savaglio et al. 2005; Maier et al. 2005; Zahid et al. 2011; Pérez-Montero et al. 2013) and high redshifts (Erb et al. 2006; Maiolino et al. 2008; Mannucci et al. 2009; Laskar et al. 2011; Kulas et al. 2013).

The association that galaxies are less enriched at higher redshifts when held to the same stellar mass is still contested. Investigations of the chemical history of galaxies give substantial constraints to the star formation processes and gas fluxes in galaxy evolution models (e.g., Brooks et al. 2007; Davé et al. 2011; Zahid et al. 2012, 2014; Torrey et al. 2014). Although some studies indicate that as the universe evolves, the gas in galaxies will grow increasingly metal-rich, other researchers have stated that this relation will be flat for massive galaxies at late periods of time (Savaglio et al. 2005; Maier et al. 2005; Maiolino et al. 2008; Zahid et al. 2011, 2013). This flattening may be attributed to galaxy downsizing (Cowie et al. 1996), which is defined as the process through which star formation becomes more prominent in lower-mass systems at late periods of time.

Figure 12 represents our MZ relation for all 28 DOGs with two other samples of galaxies that show a correlation in MZ relation: Savaglio et al. (2005) for 24 galaxies from the Gemini Deep Deep Survey (GDDS) at 0.4 < z < 1.0, and Maiolino et al. (2008) for nine galaxies of the AMAZE sub-sample at z ∼ 3. Our data points in figure 12 show a kind of random distribution; this could be due to our sample size not being big enough to determine the general trend of the MZ relation of our sources.

In addition, we do not detect any significant redshift evolution for the total sample of 28 DOGs, from z = 0.47 to 1.63, as shown in figure 13.

5 Comparing the DOGs to other galaxy samples

The finding of a link between SFR and the stellar mass is referred to as the “main sequence” (MS) of SF galaxies at z ∼ 0 (Brinchmann et al. 2004) to z = 7 (e.g., Stark et al. 2009) and throughout z ∼ 1 (Noeske et al. 2007; Elbaz et al. 2007), z ∼ 2 (Daddi et al. 2007; Pannella et al. 2009), and z = 3–4 (Daddi et al. 2009). The MS of SF galaxies is an essential stage in directing the galaxy evolution. Nevertheless, several studies conclude that the normalization, slope, and form of the MS are significantly different. These differences are mostly attributable to the varied selection criteria that were used to isolate SF galaxies and some galaxies with a certain SFR might be included or excluded from the analysis depending on the criteria. The strength of this relation is not compatible with frequent random bursts generated by processes like major mergers of gas-rich galaxies and thus promotes more stable star formation histories (Noeske et al. 2007). In addition, systematic examinations of the dust characteristics of the “average galaxy” at various redshifts demonstrate that LIRGs at z = 1 and ULIRGs at z = 2 are closely similar to conventional SF galaxies at z = 0.

The MS correlation outliers, which are often referred to as “starbursts,” exhibit indicators of varied dust characteristics, including the following: A geometry that is more compact (Rujopakarn et al. 2010), an excess of IR8 ≡ L_IR(8–1000 μm)/L_ν(8 μm) (Elbaz et al. 2011), and a shortage of [C ii] (Díaz-Santos et al. 2013). Furthermore, these starburst galaxies have been shown to have an elevated effective dust temperature (Elbaz et al. 2011) and a polycyclic aromatic hydrocarbon (PAH) deficit (Nordon et al. 2012), suggesting that they are real analogs of local LIRGs and ULIRGs. According to this paradigm, the features of these galaxies are not very tied to their rest-frame bolometric luminosities, but rather to their higher SFR compared to the MS (Schreiber et al. 2015).

In our study, we used the MS of the SF galaxies diagram of Schreiber et al. (2015) and located our sample on it, as appears in figure 14. Our sample is clearly shown to be located far above the MS of SF galaxies, which means they are consistent with starbursts.

The paper of Elbaz et al. (2011) presents the Herschel-derived FIR SEDs of SF galaxies that are located in the GOODS fields. It was discovered that the majority of them, including the z = 1–2 LIRGs and the ULIRGs, are thought to follow an IR MS which they define depending on the IR8 parameter.

In a predictable manner, PAH strength follows L_IR for the majority of SF galaxies in the local Universe, e.g., IR8 ∼ 4 (Elbaz et al. 2011). These typical SF galaxies exhibit a narrow range of specific SFRs and describe the IR MS (Elbaz et al. 2007; Noeske et al. 2007). On the other hand, PAH strength does not follow L_IR with an increase in IR8 for galaxies undergoing a rapid starburst. ULIRGs in the local Universe are often found outside of the IR MS, and their IR8 ≫ 4.

Calculating the IR8 values of our selected galaxies based on the Herschel observations that provide the best-fitting SEDs allows us to compare our DOGs directly with the other samples. The bulk of IR8 values of our PL DOGs lie in the 1 < IR8 < 8 range, with a mean value of 5.4. In comparison, the seven bump DOGs of our sample (DOG 6, 13, 16, 20, 26, 27, 28) have a broad variety of IR8 values with a high mean value ∼14.2 (see table 5).

These results are consistent with another sample of DOGs studied by Melbourne et al. (2012). Their PL DOGs have a narrow range of IR8 values with a mean IR8 ∼ 6, while their bump DOGs exhibit a broad range of IR8 values of 2 < IR8 < 23.

According to Elbaz et al. (2011), IR8 values rise when star formation is taking place in morphologically compact areas. These extremely compressed star-forming regions are generally the result of major mergers that move gas to the cores of these systems when they occur in the local ULIRGs. However, regarding the DOGs, whether the same merger-related procedures contribute to the high values of IR8 is not very clear. Even though the DOG samples obviously show some merging, the fractions with clear major merger signs are still quite small, less than |$30\%$| (Melbourne et al. 2009).

6 Statistical properties of the fitted parameters

Fitting the models to the observed spectra yielded many parameter values for each galaxy in the sample. Six of them were chosen for PCA, namely redshift, L_AGN [L|${}_\odot $|], L_dust [L|${}_\odot $|], M_dust [M|${}_\odot $|], M_* [M|${}_\odot $|], and FUV_att [mag]. The reason for this choice is that some other parameters have discrete values, and that is not allowed in PCA. In addition, introducing two very similar quantities in the analysis like FUV and V attenuation (also SFR and dust luminosity) can perhaps negatively impact the way that the PCA algorithm performs, therefore they were also excluded.

6.1 Performing PCA of our sample

Although the six obtained PCs reproduce the observed variables completely, their actual value may be determined by the physical parameters of a much smaller number. Therefore, it could happen that a smaller number of PCs is sufficient to reproduce the basic characteristics of SEDs. The magnitudes of λ eigenvalues give information on the importance of the given PC in reproducing the original observed variable. The λ eigenvalues give the variances (squared standard deviations) of the PCs. The direction of the first PC (PC1) in the 6D parameter space of the observed DOGs corresponds to the greatest variance in the point pattern representing the galaxies in the studied sample. We display the obtained eigenvalues in figure 15, in descending order.

If there is a correlation between the variables, the non-diagonal matrix elements differ from zero. Then the λ eigenvalues obtained as a solution of the eigenvalue equation already differ from each other, although their sum does not change, i.e., in our case it remains 6. As the sum of λ eigenvalues remains unchanged, some of them will exceed 1 and some will be smaller in the case where there are correlations between variables. The λ divided by the sum of the eigenvalues (in our case, by 6) gives the fraction a given λ represents in the sum of eigenvalues. According to the Kaiser (1960) criterion, PCs with eigenvalues greater than 1 contribute significantly to the observed variables. The dashed horizontal line in figure 15 shows the Kaiser criterion level. As one can see in figure 15, there are only two eigenvalues exceeding this level. It means there are only two PCs (PC1 and PC2) contributing significantly to the observed shape of the DOGs’ SEDs.

Dividing the eigenvalues with their sum, the total variance, we get what percentage of the total variance the given eigenvalue represents. These fractions are displayed in figure 16. Apparently, the two eigenvalues above the Kaiser limit (dashed line) explain |$72.7\%$| of the total variance. The eigenvectors belonging to these two eigenvalues stretch a 2D subspace in the 6D parameter space of the observed variables. The points representing the original SEDs in the 6D parameter space tend to concentrate on this 2D subspace. In the following, we discuss the statistical properties of the projected positions of the original points in this subspace.

The elements of the eigenvectors (called loadings) obtained from the PCA represent the correlations with the respective observed variables. Figure 17 shows a biplot of the correlations of the observed variables with PC1 and PC2 along with the position of the DOGs in this coordinate system. The arrows are directed to the point at which the coordinates are proportional to the correlations of the corresponding variable with PC1 and PC2, respectively. Accordingly, the vector length indicates the strength of correlation.

6.2 Results obtained from PCA

Figure 17 indicates that the observed variables, except for FUV_att, direct closely to Dim1 (PC1), the redshift vector has almost the same direction. The similar directions of the arrows, representing the observed variables, indicate a strong correlation between them. The length of the vector for each variable means the extent to which PC1 and PC2 reproduce the variable itself. The longest vector belongs to the variables L_dust, followed by the Redshift and FUV_att variables. The closest to PC1 is the Redshift variable. In this, the PC2 component is the smallest. However, due to the small sample size, this minimum value does not necessarily mean that this variable is actually more closely related to PC1 than other variables that significantly affect the intensity of FIR radiation (AGN luminosity and dust luminosity) and have a larger PC2 but not a significant component. (Non-significant values are indicated by a blank field in figure 18.) Figure 18 shows the correlations between the observed variables in color-coded form. The blue code represents positive correlations and the red code represents negative correlations. The size and darkness of circles indicate the strength of the correlations. The basic equation of PCA is the eigenvalue equation of the matrix made up of linear (Pearson) correlations between the variables. These linear correlations express the relationship between variables and also linear correlations yield their relationship with PCs. The linear correlation between variables is shown in figure 18.

Similar conclusions can be drawn from figure 19, which shows the correlation between the observed variables and the two strongest PCs and two more PCs in color-coded form. The color coding of the figure is similar to that of the previous one.

Figure 18 perhaps even more strikingly shows than figure 17 that four variables correlate strongly. In contrast, ultraviolet attenuation is not involved in this. A further peculiarity of this variable is that it is negatively correlated with the other variables, as shown in figure 18; it shows a significant negative correlation with L_AGN and Redshift variables. In subsection 7.2, we return to a possible explanation of this fact.

6.3 Partial correlation between observed variables

In the previous subsection, we were dealing with performing PCA on observed variables. PCA is an eigenvalue problem of the matrix obtained from the mutual correlations between variables obtained from the parameters used for fitting SEDs. When computing the mutual correlation of several variables, it may happen, however, that two of them show a correlation only because both of them correlate with a third one. The partial correlation is a procedure that simply examinines whether the correlation between the variables is maintained even if the effect of this third is removed.

In our case, the correlation between quantities describing IR radiation and masses may be the result of their correlation with redshift, and this correlation may disappear in the co-moving system. Given such a problem, we estimated the value of the correlation (the partial correlation) provided that some of it is given by their correlation with redshift.

To calculate the partial correlation, the pcor() procedure in the ppcor library of R was used. The results are shown in table 1. The comparison between the PCA and partial correlations for some relations is discussed in subsection 7.2.

6.4 Grouping of points in the parameter space

The classification of a sample into groups of similar objects requires a distance or (dis)similarity measurement among each pair of its components. A dissimilarity or distance matrix is the outcome of this computation.

The usual methods for distance measures are Euclidean (square root of the sum of squares of coordinate differences) and Manhattan (sum of absolute values of coordinate differences) distances. In both cases, the coordinate axes are perpendicular to each other. If the Euclidean distance is used, observations with high similarity values will be clustered together. The same is true for observations with low values of distances. The grouping was carried out in the {PC1, PC2} plane in our case to ensure that the coordinate axes were orthogonal.

We used the NbClust library in R for estimating the best number of clusters in our data set. This library includes 30 indices for determining the optimal number of clusters; it also proposes the use of the best clustering scheme based on the various results achieved by varying all combinations of a number of clusters, distance measures, and clustering methods. In the NbClust library, we used the Euclidean distance to compute the dissimilarity matrix, and “ward.D2” as the cluster analysis method.

The Ward method is an agglomeration procedure. To begin, all the sample elements are considered as individual clusters. During the agglomeration procedure, the elements are merged into larger clusters until the whole sample represents a big one. The whole agglomeration procedure can be displayed in a merger tree, the dendrogram, showing the levels where two clusters were merged into a bigger one. The Ward method helps to reduce the overall variation that exists inside the cluster. For every stage, a couple of clusters with the shortest cluster distance are combined. In order to put this strategy into practice, at each stage we look for the pair of clusters that result in the least overall increase in the total within-cluster variation after being merged. The resulted dendrogram is shown in figure 20.

According to the result of the nbclust procedure of the NbClust library, there are four clusters which are represented in four red boxes in figure 20, showing the dendrogram of the data clustering. Each box contains the row number of our DOGs. With the help of discriminant analysis, we can find the directions along which the four obtained clusters are best separated. Discriminant analysis was performed using the lda procedure in the R’s mass library. The separation of the four clusters is the strongest along with the two obtained orthogonal directions (Component 1 and Component 2; figure 21).

Now, we can obtain the matrix plot indicating the four obtained clusters as shown in figure 22. Apparently, the redshift, L_dust, and M_* separate cluster 4 from clusters 1, 2, and 3 best. This is more or less true for the other variables that correlate strongly with it. The clusters are sharply separated in the variable L_dust with M_* and in the variable FUV_att with other parameters. Cluster 1 is between other clusters. Consequently, the light attenuation and AGN luminosity closely relate to variable Component 1, perpendicular to the variable Component 2.

Due to the larger redshift, cluster 2 represents the earlier state of the Universe versus clusters 3, 1 and 4, respectively. The observed bifurcation in light attenuation appears at a later age of the Universe. In our sample, light attenuation does not correlate with the dust luminosity and stellar mass, and does so weakly with the other two variables. The strength of attenuation depends on the amount of interstellar dust in the line of sight of the stellar component of DOGs. However, it depends much less on the absolute mass of the dust responsible for the IR radiation.

Accordingly, the observed bifurcation in light attenuation is a geometric effect dependent on the angle of dust distribution to the line of sight. In an earlier age of the Universe, the difference in the flattening of the stellar mass and of the dust is much less pronounced. At a later age, however, it could be much more advanced and resulted in the observed bifurcation of the UV light attenuation. Let us assume a strongly flattened dense dust layer, embedded in a much less flattened stellar component. When the observer sees this dust layer edge-on, the attenuation is minimal (except for a thin visible layer, the stellar component is little attenuated) and it is maximal when it is face-on to the observer (about half of the stars are then heavily attenuated). The magnitude of the effect can be estimated with a simple model. In our Galaxy, the disk has a scale height of 300 pc, while that of dust is 100 pc. If we see the galaxy edge-on, the dust covers about |$10\%$| of the radiating stars. As soon as we see the Galaxy face-on, the radiation of that part of the stars facing us is not attenuated as the dust disk is behind them. The other part, on the other hand, is placed behind the disc and the degree of attenuation depends on the column density of the light-scattering particles that make up the dust disc. Thus, the degree of light attenuation, even with the same physical properties, depends on the angle of the galaxy with respect to the line of sight.

7 Summary and conclusions

We analyzed the physical and statistical properties of a sample of 28 bright IR selected DOGs detected at 22 μm. These sources were discovered by Noboriguchi et al. (2019) by combining the WISE all-sky data and the deep optical imaging data obtained with the Subaru HSC. We have added sub-mm data points from the Herschel SPIRE point source catalog (FIR).

The software CIGALE v.11.0 (Burgarella et al. 2005; Noll et al. 2009; Roehlly et al. 2011) was utilized for fitting spectral energy distribution in order to estimate the fundamentally important physical characteristics of the DOGs. These fundamentally important physical characteristics include SFR, stellar mass, dust luminosity, AGN fraction of IR emission, and the AGN’s viewing angle with respect to the line of sight. All 28 galaxies of our sample have SEDs fitted with photometric redshifts. The accuracy of the fits was satisfied, and the physical parameters found were thought to be reliable.

By PCA analysis, the six dimensions of variables are reduced to just the two most important dimensions. The redshift, dust luminosity, stellar and dust masses and AGN luminosity influence each other, while do not influence FUV attenuation, and vice versa. Then the partial correlation between observed variables is used to test the real physical correlations without the influence of redshift.

7.1 Notes on the derived physical properties

Our DOG sample shows very high values in stellar masses, log (M_*/M|${}_\odot $|⁠) = 11.49 ± 1.61, compared with other samples of DOGs, as shown in subsection 4.4. We could say that our galaxies are “overweight” DOGs (hereinafter ODOGs). Bussmann et al. (2012) found that PL DOGs and bump DOGs are typically 2 and 1.5 times more massive than SMGs, respectively (median and interquartile values of M_* for SMGs, bump DOGs, and PL DOGs are |$\log (M_{*}/M_{\odot }) = 10.42^{+0.42}_{-0.36}$|⁠, |$10.62^{+0.36}_{-0.32}$|⁠, and |$10.71^{+0.40}_{-0.34}$|⁠, respectively). During a major galaxy merger, it is believed that DOGs are transitioning between the starburst-dominated and AGN-dominated phases. According to this scenario, SMGs evolve into DOGs, then (optical) quasars, and at the end become “red and dead” elliptical galaxies (e.g. Narayanan et al. 2010; Bussmann et al. 2012).

Based on Herschel data, Riguccini et al. (2015) showed that faint DOGs (with 8 μm luminosity below 10¹² L|${}_\odot $|⁠) are dominated by star formation, while the AGN activity is typical for brighter DOGs. Our results in subsection 4.1 show that the luminosities in our DOG sample are larger than 10¹² L|${}_\odot $|⁠. Moreover, subsection 4.6 clearly presents that |$71\%$| of our DOGs are Type 2 AGNs. This indicates a large contribution to the emission of MIR for galaxies with high dust luminosities, and the high luminosity of the dust part of the spectra is related to AGN activity.

DOGs look like a transitional point between the dusty starburst-heated and the AGN-dominated phases of the massive galaxies evolution (e.g., Bussmann et al. 2012; Riguccini et al. 2015), where AGN feedback is completely active. The same holds true for the tiny subclass of massive, enormously luminous DOGs. Discovering a sub-class of DOG objects like ODOGs is promising. ODOGs may be an indication that the DOG phase has come to a close, after that they transform into visible quasars.

7.2 Notes on the statistical properties

The PCA analysis in this work shows a negative correlation with L_AGN and redshift in subsection 6.2. The magnitude of the light attenuation depends not only on the dust mass but also on the average size of the light-scattering dust particles. Interpreting the dependence of the attenuation effect on redshift as an evolutionary effect over time, we find that although the total dust mass decreases, the frequency of small, visible, ultraviolet light-scattering dust particles increases. The negative correlation with L_AGN can be explained by the fact that small particles, active in light scattering, are destroyed during activity and attenuation decreases.

It is important to note that redshift is among the three highly correlated variables, although it is independent of the DOGs’ internal physical conditions. It would be trivial to interpret this fact as a selection effect, as all observed samples are limited by the apparent brightness detection and require greater absolute brightness at a larger redshift. Provided that the DOGs’ absolute brightness distribution is independent of redshift and correlates with the stellar mass of the galaxy.

Comparing the (dust mass–dust luminosity) correlation values in table 1 with the correlation between the variables given in figure 18 by PCA, it can be realised that they no longer correlate with each other. This means that this correlation is just a common background variable relationship. On the other hand, (stellar mass–dust luminosity) and (attenuation–AGN luminosity) correlations remained strong, meaning that they have a real physical relationship. That confirms our correlation between the dust luminosity and the stellar mass in figure 9.

Acknowledgements

We are grateful for the expert assistance of the staff of the Astronomy Department at Eötvös Loránd University and Konkoly Observatory in Hungary. In particular, we thank Krisztina Perger and Sándor Pintér. This work is based in part on observations made with the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center at the National Astronomical Observatory of Japan. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication has made use of data from the VIKING survey from VISTA at the ESO Paranal Observatory, program ID 179.A-2004. Data processing has been contributed by the VISTA Data Flow System at CASU, Cambridge and WFAU, Edinburgh. This research has made use of the HPSC which is available on NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work has been supported by the Stipendium Hungaricum scholarship by the Tempus Public Foundation, the JSPS KAKENHI grant Nos. 18J01050 and 19K14759 (Y.T.), and the Hungarian National Research, Development and Innovation Office grant No. OTKA K134213 (S.F.). Facilities: HSC, VISTA VIKING, WISE, Herschel (SPIRE). Software: CIGALE (Noll et al. 2009), Python 3 (Van Rossum & Drake 2009).

Funding

The publishing of this paper was funded by the Hungarian institutions and the Electronic Information Service National Program (EISZ).

Conflict of interest

All authors declare that they have no conflicts of interest.

Appendix 1. The SED fits

This appendix contains the 28 SED fits for our sample that are resulted by using CIGALE, they are presented in figures 23–27.

Appendix 2. The input parameters

The input parameters of the used modules in the CIGALE code are shown in table 2.

Appendix 3. The physical properties

Table 3, 4, and 5 show the physical parameters of our sample that are obtained based on the CIGALE SED-fitting code.

Footnotes

References