Local Analogs to High-redshift Galaxies. I. Characterization of Dust Emission and Star Formation History

The Brown et al. (2014) templates are used to fit to the photometry of high-redshift galaxies from CANDELS. CANDELS is a treasury program on Hubble Space Telescope (HST) and provides deep multiwavelength imaging data in five legacy deep fields: GOODS-South, GOODS-North, UDS, COSMOS, and EGS (Grogin et al. 2011; Koekemoer et al. 2011). In this study we use H-band-selected catalogs from all five fields. Details on the data and photometry at different wavelength bands are given for GOODS-S (Guo et al. 2013; Santini et al. 2015), GOODS-N (Barro et al. 2019), UDS (Galametz et al. 2013; Santini et al. 2015), COSMOS (Nayyeri et al. 2017), and EGS (Stefanon et al. 2017).

The CANDELS catalogs contain photometric data from the UV to NIR and MIR wavelengths, observed in broadband, medium-band, and, in some cases, narrowband filters. The catalogs also contain derived physical parameters, such as photometric redshift, stellar mass, and SFRs, compiled using several different teams applying different SED methods, SFHs, and templates (e.g., Dahlen et al. 2013; Mobasher et al. 2015; Santini et al. 2015). The optical (HST/ACS) and NIR (HST/WFC3) data are consistently combined with the MIR data (Spitzer/IRAC) and ground-based observations (UV and K band). All of the photometric catalogs were selected in the HST/WFC3 F160W band using SExtractor (Bertin & Arnouts 1996). The catalogs for the five different CANDELS fields contain a total of 185,956 galaxies. The number of photometric bands ranges from 18 (GOODS-S) to 42 (CANDELS COSMOS). The latter catalog includes extensive narrow- and medium-band photometric bands. The MIR data are from Spitzer/IRAC (Ashby et al. 2013).

All 129 of the Brown et al. (2014) templates are used in the SED fitting of the 185,956 CANDELS galaxies. The fits are done in two different modes: keeping the extinction equal to the intrinsic dust extinction of the templates, and allowing additional extinction on top of the intrinsic. The additional extinction is parameterized through the E(B − V) color index. When treated as a free parameter, we allowed additional dust extinction up to E(B − V) = 0.6. We use the Calzetti et al. (2000) attenuation law for the additional extinction. We do not include "negative extinction" in the SED fitting. The Brown et al. (2014) sample contains galaxies with a wide range of metallicities and dust content and should be representative for galaxies with metallicities down to the levels expected for moderately massive galaxies at high redshifts. As shown below, the results in terms of morphology fraction seem to be insensitive to the inclusion of additional extinction (Figure 1). An initial test was done keeping the photometric redshift as a free parameter. The redshifts from the SED fits were similar to those presented in the CANDELS catalogs. For the reminder of the analysis, we keep the redshift fixed at the CANDELS values, which includes spectroscopic redshifts whenever available.

Zoom In Zoom Out Reset image size

The SED fitting procedure involved redshifting the template to the redshift of each galaxy, and then each galaxy template is convolved with the filter response functions appropriate for a given CANDELS data set. The template is then fitted to the observed flux densities using χ² minimization. If additional extinction is included, we repeat this process for E(B − V) values ranging from 0.0 to 0.6, in steps of 0.01.

We divide the CANDELS galaxies into 13 redshift bins between 0.3 ≤ z ≤ 5 (Δz = 0.36). The total number of galaxies in each redshift bin ranges from ∼2000 at z ∼ 0.3 to ∼10,000 at z ∼ 2. For each redshift bin we derive the number of best-fit SEDs for each morphological group as defined above. In Figure 1 we plot the fraction of best-fit SEDs for the four morphological classes as a function of redshift. From the figure it is clear that the irregular class (6 ≤ T ≤ 11) provides about 35% of the best-fit templates at z ∼ 0.3 and increases to dominate at z ≳ 1. At z ≳ 2, these templates account for 85%–90% of all best-fit templates. Of the 20 templates in the 6 ≤ T ≤11 morphological class, 11 of them provide the most frequent best-fit SED (dotted line in Figure 1). These 11 galaxies are selected as our sample of potential local analog candidates.

In order to test the robustness of the morphological fractions shown in Figure 1, we selected the template galaxy with ${\rm{\Delta }}{\chi }_{\nu }^{2}=3$ relative to the best-fit SED. The resulting morphological fraction as a function of redshift is shown in Figure 2. At z ≳ 1.5, the fraction of the 6 ≤ T ≤ 11 morphological class decreases slightly, while the fraction of star-forming spiral disks increases. However, the overall dominance of irregular galaxies as the best templates for high-redshift galaxies remains. This shows that the result is robust, and local templates with T-class morphologies corresponding to irregular star-forming galaxies are best at fitting the SED of high-redshift galaxies at z ≳ 1.

Zoom In Zoom Out Reset image size

This result is not surprising, as galaxies in the 6 ≤ T ≤11 morphological class share many properties with star-forming galaxies at high redshift (e.g., Elmegreen et al. 2009, 2012). The general properties of the 11 candidates are summarized in Table 1; they have small physical sizes (1–2 kpc), they are actively forming stars, and they typically have a high gas mass fraction. Their gas-phase metallicity is low, with an average metallicity of 12 + logO/H = 8.11 ± 0.20, or ∼0.26 Z_⊙. In Figure 3 we show the Sloan Digital Sky Survey (SDSS) i-band images of the 11 candidates. From these images it is clear that they have irregular morphology with prominent star-forming clumps.

Table 1. Properties of the Local Analogs

No.	Galaxy	R.A.	Decl.	D	$\mathrm{log}{M}_{* }$ a	$\mathrm{log}{M}_{{\rm{H}}\,{\rm\small{I}}}$ b	fgas c	Metallicity d	Alternative
	Name	J2000.0		(Mpc)	(M⊙)	(M⊙)	M⊙	$12+\mathrm{log}({\rm{O}}/{\rm{H}})$	Name
1	NGC 2537	08:13:14.4	+45:59:13	8.6	9.02	8.43	0.17	8.19	Arp 6; "Bear Paw"
2	Mrk 140	10:16:28.3	+45:19:18	27.6	8.60	8.90	0.56	8.30	Mrk 140
3	Haro 02	10:32:31.9	+54:24:02	23.7	9.21	8.73	0.25	8.45	Mrk 033
4	Mrk 1450	11:38:35.6	+57:52:27	14.7	7.27	7.35	0.54	7.96
5	UM 461	11:51:33.1	−02:22:22	20.7	7.37	8.47	0.96	7.78
6	Mrk 1307	11:52:37.4	−02:28:09	21.0	8.07	8.73	0.81	7.96	UM 462
7	Haro 06	12:15:18.4	+05:45:39	35.1	8.58	8.83	0.69	8.18
8	NGC 4670	12:45:17.1	+27:07:31	20.0	9.38	9.02	0.37	8.30	Arp 163
9	Mrk 475	14:39:05.5	+36:48:21	10.9	6.95	6.62	0.56	7.93
10	UGCA 410	15:37:04.2	+55:15:48	10.5	7.26	7.58	0.57	8.10	Mrk 487
11	Mrk 930	23:31:58.6	+28:56:50	77.5	8.89	9.51	0.74	8.08

Notes.

^a Stellar mass M_*; this work. ^b H i data from Paturel et al. (2003), except for Mrk 1450 (van Driel et al. 2016) and Haro 02 (Bravo-Alfaro et al. 2004). ^c Gas fraction is defined as f_g = M_g /(M_g + M_*). ^d Metal Abundances UM 461 and UM 462 (Campos-Aguilar et al. 1993), Haro 02 (Davidge 1989).

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

All 11 galaxies fall under the classification of BCDGs (e.g., Sargent & Searle 1970; Kunth et al. 1988; Thuan & Izotov 2005). BCDGs have previously been proposed to be analogs of high-redshift galaxies (e.g., Hoopes et al. 2007). As a group they share many properties, such as small sizes, clumpy, and intense star formation. They show strong optical emission lines, often exhibiting high excitation (Thuan & Izotov 2005). There are, however, differences among the BCDGs; some show the star-forming regions superposed on faint extended stellar emission, suggesting star formation extended over longer time periods (e.g., Noeske et al. 2003), while others seem to be truly young systems (e.g., Papaderos et al. 2008; Hunt et al. 2014). The gas-phase metallicities can range from extremely low, like SBS 0335−052E and I Zw 18 with Z ≈ 0.03 Z_⊙ (Izotov et al. 2005), to Z ≈ 0.3 Z_⊙.

One feature that most BCDGs share is the presence of dust. Even the low-metallicity systems SBS 0335−052E and I Zw 18 show thermal dust emission, even when detectable CO emission is absent (Hunt et al. 2014). Another feature shared among the BCDGs is that they are gas-rich. All of the galaxies discussed in this paper have neutral hydrogen gas with masses similar to or exceeding their stellar mass. Their gas fractions range from 0.17 to 0.96 (Table 1). The H i data are mostly based on single-dish observations. When higher angular resolution observations are available, it is found that the H i usually extends beyond the stellar component of the galaxies. Despite several dedicated surveys, only a small number of BCDGs have a detectable molecular gas component (e.g., Sage et al. 1992; Barone et al. 2000; Meier et al. 2001; Bravo-Alfaro et al. 2004; Israel 2005; Beck et al. 2020). The molecular gas mass is typically just a few percent of the total H i mass. One exception is Haro 02, where ${M}_{{{\rm{H}}}_{2}}/{M}_{{\rm{H}}\,{\rm\small{I}}}\sim 0.2$.

In the reminder of this paper we will characterize the dust emission and the SFH of the 11 local analog candidates defined above, to see whether they can be viewed as scaled-down versions of high-redshift galaxies.

SOFIA (Young et al. 2012; Temi et al. 2014) is a partnership between NASA and the German Aerospace Center (DLR), consisting of an extensively modified Boeing 747SP aircraft carrying a 2.7 m (106-inch) reflecting telescope (with an effective diameter of 2.5 m, or 100 inches). We were granted ∼8 hr of observing time with SOFIA, on Cycle 6 under the program 06-0222 (PI: T. Wiklind). We observed five galaxies of our sample, using the High-resolution Airborne Wideband Camera-plus instrument (HAWC+; Vaillancourt et al. 2007; Dowell et al. 2010; Harper et al. 2018), covering all available filters to derive the FIR SED.

The HAWC+ instrument has five bands identified as A, B, C, D, and E. We observed in bands A, C, D, and E at 53, 89, 154, and 214 μm, which have beam sizes (FWHM) of 485, 780, 136, and 182 respectively. Band B was not available during Cycle 6. The galaxies Mrk 1450, NGC 4670, and Haro 06 were observed in all four bands available. Due to time constraints, Mrk 1307 and NGC 2537 were only observed at 53, 89, and 155 μm bands. The data were reduced to obtain the total intensity maps using the HAWC+ Data Reduction Pipeline v1.3.0 as detailed in Harper et al. (2018).

3.1.1. SOFIA/HAWC+Photometry

The HAWC+ data are presented in five levels: level 0 is raw data, level 1 is the demodulated (chop not subtracted) data, after instrument corrections (flat-fielded, dark subtracted, bad pixel removed) becomes level 2, level 3 is obtained after flux calibration and telluric correction, and the combined observations from level 3 produce the level 4 data. More details on the pipeline process can be found in the data handbooks (Gordon et al. 2018).

We started our analysis with level 4 maps in units of Jy pixel⁻¹ and performed annular photometry using the HIPE environment. A circular aperture, centered on the source emission, gives the flux + background emission; an annular region, free of emission and centered on the inner aperture, provides an estimate of the background. To determine the photometric errors, we placed six circular apertures of the same radius as the inner source aperture, at equidistant, flux-free areas around the source. The average flux of the six apertures (Bg_ave) is then subtracted from the source flux:

$$\begin{eqnarray}&&{F}_{\mathrm{source}}=({f}_{\mathrm{tot}})-{\mathrm{Bg}}_{\mathrm{ave}}.\end{eqnarray} \tag{ 1 }$$

As a further test, we estimated the integrated flux using Casa ³ and the built-in two-dimensional Gaussian procedure. This produced integrated fluxes very similar to the ones obtained using aperture photometry.

To determine the photometry errors, we used six photometry apertures around the source, taking care to avoid any other source in the background, and then uncertainty of the measurement was defined as the standard deviation of the background apertures. The total error is computed as the addition on quadrature of the uncertainty of the measurement Bg_std plus the absolute calibration error; the absolute calibration percentages are listed in Table 3.

$$\begin{eqnarray}&&{F}_{\mathrm{unc}}=\sqrt{{({F}_{\mathrm{source}}^{* }{\mathrm{Abs}}_{\mathrm{cal}})}^{2}+{{\mathrm{Bg}}_{\mathrm{std}}}^{2}}.\end{eqnarray} \tag{ 2 }$$

Table 3. Fractional Calibration Error

Telescope	Instrument	Abs. Calibration
SOFIA	HAWC+	15%
Herschel	PACS	10%
Herschel	SPIRE	15%
WISE	W1-W4	6%
AKARI	N:60,160,WIDE:S,L	10%
Spitzer	IRS	5%
Spitzer	IRAC	5%
Spitzer	MIPS	5%

Download table as:  ASCIITypeset image

The SOFIA/HAWC+ and photometry measure and uncertainties are shown in Table 4 (the full table is available in machine-readable form in the digital version of this paper).

Table 4. UV–FIR Photometric Data ^a

Galaxy	HAWC+A	MIPS	HAWC+C	PACS	HAWC+D	HAWC+E	SPIRE
Name	55 μm	60 μm	89 μm	100 μm	155 μm	214 μm	250 μm
NGC 2537	1302.4 ± 255.0	4040.0 ± 286.4	4893.3 ± 914.1	⋯	6301.9 ± 1085.0	⋯	⋯
Mrk 140	⋯	⋯	⋯	563.2 ± 44.8	⋯	⋯	172.0 ± 26.3
Haro 02	⋯	⋯	⋯	5051.3 ± 254.6	⋯	⋯	911.9 ± 137.2
Mrk 1450	274.8 ± 41.6	⋯	317.0 ± 47.9	272.0 ± 21.8	168.2 ± 27.4	52.5 ± 17.3	48.7 ± 8.4
UM 461	⋯	⋯	⋯	132.0 ± 12.7	⋯	⋯	23.1 ± 7.1
Mrk 1307	724.7 ± 119.9	1000.0 ± 71.6	1277.6 ± 198.6	⋯	616.5 ± 115.4	⋯	⋯
Haro 06	646.3 ± 97.0	⋯	849.4 ± 127.4	742.3 ± 39.2	754.3 ± 115.9	309.8 ± 48.7	197.3 ± 30.0
NGC 4670	1428.2 ± 246.1	3470.0 ± 245.7	3505.9 ± 530.7	⋯	3531.9 ± 713.4	2486.7 ± 428.9	⋯
Mrk 475	⋯	110.0 ± 16.0	⋯	⋯	⋯	⋯	⋯
UGCA 410	⋯	⋯	⋯	⋯	⋯	⋯	⋯
Mrk 930	⋯	⋯	⋯	1140.0 ± 91.9	⋯	⋯	210.4 ± 31.9

Note.

^a This is an abbreviated version of the complete table of flux values used for the analysis. The complete table is available in machine-readable form in the digital version of this paper. The full table contains 49 columns, indicating the telescope and instrument/band for each measurement and the uncertainty. Fluxes are given in mJy.

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

Download table as:  DataTypeset image

The Brown et al. (2014) templates contain SEDs covering a wavelength range from UV to the MIR for 129 nearby galaxies. As described in Section 2, the SEDs of these templates have been normalized and verified with matched-aperture photometry. They provide accurate photometric data, and since our LAHz sample is drawn from these templates, we use the photometric data as presented in Brown et al. (2014).

Some galaxies in our sample were observed with the Herschel Space Observatory, using PACS at 70, 100, and 160 μm, with a field of view (FOV) of $1\buildrel{\,\prime}\over{.} 75\times 3\buildrel{\,\prime}\over{.} 5$, and SPIRE at 250, 350, 500 μm, with an FOV of $4^{\prime} \times 8^{\prime}$.

The Herschel Dwarf Galaxy survey (DGS, Madden et al. 2013; Rémy-Ruyer et al. 2013, 2015) provides photometric and spectroscopic FIR data on 50 dwarf galaxies. This survey includes data for some of the galaxies in our survey: Mrk 930, Mrk 1450, and UM 461, which were observed with all bands of both PACS and SPIRE instruments. A limited set of Herschel photometric data is available for UGCA 410 (also knowns as Mrk 140) and Haro 02 (Herschel Science 2013). See the full Table 4 for a complete listing of the FIR photometric data for our sample. For consistency in our flux estimates, we reduced the available Herschel PACS and SPIRE data using the Herschel pipeline for total intensity photometric measurements. We used the Herschel Interactive Processing Environment (HIPE; Ott 2010) with version 7 of the photometric calibration.

3.3.1. Herschel Photometry

The simplest characterization of dust emission is to fit it with a blackbody F_ν ∝ B_ν (T). In a more realistic situation, this is modified by taking the transfer of radiation into account. This can be done by using a so-called modified blackbody curve, ${F}_{\nu }\propto (1-\exp (-\tau )){B}_{\nu }(T)$. We will use the expression

$$\begin{eqnarray}&&{mBB}(T)={\rm{\Omega }}{B}_{\nu }(T)\left(1-\exp \left[-{\left(\displaystyle \frac{{\lambda }_{0}}{\lambda }\right)}^{\beta }\right]\right),\end{eqnarray} \tag{ 3 }$$

where B_ν (T) is the Planck function and Ω is the solid angle extended by the emission. The opacity is parameterized as $\tau ={({\lambda }_{0}/\lambda )}^{\beta }$, where λ₀ is the wavelength where the optical depth equals unity (e.g., Draine 2006), and the coefficient β represents the opacity spectral index (Kovács et al. 2010). Observationally, β has been derived to be between 1 and 2 (see Hildebrand 1983; Casey 2012). Data from the Planck mission provide a robust derivation of β in the Milky Way (MW), leading to β = 1.8 ± 0.1 (Planck Collaboration et al. 2011). For extragalactic objects it is common to use β in the range 1.3–2.0. Here we will adopt the commonly used value β = 1.5 and keep it constant in our fits. We did try fitting with different β values, but we did not find that it changed our results in any significant way. We also kept the wavelength where the optical depth is unity constant at λ₀ = 100 μm. The two parameters β and λ₀ affect the flux on the long-wavelength Rayleigh–Jeans part of the dust SED relative to the short-wavelength part; for a fixed dust temperature, an increasing λ₀ leads to a relative increase in the Rayleigh–Jeans part, while an increase in β leads to a decrease in the Rayleigh–Jeans part relative to the Wien part of the SED.

When fitting a modified blackbody curve to flux measurements at λ ≳ 60 μm, most galaxies exhibit an MIR excess at wavelengths shorter than 60 μm (e.g., Casey 2012). This excess can be attributed to the presence of a warmer dust component than the one characterizing the longer wavelengths, or an additional heating source, such as an AGN (see Scoville & Kwan 1976). The observed MIR emission can be modeled as a power law, with a cutoff wavelength λ_c (see Casey 2012), or as a second dust component characterized by a warmer dust temperature. Since we do not expect significant AGN activity in the dwarf galaxies making up our sample, we will adopt the latter case and use two modified blackbody curves, with the same β and λ₀, but with different dust temperatures. Our dust SED to be fitted to the observed data is then

$$\begin{eqnarray}\begin{array}{rcl}{mBB} & = & {{\rm{\Omega }}}_{1}{B}_{\nu }({T}_{1})\left(1-\exp \left[-{\left(\displaystyle \frac{{\lambda }_{0}}{\lambda }\right)}^{\beta }\right]\right)\\ & & +\,{{\rm{\Omega }}}_{2}{B}_{\nu }({T}_{2})\left(1-\exp \left[-{\left(\displaystyle \frac{{\lambda }_{0}}{\lambda }\right)}^{\beta }\right]\right).\end{array}\end{eqnarray} \tag{ 4 }$$

With the opacity spectral index kept fixed at β = 1.5 and the critical wavelength at λ₀ = 100 μm, the modified blackbody curve contains four free parameters: two dust temperatures T₁ and T₂, and two normalization constants containing the solid angle of the source Ω. We use the Levenberg–Marquardt fitting algorithm. We explored the fitting using values for the dust opacity coefficient β ranging from 1.2 to 2.0 and for λ₀ in the range 50–200 μm. We found no significant variation in the derived dust temperatures, and we conclude that a fixed β = 1.5 with λ₀ = 100 μm is sufficient and that the models will be characterized by the temperatures and normalization coefficients.

We employ Bayesian statistics inference in multiparameter spaces to determine the most probable set of model parameters ϕ of a model mBB_i in comparison to the data D_i . Bayesian inference is based on Bayes's theorem, which establishes the conditional probability of the parameters ϕ given the data D_i as

$$\begin{eqnarray}&&p(\phi | {D}_{i})=\displaystyle \frac{p({D}_{i}| \phi )* p(\phi )}{p({D}_{i})},\end{eqnarray} \tag{ 5 }$$

where the likelihood of the model is p(D_i ∣ϕ). Then, p(ϕ∣D_i ) is the posterior probability, p(ϕ) are the priors or initial knowledge of the parameters. p(D_i ) is used to normalize the posterior probability. The likelihood is assumed to be

$$\begin{eqnarray*}&&p\propto \exp \left(-{\chi }^{2}/2\right).\end{eqnarray*}$$

By maximizing the posterior probability, we will determine each ϕ parameter's best values for the modified blackbody models mBB_i . Since we are assuming a fixed β = 1.5 and λ₀ = 100 μm for both components, the ϕ parameters are the temperatures T₁ and T₂ and normalization constants Ω₁ and Ω₂. The priors are defined by allowing the ϕ parameters to vary through a very wide range initially, to get a posterior probability, and then Bayesian analysis allows us to update the knowledge of the probability distribution of the parameters and apply new priors multiplying by flat priors that reduce the range of the parameter space. This allows us to reduce the region covered by the parameter and obtain a better resolution of the parameter grid. Inspection of the marginalized probability distribution of each parameter space allows us to determine the appropriate range in which the best set of solutions are located, as well as the most likely value. We can also define the confidence interval for each parameter.

4.1.1. Dust Masses

The dust mass can be derived from the flux density of optically thin emission:

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\displaystyle \frac{{S}_{\nu }\,{D}^{2}}{{\kappa }_{\nu }\,{B}_{\nu }({T}_{d})},\end{eqnarray} \tag{ 6 }$$

where S_ν is the observed flux density, B_ν (T_d ) is the blackbody emission at a given temperature, D is the distance, and κ_ν is the dust absorption coefficient for dust grains at frequency ν. The dust absorption coefficient is a parameter encompassing both the dust grain composition and the grain size distribution and is difficult to derive observationally (e.g., Gordon et al. 2014).

Different values of the ratio of visual extinction to reddening R_V are used to define κ_ν . R_V = 3.1 is the standard value for the diffuse galactic ISM, also used for the Large Magellanic Cloud by Fitzpatrick (1986). Prevot et al. (1984) and Bouchet et al. (1985) found R_V = 2.72 ± 0.21 for the Small Magellanic Cloud (SMC), and Calzetti et al. (2000) found R_V = 4.05 ± 0.8 for starburst galaxies. More recently, Weingartner & Draine (2001) derived κ_ν for different values of R_V for MW dust and LMC and SMC extinction laws. For the dust masses presented in this paper, we use the updated MW extinction laws ⁵ κ_ν presented in Weingartner & Draine (2001), and we calculate κ_ν for the average of R_V = 3.1, 4.0, and 5.5 at λ = 250 μm. This yields κ₂₅₀ = 0.41 m² kg⁻¹. Using the MW, we may be underestimating the dust mass; using the SMC and LMC dust extinction laws produces κ₂₅₀ = 0.52 m² kg⁻¹ and κ₂₅₀ =0.50 m² kg⁻¹, which would yield dust masses ∼20% larger. We do not have observed flux values at λ = 250 μm for all of our galaxies. Instead, we use the flux density provided by our mBB fits to the overall dust SED. Hence, the dust mass, using the flux density at 250 μm, can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{dust}}}{{M}_{\odot }}=188.4\left(\displaystyle \frac{{S}_{\nu }}{\mathrm{Jy}}\right){\left(\displaystyle \frac{D}{\mathrm{Mpc}}\right)}^{2}\,\displaystyle \frac{\left(\exp \left(57.55/T\right)-1\right)}{{\kappa }_{250}}.\end{eqnarray} \tag{ 7 }$$

The dust masses derived in this way are presented iin Section 5.1 and discussed in Section 5.3.

We performed SED fitting using Lightning (Eufrasio et al. 2017, ASCL:1711.009 ⁶ ), which fits nonparametric SFHs in discrete, fixed-width, stellar age bins. We made use of the most recent update to Lightning, which uses an adaptive Markov Chain Monte Carlo (MCMC) algorithm (Doore et al. 2021). For this work we have used SFH bins of 0−10 Myr, 10−100 Myr, 100 Myr−1 Gyr, 1−5 Gyr, and 5−13.3 Gyr and assumed a Kroupa initial mass function (Kroupa 2001) and a metallicity of 0.2 Z_⊙. We used PEGASE (Fioc & Rocca-Volmerange 1997, 1999) stellar populations. For intrinsic attenuation, we adopted a Calzetti et al. (2000) extinction law, modified as in Noll et al. (2009) to include a UV bump at 2175 Å and a parameter δ to control the slope of the extinction curve. For instance, δ = 0 produces the Calzetti curve, a lower δ corresponds to a higher slope of the extinction curve, and a higher δ corresponds to a flatter extinction curve. Draine & Li (2007) MW dust emission was used. For a more thorough description of the prescriptions adopted by Lightning, we refer the reader to Eufrasio et al. (2017).

The IR photometry previously described described was combined with literature SEDs from Brown et al. (2014), producing an SED with 45 data points from the FUV to the FIR (including GALEX, Swift UVOT, SDSS, 2MASS, WISE, Spitzer, SOFIA, IRAS, AKARI, and Herschel). Our fits explore a 13-dimensional parameter space, with 5 parameters regarding the SFH intensities ψ_i , 3 parameters $({\tau }_{V}^{\mathrm{DIFF}}$, δ, and ${\tau }_{V}^{\mathrm{BC}}$) controlling the dust extinction, and 5 parameters controlling dust emission (α, U_min, U_max, γ, and q_PAH).

The range allowed for each parameter is shown in Table 5. This range of δ allows the dust extinction curve to be steeper or flatter relative to the Calzetti extinction curve (for more details, see Eufrasio et al. 2017). Also, the five dust emission parameters explore the full space of the Draine & Li (2007) models. An illustration of how these parameters affect the SED is shown in Leja et al. (2017). We run the adaptive MCMC algorithm for 3 × 10⁵ steps and derive statistics from the last 5000 steps of the chain. All the derived quantities quoted here (SFHs, SFRs, M_⋆, M_dust) are derived from these last 5000 steps.

The blackbody modeling described in Section 4.1 shows that a reasonable description of the dust SED can be achieved by using a two-component modified blackbody fit. The warm component is less well constrained than the cold component owing to a smaller number of photometric data points at wavelengths λ ≲ 60 μm. Nevertheless, the model SED joins smoothly with the NIR and MIR fluxes of the Brown et al. (2014) templates.

In Figure 5 we show multiparameter-space diagnostic plots for two galaxies in our sample: Haro 06 and Mrk 1450. We can see the marginalized distribution for each ϕ parameter in the upper frame of each column; the lower frames show the marginalized probability in blue and the correlation values between the two parameters. The value that maximizes the χ² distribution (likelihood) for each parameter is indicated with a red star; this is also shown as a red dashed line in the upper frames for each parameter. The confidence interval is shown in light red and indicates the 16th and 84th percentiles in the probability distribution, which imply that 68% of the probability is inside that region.

Zoom In Zoom Out Reset image size

This visualization of the probability distributions allows us to understand how the different parameters relate. For instance, we can see how the temperature T₁ and the constant Ω₁ are anticorrelated, while the two constants Ω₁ and Ω₂ are positively correlated for both galaxies. In the upper right box we see the best-fit model indicating the two blackbody components representing a warm and a cold dust component. The corresponding Brown et al. (2014) template is shown as a reference, but it is not part of the fit. The resulting best-fit model agrees smoothly with the Brown templates.

In Figure 6 we present the SEDs produced by the two-component blackbody models for all of our LAHz galaxies. In a few cases the warm and the cold components form two humps or peaks in the SED, often found for star-forming galaxies. However, in several cases, the warm and cold components blend together and form a broad single-peaked SED. In the figures we also include the corresponding NIR and MIR spectra from Brown et al. (2014). These spectra were not part of the fit, but they merge smoothly with the dust SED. In three cases (NGC 2537, Haro 02, and NGC 4670), the MIR spectra suggest that there might be additional flux at wavelengths between 24 and 70 μm, compared to what is indicated in the two-component blackbody fits.

Zoom In Zoom Out Reset image size

5.1.1. Dust Temperatures

We obtained the dust temperatures characterizing the warm and cold components from the modified blackbody fits. For comparison, we use the DGS presented by Rémy-Ruyer et al. (2013), which includes local dwarf galaxies, and in order to provide a comparison with more evolved galaxies, we use a subset of 58 galaxies from the Key Insights on Nearby Galaxies: a Far-Infrared Survey with Herschel (KINGFISH; Kennicutt et al. 2011), where we exclude the irregular dwarf galaxies and early-type galaxies from the original set of 109 galaxies.

In order to compare our results with the galaxies in the DGS and KINGFISH samples, we reanalyzed the FIR photometric data from these surveys using the same modified blackbody method as for our LAHz galaxies. For each galaxy we derive the temperatures characterizing the warm and cold components, as well as the temperature corresponding to the peak of the combined SED, defined as

$$\begin{eqnarray}&&{T}_{\mathrm{peak}}=\displaystyle \frac{b}{{\lambda }_{\mathrm{peak}}},\end{eqnarray} \tag{ 8 }$$

where b is Wien's displacement constant. All the galaxies were characterized using two dust components. The average dust temperatures for the sample of LAHz are T_d,warm ∼ 124 K and T_d,cold ∼ 42 K. This is significantly higher than for star-forming disk galaxies in the KINGFISH sample, with T_d,warm ∼ 68 K and T_d,cold ∼ 23, and higher than for the DGS sample, with T_d,warm ∼ 74 K and T_d,cold ∼ 24. The range of temperatures for the DGS survey for the cold component is, however, quite large, ranging from 18 to 41 K (Rémy-Ruyer et al. 2013). The results are shown in Figure 9.

As described in Section 4.2, the Lightning models fit the whole UV–FIR SED, while self-consistently applying an extinction correction. The FIR is modeled using Draine & Li (2007) templates. The resulting SEDs for our LAHz galaxies are shown in Figure 7. The black data points represent the observed SED; the blue line is the spectra without dust attenuation; the red line is the SED with attenuation; the orange line is the best-fit dust emission. In Table 7 we show luminosities, dust masses, SFR averaged over the last 100 Myr, and sSFRs derived from the Lightning models.

Zoom In Zoom Out Reset image size

We derive dust masses and L_IR using the SEDs produced with the two different modeling methods described in Section 4. The dust masses and L_IR from the mBB models and the Lightning method are presented in Tables 6 and 7, respectively. The L_IR is derived by integrating the SED from 8 to 1000 μm. The LAHz galaxy luminosities cover a wide range of $\mathrm{log}({L}_{\mathrm{IR}}/{L}_{\odot })=[7.4,10.2]$ for the blackbody models; similarly, the Lightning method luminosities range between $\mathrm{log}({L}_{\mathrm{IR}}/{L}_{\odot })=[7.5,10.3]$.

The derived dust masses using the blackbody models are in the range of $\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\odot })=[3.1,6.3]$. As expected, the cold dust component dominates the total dust mass. We also derive the dust mass from the Draine & Li best fits in the Lightning models: $\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\odot })=[3.2,6.3]$. Some of the galaxies in our sample have published dust mass values, derived using a subset of the photometric data used here (e.g., Engelbracht et al. 2008; Rémy-Ruyer et al. 2013), which agree very well with the results presented here. For the discussion presented in the rest of this paper, we use the dust masses derived with the mBB method, since this allows the comparison with dust masses derived for survey samples like KINGFISH and DGS using the same methodology.

We compare the dust masses and FIR luminosities for our LAHz galaxies, derived with the modified blackbody and Lightning methods. The latter method uses the Draine & Li (2007) templates to model the FIR SED. The dust masses derived from the mBB and the Draine & Li models show dispersion of ∼0.3 dex but show no systematic difference. The FIR luminosities, however, do show a systematic offset, where the Draine & Li (2007) based models give ∼20% higher L_FIR than those derived by integrating over the modified blackbody models. This is expected since the Draine & Li (2007) models include PAH emission in the NIR and better account for the NIR to MIR emission than the blackbody models. In the rest of the paper we will use the mBB-derived dust masses for comparison with the DGS and KINGFISH surveys, as well as the Lightning FIR luminosities for discussions and conclusions.

In Figure 8 we show the SFH for the LAHz galaxies in our sample. For comparison, we also show the SFH of M51 from Eufrasio et al. (2017). The SFH was derived using the Lightning fitting package, and it is presented in five look-back time bins: 0−10 Myr, 10–100 Myr, 0.1–1 Gyr, 1–5 Gyr, and 5–13.3 Gyr. The Lightning modeling treats the SFR as constant over each time bin. The solid red lines denote the median SFR for each time bin from the simulations, and the yellow region denotes the 16th and 84th percentiles, enclosing the 68% confidence intervals. The SFHs of the LAHz are for the most part very different from those of more massive galaxies. For instance, the SFH of M51, a large spiral galaxy gravitationally interacting with a neighbor galaxy, shows a more complex SFH, with a significant star formation contribution in bins corresponding to ages older than 5 Gyr.

Zoom In Zoom Out Reset image size

The SFH temporal resolution is nonlinear, with longer timescales for higher look-back times. Even a very low SFR in the older bins can significantly contribute to the total stellar mass. However, for six galaxies we recover SFH solutions compatible with null stellar masses older than 1 Gyr: Mrk 1450, UM 461, Mrk 1307, Mrk 475, UGCA 410, and Mrk 930. These galaxies do not present statistically significant evidence of a stellar population older than 1 Gyr. This suggests that they are the youngest galaxies in our sample. Two of the galaxies, NGC 2537 and NGC 4670, show star formation activity in all five time bins, showing a statistically significant presence of an old stellar population. Hence, despite being designated as BCDGs, the stellar masses of these two galaxies are dominated by old stars. The situation is less clear for three of the galaxies in our sample, Mrk 140, Haro 02, and Haro 06, where the Lightning results indicate that it is probable that they had a low level of star formation activity at times ≳1 Gyr (see Figure 8).

We consider the six LAHz galaxies that have SFHs consistent with being younger than ∼1 Gyr as our prime candidates for being local analogs. NGC 2537 and NGC 4670 are clearly not young galaxies. Mrk 140, Haro 02, and Haro 06 are likely to have an old stellar population, although Haro 06 is a borderline candidate for being a local analog.

When fitting a two-component modified blackbody curve to our LAHz galaxies, the DGS, and spiral galaxies in the KINGFISH sample, we find that the DGS and LAHz galaxies show a higher T_peak than the more massive KINGFISH galaxies (Figure 9). While the average T_peak is 33 and 35 K for the LAHz and DGS samples, respectively, the KINGFISH galaxies have an average of T_peak = 23 K. When we compare the warm/cold dust component separately, we see that the LAHz galaxies are characterized by higher dust temperatures than both the KINGFISH and DGS galaxies. These results are consistent with a general trend of warmer dust in metal-poor galaxies compared with the more metal-rich and more massive counterparts (e.g., Rémy-Ruyer et al. 2013), as well as for high-z galaxies (e.g., Faisst et al. 2017; Sommovigo et al. 2020).

Zoom In Zoom Out Reset image size

In Figure 10 (top left panel), we show the dust-to-stellar mass ratio as a function of stellar mass for our sample, together with those of the DGS and KINGFISH spiral galaxies. While the KINGFISH galaxies seem to have a weak dependence on stellar mass, no such correlation is apparent for our LAHz galaxies or for galaxies in the DGS sample. In fact, the dust-to-stellar mass ratio for the LAHz galaxies seems to be approximately constant, with a median ratio $\mathrm{log}({M}_{\mathrm{dust}}$/M_⋆) = (8.5 ± 5) × 10⁻⁴. A similar argument can be made for the DGS sample, but with a slightly lower average dust-to-stellar mass ratio. The dashed line shown in Figure 10 shows the anticorrelation between dust-to-stellar mass ratio and stellar mass found by Clemens et al. (2013). This anticorrelation was obtained for galaxies from the Planck Survey, and the fit is dominated by galaxies in the stellar mass range $\mathrm{log}{M}_{* }/{M}_{\odot }\sim 10-11$. The anticorrelation is surprising in view of the fact that the gas-phase metallicity scales with stellar mass (e.g., Tremonti et al. 2004) and one would expect a lower M_dust/M_⋆ ratio for low-mass galaxies. The Clemens et al. (2013) relation has been interpreted as being due to a rapidly increasing gas fraction for lower-mass galaxies, more than offsetting the metallicity effect. However, the results given here, for both our sample and that of DGS, show that this relation cannot be extended to low-mass galaxies.

Zoom In Zoom Out Reset image size

The mass–metallicity relation is established down to $\mathrm{log}{M}_{* }/{M}_{\odot }\sim 8.5$. Below this stellar mass, depletion of metals through galactic winds starts to become important (e.g., Tremonti et al. 2004). Our LAHz galaxies have $\mathrm{log}{M}_{* }/{M}_{\odot }$ in the range 7.0–9.4, straddling the $\mathrm{log}{M}_{* }/{M}_{\odot }=8.5$ value. In Figure 10, top right panel, we show the dust-to-stellar mass ratio as a function of metallicity. In the figure we also include the DGS and KINGFISH spiral samples. Based on the latter two samples, Rémy-Ruyer et al. (2013) defined a correlation between the dust-to-stellar mass ratio and gas-phase metallicity, shown as a dashed line in Figure 10. We find no obvious correlation between the dust-to-stellar mass ratio and metallicity for our LAHz sample, and the overall dispersion of dust-to-stellar mass ratios is large. As discussed in Rémy-Ruyer et al. (2013), this could at least in part be due to observational uncertainties. However, for the metal-poor and low-mass galaxies, this could also be due to the impact of individual SFHs, where the interstellar gas can be promptly enriched through recent star formation or depleted owing to galactic winds. Evidence for both effects has been seen in these types of galaxies; Lagos et al. (2018) found spatial variations in the gas-phase metallicity in UM 461, attributed to recent infall of low-metallicity gas, and Haro 02 shows clear evidence of a galactic gas outflow (e.g., Meier et al. 2001; Beck et al. 2020).

For comparison, we also include two of the lowest-metallicity galaxies known, SBS 0335−052E and I Zw 18, in Figure 10. Both of these galaxies have very low gas-phase metallicity, 12 + logO/H ∼ 7.2, but even at these extremely low metallicities they still show dust continuum emission (SBS 0335−052E, Hunt et al. 2014; I Zw 18, Wu et al. 2007). These two particularly low-metallicity galaxies were studied by Hunt et al. (2014) and found to have dust-to-stellar mass ratios within the range of values found for the galaxies presented in Figure 10. The dust-to-stellar mass ratio of I Zw 18 is $\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\star })=-3.60$, which is consistent with the LAHz average of $\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\star })=-3.17$. SBS 0335−052E, however, has $\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\odot })\approx -2.15$, which is significantly higher than the LAHz and DGS galaxies, comparable to metal-rich spiral galaxies. While this wide disparity among the most metal-poor galaxies could be due to observational uncertainties, it could also be an effect of the above-mentioned prompt enrichment and galactic outflow/inflow processes affecting low-mass galaxies.

The bottom right panel of Figure 10 shows the L_IR/M_dust ratio as a function of metallicity for the LAHz galaxies (red squares). For comparison, we also show the DGS and KINGFISH galaxies. The distribution of L_IR/M_dust is indicated on the side for each group. The L_IR/M_dust of the LAHz galaxies is distributed in the same range as the DGS (Rémy-Ruyer et al. 2013). This ratio tends to increase for lower metallicities, and when we examine the left region in the range of metallicities, the LAHz galaxies exhibit a higher L_IR/M_dust ratio than the DGS or KINGFISH galaxies for the corresponding metallicity. It is not possible with the present data to determine whether this is due to different grain properties or is simply an effect of the higher dust temperatures characterizing the dust emission in these galaxies. As shown in Figure 10, this trend extends to the extremely metal-poor galaxies SBS 0335−052E and I Zw 18, covering more than an order of magnitude in both metallicity and L_IR/M_dust. These two galaxies are dominated by warm dust; their single-component dust SEDs are characterized by T_dust = 59 and 67 K, respectively (Hunt et al. 2014). While these temperatures are high, they are within the range of the temperatures characterizing our LAHz galaxies.

The high temperatures characterizing the warm dust component of the LAHz galaxies, as well as the high MIR fluxes relative to the peak FIR flux of several of the LAHz galaxies, could potentially be due to the presence of a nonstellar heating source, such as AGNs. AGNs in low-mass galaxies are more difficult to detect than their counterparts in massive galaxies. None of our LAHz galaxies are, however, believed to harbor an active AGN. Reines et al. (2020) used the Very Large Array to show that previously identified radio sources that could potentially be AGNs in Haro 02 and Mrk 1307 are most likely associated with star formation activity.

We explore the possibility of AGN, LINER, or Seyfert activity in the LAHz galaxies by plotting them in a BPT diagram (Baldwin et al. 1981). This diagnostic uses the Hα/[Nii] ratio together with the [O iii] λ5007/Hβ ratio to distinguish between star formation and AGN activity. The spectral line observations presented in Moustakas & Kennicutt (2006) include all of our LAHz galaxies, and the BPT diagram of the galaxies in Moustakas & Kennicutt (2006) is shown in Figure 11, with our LAHz galaxies marked as red circles. All LAHz galaxies, except Mrk 140, fall in the star formation part of the diagram. Mrk 140 falls close to the line depicting maximum star formation activity as defined by Kewley et al. (2001). Since there are no other indications suggesting that Mrk 140 harbors an AGN, we will treat it as a pure star-forming system.

Zoom In Zoom Out Reset image size

In Figure 12 we show the correlations between M_dust/SFR and M_dust/M_⋆ as a function of the sSFR (sSFR = SFR/M_⋆). As shown in da Cunha et al. (2010), the dust-to-stellar mass ratio strongly correlates with the sSFR, confirming that dust mass and star formation are tightly related. In Figure 12 we color-code the stellar mass of each galaxy. More massive galaxies exhibit a lower sSFR and higher M_dust/SFR ratio, contrary to the lower-mass galaxies, which show higher sSFRs, tracing a clear decreasing relation between $\mathrm{log}({M}_{\mathrm{dust}}/\mathrm{SFR})$ and sSFR. In the right panel we plot M_dust/M_⋆ ratio versus sSFR. Here we see an increasing trend, where M_dust/M_⋆ ratio increases as sSFR increases. Interestingly, the less massive galaxies seem to follow a parallel path to the more massive ones. In both cases, the LAHz candidates fall in the right side, showing that they have higher sSFR than more massive and dusty galaxies.

Zoom In Zoom Out Reset image size