Extinction Maps and Dust-to-gas Ratios in Nearby Galaxies with LEGUS

Dust has large and varied impacts on studies of star formation, chemical evolution, and galaxy evolution. The most obvious of these impacts is the absorption and scattering of optical and ultraviolet (UV) light by dust. Typically, up to 50% of the total stellar energy in a galaxy is attenuated by dust, hampering the interpretation of galaxy spectral energy distributions (SEDs) for their fundamental parameters, such as age, stellar population mix, star formation rates, and stellar initial mass functions (Calzetti 2001). This, in turn, impacts our ability to obtain star formation histories and constrain theories of galaxy evolution and star formation.

Dust also plays a major role in star formation, both as a means of radiative feedback and as a catalyst for the formation of molecular hydrogen and other molecules (Mathis 1990; Draine 2003; Caxaux & Spaans 2009; McKinnon et al. 2016). Finally, dust and the dust-to-gas ratio have effects on galaxy and chemical evolution models, since dust acts as a sink for many metals (McKinnon et al. 2016; Aoyama et al. 2017).

Extinction, as measured toward individual stars, includes the effects of absorption and scattering. For an assumed dust population, the extinction is linearly proportional to the column density of that population through the proportionality of optical depth to dust column density. For the Milky Way (MW), Bohlin et al. (1978) and Rachford et al. (2009) gave the relation

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{{\rm{H}}}}{E(B-V)}=5.8\times {10}^{21}\,{\rm{H}}\,{\mathrm{cm}}^{-2}\,{\mathrm{mag}}^{-1},\end{eqnarray} \tag{ 1 }$$

where ${N}_{{\rm{H}}}={N}_{{\rm{H}}{\rm{I}}}+2{N}_{{{\rm{H}}}_{2}}$ is the total neutral hydrogen column density and E(B − V) is the reddening in the B − V color. Here N_{H i} and ${N}_{{{\rm{H}}}_{2}}$ are the column densities of atomic and molecular hydrogen, respectively.

For a Galactic extinction law $\left({R}_{V}=\tfrac{{{\boldsymbol{A}}}_{V}}{{\boldsymbol{E}}({\boldsymbol{B}}-{\boldsymbol{V}})}=3.1\right)$, Equation (1) gives

$$\begin{eqnarray}\displaystyle \frac{{A}_{V}}{{N}_{{\rm{H}}}} & = & \displaystyle \frac{3.1}{5.8\times {10}^{21}\,{\rm{H}}\,{\mathrm{cm}}^{-2}\ {\mathrm{mag}}^{-1}}\\ & = & 5.3\times {10}^{-22}\,\mathrm{mag}\,{\mathrm{cm}}^{2}\,{{\rm{H}}}^{-1},\end{eqnarray} \tag{ 2 }$$

where A_V is the extinction in the V-band magnitude. Other galaxies with similar metallicities are expected to have a similar relation.

Measurements of the dust-to-gas mass ratio are often done through direct measurements of the mass of neutral gas and dust from radio and IR observations, such as in the analysis of Draine et al. (2007), who gave a dust-to-gas mass ratio that varies with oxygen abundance,

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{dust}}}{{M}_{{\rm{H}}}}=0.010\displaystyle \frac{({\rm{O}}/{\rm{H}})}{{({\rm{O}}/{\rm{H}})}_{\odot }},\end{eqnarray} \tag{ 3 }$$

where $\tfrac{{M}_{\mathrm{dust}}}{{M}_{{\rm{H}}}}$ is the ratio of dust mass to neutral gas mass and $\tfrac{({\rm{O}}/{\rm{H}})}{{({\rm{O}}/{\rm{H}})}_{\odot }}$ is the ratio of the metallicity to solar metallicity. More recent work suggests the relationship in the MW might be closer to $\tfrac{{M}_{\mathrm{dust}}}{{M}_{{\rm{H}}}}=0.0091\tfrac{({\rm{O}}/{\rm{H}})}{{({\rm{O}}/{\rm{H}})}_{\odot }}$ (Sofia & Parvathi 2009; Draine 2011).

Equation (3) is consistent with the observed values of the dust-to-gas mass ratio to within a factor of 2 for each galaxy in the Draine et al. (2007) sample. New calibrations by Faniullo et al. (2015) and the Planck Collaboration et al. (2016) indeed suggest that the Draine & Li (2007) models used in Draine et al. (2007) overpredict dust masses by about a factor of 2 regardless of metallicity. However, Draine et al. (2007) had no galaxies with $12+\mathrm{log}\,{\rm{O}}/{\rm{H}}\lt 7.6$, and Galametz et al. (2011) and Rémy-Ruyer et al. (2014) found that low-metallicity galaxies have a lower-than-expected dust-to-gas mass ratio. Galametz et al. (2011) suggested that perhaps the lack of submillimeter maps on some dwarf galaxies may lead to an underprediction of the dust, while Rémy-Ruyer et al. (2014) suggested that the decreased dust-to-gas mass ratio is real and proposed a broken power law for the relationship between galaxy metallicity and dust-to-gas mass ratio.

All of the above studies require an assumption of the CO-to-H₂ conversion factor X_CO, which also appears to have a dependence on metallicity for non-MW-like environments. Draine et al. (2007) used a constant X_CO = 4 × 10²⁰ cm² K⁻¹ (km/s)⁻¹. More recent studies tend to use the lower value of X_CO = 2 × 10²⁰ cm² K⁻¹ (km/s)⁻¹ from Leroy et al. (2009). Rémy-Ruyer et al. (2014) compared the Leroy et al. (2009) fixed X_CO with one that varies with metallicity using the Schruba et al. (2012) scaling (O/H)⁻² for a sample of 126 galaxies over a 2 dex metallicity range. This results in a slight decrease in the dust-to-gas mass ratio for a low-metallicity galaxy. Brinchmann et al. (2013) found that using the metallicity-dependent $\mathrm{log}{X}_{\mathrm{CO}}=-1.01\mathrm{log}\,{\rm{O}}/{\rm{H}}+29.28$ from Boselli et al. (2002) is required to achieve good agreement between dust-to-gas mass ratios obtained from CO and H i and their X_CO-independent method using optical spectroscopy.

The ability to correct photometry for dust is vital for extragalactic studies of the interstellar medium (ISM) and stellar populations. Extinction maps can be used to correct optical (particularly Hα) and UV images for the effects of dust in order to better determine star formation rates and histories in nearby galaxies. This is particularly important for studies of the initial mass function (e.g., Ashworth 2017) and mass-to-light ratios, as classifying young massive stars often relies heavily on the UV, where extinction is severe.

The "pair method" is a common way of obtaining reddening corrections for individual stars by comparing a reddened spectrum or observed multiband photometry to an unreddened spectrum from a star of the same or similar spectral type (Elmegreen 1980; Cardelli et al. 1992). However, this is difficult for crowded regions, as additional stars will contaminate the spectrum of an individual star (Maiz-Apellaniz et al. 2014). Previous studies, such as the Hubble Tarantula Treasury Project (Sabbi et al. 2013; De Marchi et al. 2016) and the Panchromatic Hubble Andromeda Treasury (PHAT) survey (Dalcanton et al. 2015), have used near-IR photometry and isochrone matching of red clump stars to measure extinctions to circumvent this problem. This is difficult for stars in more distant galaxies, as the lower end of the red giant branch is too faint to obtain accurate photometry without extremely deep imaging, and crowding of old stars becomes an issue. Zaritsky et al. (2002) and Zaritsky (1999) fit four-band optical photometry of stars in the Magellanic Clouds to stellar atmosphere models to obtain extinction corrections and generate extinction maps. These studies relied on ground-based photometry, so the SED-fitting method they pioneered has not seen much use outside the Local Group until recently. Gordon et al. (2016) used SED fitting of the PHAT data to derive extinction corrections for stars in M31 and very recently made their fitting code publicly available.

Reddening-free Q parameters have been used to identify the spectral type of massive stars and obtain reddening corrections for decades. Johnson & Morgan (1953) generated a table of $Q=(U-B)-\tfrac{E(U-B)}{E(B-V)}(B-V)$ for B- and O-type stars from photometry of nearly 300 stars within the MW. Van den Bergh & Hagen (1968) used a similar system to derive reddenings for clusters in the Magellanic Clouds. This early system worked well for B-type stars; however, the relation between Q-value and spectral type flattened for O-type stars, making it difficult to distinguish colors and temperatures of the most massive stars. Despite the difficulty in determining the spectral types of the most massive stars, reddening corrections are still accurate. More recent work from Kim et al. (2012) and Blair et al. (2015) used isochrone-matching and slightly different definitions of Q to more precisely measure intrinsic colors, improving upon older tables of spectral types and Q-values.

In this work, we use the photometry of young massive stars described in Section 2 to generate extinction maps for five galaxies using the isochrone-matching method described in Kim et al. (2012) and in more detail in Section 3. In Section 4, we discuss the subset of stars used to generate the extinction maps shown in Section 5. We calculate dust-to-gas mass ratios for each galaxy in Section 6 and examine the correlation between dust-to-gas mass ratio and metallicity, using metallicity values and neutral gas maps obtained from the literature, in Section 7.1. Finally, we compare the extinctions obtained from isochrone matching to those obtained via the SED-fitting method used in the Bayesian Extinction And Stellar Tool (BEAST) in Section 7.2.

The Legacy ExtraGalactic Ultraviolet Survey (LEGUS) has observed 50 star-forming galaxies at distances ranging from approximately 3 to 16 Mpc. For each galaxy, LEGUS has generated five-band stellar and cluster photometric catalogs with a combination of the Wide-Field Camera 3 (WFC3) and Advanced Camera for Surveys (ACS) instruments (Adamo et al. 2017; Sabbi et al. 2018). Our wavelength coverage runs from the near-UV to the I band, using the F275W (NUV), F336W (U), F438W (B), F555W (V), and F814W (I) filters for the WFC3 instrument and the F435W, F555W, and F814W filters for the ACS instrument. The survey is described in full in Calzetti et al. (2015). We achieve greater resolution with HST imaging compared to ground-based observations, providing us with the ability to obtain individual stellar photometry for the most nearby galaxies.

In general, the LEGUS observations were designed to reach a depth of ${m}_{{\rm{F}}275{\rm{W}}}=26.0\,\mathrm{mag}$ with a signal-to-noise ratio = 6, with similar depth in the other filters. Individual stellar positions and fluxes were determined using the ACS and WFC3 point-spread-function-fitting (PSF-fitting) modules provided in the DOLPHOT photometric package (Dolphin & Kennicutt 2002). Sources with a robust detection in a band do not deviate more than 3σ from the average χ² value for the PSF fit for the entire galaxy and have a photometric error of less than 0.5 mag (Sabbi et al. 2018). The version of the LEGUS stellar catalogs used in this analysis is described in Sabbi et al. (2018). We use the V2 catalogs, which contain only the stars with the highest-quality photometry, and we limit our sample to stars with robust detections in the NUV, B-, V-, and I-band filters. These catalogs are uncorrected for foreground MW or local extinction. Both of these corrections are discussed in Section 3. Table 1 lists the pointings used in this study and the instrument and filter combinations used for each pointing.

In this work, we focus on five of the 50 galaxies in the LEGUS sample: NGC 628, UGC 4305 and UGC 5139 (Holmberg II and Holmberg I), NGC 6503, and NGC 7793. Table 2 lists the properties of our galaxy sample. The five galaxies together represent a good subset of the LEGUS sample in distance, inclination, metallicity, and size. NGC 628, NGC 6503, and NGC 7793 are all spiral galaxies with a range of masses, though NGC 7793 is categorized as a flocculent. UGC 4305 and UGC 5139 are dwarfs. The three spirals are at a range of distances, with NGC 6503 roughly in the middle of the LEGUS distance range. NGC 628 is one of the most distant galaxies in the sample, and NGC 7793 is relatively nearby, at roughly the same distance as the two dwarfs. NGC 6503 is also the most inclined of the LEGUS galaxies at about 70°, while the other four can be approximated as face-on. The large sizes (∼12' and ∼14') and multiple LEGUS pointings of NGC 628 and NGC 7793 provide us with an opportunity to study radial variations within each galaxy. Both galaxies also have published metallicity gradients from Sánchez-Blázquez et al. (2014) and Stanghellini et al. (2015), providing us with a larger sample of metallicities for this study. Later studies will expand upon this sample.

Table 2.  Properties of the LEGUS Galaxy Sample

Namea	vHa	Morph.a	Tb	Inclin.a	Dist.c	Methodd	12+log(O/H)e		SFR(UV)f	M*g	M(H i)h
	(km s−1)			(deg)	(Mpc)		(PT)	(KK)	(M⊙ yr−1)	(M⊙)	(M⊙)
NGC 0628	657	SAc	5.2(0.5)	25.2	9.9	SNII	8.35	9.02	3.67	1.1E10	1.1E10
NGC 6503	25	SAcd	5.8(0.5)	70.2	5.27	TRGB	8.69		0.32	1.9E09	1.3E09
NGC 7793	230	SAd	7.4(0.6)	47.4	3.44	Ceph	8.31	8.88	0.52	3.2E09	7.8E08
UGC 4305	142	Im	9.9(0.5)	37.1	3.05	Ceph	7.92		0.12	2.3E08	7.3E08
UGC 5139	139	IABm	9.9(0.3)	33.6	3.98	TRGB	8.00		0.02	2.5E07	2.1E08

Notes. Reproduced from Calzetti et al. (2015).

^aGalaxy name, recession velocity, and morphological type as listed in the NASA Extragalactic Database (NED). Inclination in degrees derived from the sizes listed in NED. ^bRC3 morphological T-type as listed in Hyperleda (http://leda.univ-lyon1.fr) and discussed in Kennicutt et al. (2008) for the Local Volume Legacy Survey (LVL) galaxies, from which the LEGUS sample is derived. In that paper, T-type = 11 is adopted for galaxies misclassified as early types while being compact irregular or blue compact galaxies (BCGs). Uncertainties on the morphological classification are in parentheses. Some of the galaxies have large uncertainties, and they may be misclassified. ^cRedshift-independent distance in megaparsecs. ^dMethods employed to determine the distances. In order of decreasing preference: Cepheids (Ceph), tip of the red giant branch stars (TRGB), and supernova Type II plateau (SNII). ^eCharacteristic oxygen abundances of the galaxies. For NGC 628 and NGC 7793, this is the globally averaged abundance (Moustakas et al. 2012). The two columns, PT and KK, are the oxygen abundances on two calibration scales: the PT value, in the left column, is from the empirical calibration of Pilyugin & Thuan (2005); the KK value, in the right column, is from the theoretical calibration of Kobulnicky & Kewley (2004). The dwarves and NGC 6503 have only one metallicity available (via the "direct" method from Croxall et al. 2009 and Tikhonov 2014), so the value straddles the two columns. ^fStar formation rate (M_⊙ yr⁻¹), calculated from the GALEX far-UV, corrected for dust attenuation, as described in Lee et al. (2009). ^gStellar masses (M_⊙), obtained from the extinction-corrected B-band luminosity and color information, using the method described in Bothwell et al. (2009) and based on the mass-to-light ratio models of Bell & de Jong (2001). ^hH i masses, using the line fluxes listed in NED, applying the standard formula: $M({\rm{H}}\,{\rm{I}})[{M}_{\odot }]=2.356\times {10}^{5}\,{D}^{2}S$, where D is the distance in Mpc and S is the integrated 21 cm line flux in units of Jy cm s⁻¹.

Download table as:  ASCIITypeset image

Since the LEGUS stellar photometric catalogs are not corrected for the foreground MW extinction, we first find the foreground MW extinction from the values in Schlafly & Finkbeiner (2011), converted to the HST filters for each galaxy. We then apply those extinctions to the stars in each galaxy's photometric catalog. Next, using the method described in Kim et al. (2012), we match the five-band photometry to theoretical isochrones to obtain extinction corrections for individual stars. First, we calculate reddening-free Q-values for massive stars using Equation (4),

$$\begin{eqnarray}&&{Q}_{{NBVI}}=({\bf{NUV}}-B)-\displaystyle \frac{E({\bf{NUV}}-B)}{E(V-I)}(V-I),\end{eqnarray} \tag{ 4 }$$

where NUV, B, V, and I are the observed magnitudes in the HST filter corresponding to that band and E(NUV − B) and E(V − I) are the selective extinctions in those colors. The selective extinction ratio is obtained using a standard MW extinction law (R_V = 3.1), which corresponds to $\tfrac{E({\bf{NUV}}-B)}{E(V-I)}=1.5999$ (Cardelli et al. 1989). Table 3 shows the extinction coefficients for the five bands in the LEGUS photometry. Changes in the choice of the extinction law effectively change the slope of the reddening vector (i.e., slopes of constant Q lines), which can cause changes in the derived extinction by a factor of 2 or more in the most extreme cases. We test the effects of this shift in Section 6 on our extinction maps, finding that the change in the derived extinction for the full map is negligible, as the extreme stars are averaged out. We calculate the selective extinction ratio and Q based on the center wavelength of each bandpass. There is a slight variation of Q if it is calculated based on coefficients calculated by modeling the shape of the spectrum in each bandpass; however, this effect is negligible compared to the effect of changing the extinction law.

Table 3.  Extinction Coefficients

HST Filter	λcentral (Å)a	$\tfrac{{A}_{\lambda }}{{A}_{V}}$ b
F275W (NUV)	2750	1.9893
F336W (U)	3360	1.6328
F438W/F435W (B)c	4380	1.3319
F555W (V)d	5550	0.9886
F814W (I)d	8140	0.5777

Notes.

^aApproximate central wavelength of the filter in angstroms. ^bHere $\tfrac{{A}_{\lambda }}{{A}_{V}}$ is calculated using the approximate center wavelength of the filter. ^cThe F435W filter is the ACS filter corresponding to the same approximate wavelength range as the WFC3 F438W filter. ^dThe F555W and F814W filters share the same name and wavelength range on both the ACS and WFC3 instruments.

Download table as:  ASCIITypeset image

The reddening-free Q-value is defined so that it only depends on the spectral type of a star and can be defined differently depending on the availability of photometric bands. Binney & Merrifield (1998) used the UBV system in their definition, while Kim et al. (2012) used a similar Q to Equation (4) but substituted the U band for our NUV band, which has a higher signal-to-noise ratio in the data set used in their work than in the NUV band. In the LEGUS U-band photometry, we also find lower photometric errors than we see in the NUV band, with average errors approximately 0.02 less in the U band. However, we use the NUV band because NUV is more sensitive to extinction, allowing us to measure even small amounts of dust. Massive stars also have the peak of their light distribution further in the UV, making it easier for us to differentiate masses. For stars in the range 10–30 M_⊙, the Q-value employed by Kim et al. (2012) varies by approximately 0.2, whereas the Q-value defined in Equation (4) varies by nearly 0.3 over the same mass range.

Once Q is obtained, we then determine the extinction values for each star by matching each individual star's Q-value to the Q-value corresponding to the theoretical isochrones. For stars more massive than 2 M_⊙, the isochrones for the various luminosity classes are very similar; therefore, we do not initially need to know the luminosity class of the star in order to determine the star's intrinsic color. This region is enclosed by the black solid lines in Figure 1. To limit the analysis to these stars only, a cut of $Q\leqslant 0$ is applied. We also have a limit on Q in the high-mass end of Q ≥ −2.1, as this is the lowest Q for which the isochrones have a mass ≤100 M_⊙. Kim et al. (2012) used a shallower reddening slope for these stars to attempt to find a physical Q-value; however, since this should not affect the derived extinctions in the maps from Section 5 significantly, we exclude them for the sake of simplicity. We simultaneously match stars to the 1, 2, 5, 10, 15, 25, 50, 75, and 100 Myr isochrones, choosing the isochrone that contains the Q-value closest to that of each individual star.

Zoom In Zoom Out Reset image size

We can then obtain the intrinsic color directly from the matched isochrone; thus, the color excess for each star becomes

$$\begin{eqnarray}&&E(V-I)={(V-I)}_{\mathrm{observed}}-{(V-I)}_{\mathrm{intrinsic}},\end{eqnarray} \tag{ 5 }$$

and similarly for the (NUV − B) color. We then calculate the extinction values for each wavelength band using a standard MW extinction law (Cardelli et al. 1989) and obtain corrected magnitudes for each star. Stars that are located to the left of the isochrones before the correction are given E(V − I) = 0 if they are within the photometric error of the isochrones, otherwise they are removed from the corrected catalog. Figure 1 shows a sample extinction correction for 200 stars in NGC 628. The extinction correction for each star can be visualized as shifting the star along its constant Q_NBVI line to the locus of the isochrones in the two-color plot; the distance of the shift is proportional to the reddening in each of the two colors on the axes.

Note that the corrected stars in Figure 1 do not lie precisely on the isochrones. The stars in the LEGUS catalog each have some level of photometric uncertainty, which we must take into account in our correction. After matching a star to an isochrone, we apply a shift in color based on the photometric uncertainty in the four bands used for correction. The final colors and extinctions for corrected stars reflect the correction plus the photometric uncertainty of each star.

Unlike for the red clump method, there is no significant shift in the isochrones due to metallicity or age, so our extinction correction should not be age- or metallicity-dependent, as seen in Figure 2. For clarity, we only show the 2 and 100 Myr isochrones, as intermediate ages all lie in the same location. Note that the scale of Figure 2 is zoomed in to make the small shift in the isochrones visible; thus, the error bar appears larger than that in Figure 1.

Zoom In Zoom Out Reset image size

Reines et al. (2010) found a significant contribution from emission lines for both the F555W and F814W HST filters for young massive clusters surrounded by H ii regions. Since we focus on young massive stars, it is possible that our stars are surrounded by or near H ii regions and therefore may suffer some of this contamination. However, unless the H ii region is an unresolved point source directly on top of an individual star, it is unlikely to be classified as a star by DOLPHOT if the contamination is significant. In the worst case, Reines et al. (2010) found that about half of the flux in the F555W filter is from emission-line contributions, while F814W suffers from ∼25% emission-line contribution. Combined, this results in a shift of the observed (V − I) color of ∼0.2, which will increase the observed Q by ∼0.3. These are close to our photometric errors, and the overall effect will be to slightly lower the derived A_V for the affected stars.

It should also be noted that Mason et al. (2001) found a high incidence of binary star systems among young massive stars. Since this extinction correction relies on a star's position in color–color space, we should be largely unaffected. A binary system with one high-mass and one low-mass star will have the color of the high-mass star, since the low-mass star is too faint to significantly affect the SED of the whole system. The overall flux increases, but the color (and therefore the Q-value) remains the same. Thus, we can generate accurate extinction corrections for stars in the LEGUS catalog even if they are blends of multiple stars.

In this study, we only included stars more massive than 2 M_⊙, since our extinction-correction method suffers from the isochrone degeneracy issue for stars less massive than this limit. Theoretically, using the NUV filter allows us to correct all stars with a Q_NBVI ≤ 0; however, practically, photometric uncertainties in the LEGUS sample require us to make an additional apparent magnitude cut of m_NUV ≤ 25 and m_B ≤ 26 before generating an extinction correction.

There is an inherent completeness issue in our analysis: fainter, lower-mass stars cannot be observed at high extinction. Without applying a mass cut in addition to our apparent magnitude cutoff, we will tend to underestimate extinction, since we simply do not see the lower-mass stars that have high extinctions. Table 4 shows the minimum-mass stars that can be corrected for a galaxy given our apparent magnitude limits, assuming no extinction. Table 5 shows the minimum mass of main-sequence stars that can be detected at A_V = 0, 1, 2, and 3 for each galaxy.

Table 5 also implies that for more nearby galaxies, we are seeing a slightly different population of stars, as less-massive stars become visible. Thus, to ensure we are measuring the dust-to-gas ratio across a consistent stellar population, we apply an absolute magnitude cut of −5 in the NUV band. From the PARSEC isochrones, this corresponds to an approximate mass of 14 M_⊙ when combined with our Q_NBVI ≤ 0 cut.

This higher-mass cutoff, however, could create a geometric completeness problem, as we limit the number of usable stars significantly. In Section 6, we require a minimum of 10 stars per pixel in our extinction maps to adequately sample the extinction in that pixel. The reasons for this are twofold: first, a higher number of stars will more accurately sample the line of sight, as a higher number of stars is more likely to sample the entire dust layer. Second, a lower resolution better samples the resolution seen by the neutral gas maps, which smear out gas filaments that might be fully resolved in a higher-resolution map. Stars within the gas filaments likely have a higher extinction than those outside it. Thus, averaging over a similar resolution to the neutral gas maps helps minimize some of this potential variation, as it is more likely to find stars both in and out of gas filaments in the same area with a larger resolution. For determination of a dust-to-gas ratio, we must sample the stellar extinction over areas comparable to the resolution of the radio data. After applying a magnitude cut corresponding to approximately 14 M_⊙, we find that our analysis supports the construction of 10'' × 10'' extinction maps for comparison to the neutral gas maps and calculates the dust-to-gas ratio using 10'' × 10'' extinction maps. Figure 3 shows histograms of the number of pixels with the required number of stars in each galaxy.

Zoom In Zoom Out Reset image size

This introduces a potential bias toward lower extinctions, in that we exclude regions of high extinction where we cannot observe many stars because they are obscured. We test the effects of including pixels with fewer stars in Section 7.1.

It should also be noted that the version of the LEGUS catalogs used in this paper is conservative and includes only stars with the most accurate photometry. This restriction, discussed in Sabbi et al. (2018), excludes ∼10% of stars in each pointing, so for our requirements, we miss perhaps ∼34% of stars. This effect does not change the distribution of stars with respect to apparent magnitude, only the total number of stars, so we do not expect a significant change in our results. The conservative catalog is sufficient to support the development of the extinction maps in Section 5.

We generate extinction maps for each galaxy using the individual stellar extinctions. Stars are binned spatially, then a weighted average based on the uncertainty in an individual star's extinction derivation is calculated for a given pixel. The quality of the map depends on both the resolution and the number of stars per pixel. In principle, we want high resolution; however, due to the uneven distribution of stars in galaxies, this is not always feasible. We also want a high number of stars in each pixel to reduce the effect of stars with high uncertainties on the derived extinctions and potential differences in line-of-sight distances.

To maximize resolution while maintaining enough stars for an accurate measurement of the extinction, we employ an adaptive resolution algorithm. In Figure 4, we show 1'' sampled adaptive resolution extinction maps for each of our fields. Starting in the upper left, we define the smallest possible region that contains 10 stars. The smallest region we consider is 1'' × 1''. If a region does not contain 10 stars, we move on to a 2'' × 2'', 5'' × 5'', and 10'' × 10'' region, stopping when we have at least 10 stars and using the starting 1'' × 1'' area as the center. If there are not 10 stars available even in the 10'' × 10'' case, we do not calculate an extinction for that pixel. We then repeat this for the next pixel position, moving across the row. At the end of the row, we move 1'' down and repeat the process. For example, in the most crowded regions, regions with an adequate number of stars have a pixel size of 1'' × 1''; however, on the outer edge of the galaxy, we may be required to sample stars upward of 5'' from the central point of a pixel in order to meet the minimum required number of stars. The effective resolution of these regions is then 10'' × 10''. This sometimes results in counting the same star in multiple pixels in the sparsest regions, though the effect is to spread the extinction out over a larger area as opposed to skewing the extinction to a higher value. This adaptive map is smoothed using a 3'' Gaussian to remove the sharper discontinuities in some regions.

Zoom In Zoom Out Reset image size

In general, the dwarfs UGC 5139 and UGC 4305 have much lower extinctions than NGC 628 or NGC 6503. We also see good correlation between areas of high extinction in the extinction maps and visible patches of dust in the optical image, such as the center of NGC 628, and the edges of visible dust lanes in the spiral arms of NGC 628 and the edge of NGC 6503. However, we also see large empty regions in the middle of the maps, particularly for NGC 628. These are regions where we have a lack of coverage of stars, not necessarily regions with no extinction. Some, like the very large dust lane in NGC 628, may be regions of extremely high extinction, where stars are dimmed below our ability to generate accurate photometry. Young massive stars tend to be observed offset from these regions, as density-wave progression implies that these stars appear most often downstream from the main dust lanes in spiral arms (Wielen 1974). Elmegreen et al. (2014) did find that there are massive stars in dust lanes; however, they cannot be observed in the optical/NUV, as they are highly obscured. We also have the problem of crowding. For the most crowded regions of the galaxies, point sources are blended together and therefore not resolved. These sources then cannot be extracted by DOLPHOT and are not in our catalogs.

Table 5 implies that for more distant galaxies, we are unable to measure extinctions of more than 2 mag in the V band. However, as seen in Figure 4, we do see extinction map pixels for NGC 628 in the central field with A_V > 2. These regions are near the crowded center of the galaxy, so it is possible that we are seeing a high percentage of blends in this particular case. Blends in more distant galaxies may still appear as point sources and are extracted by DOLPHOT as individual stars when there may be two, three, or even more stars contributing to the light seen from a source. Some may even be very low-mass clusters, leading to a higher-than-expected magnitude for a given color. The extinctions in this crowded region are likely to be unreliable and should therefore be treated with the appropriate caution.

In principle, these maps could be used for correcting Hα and optical images for extinction. They are also useful for qualitative comparisons to emission maps of dust from IR and FIR imaging. In both cases, we recommend caution, as these maps are likely lower limits for the extinction in these galaxies. We show them here for illustrative purposes; however, the analysis in Sections 6 and 7 relies on 10'' × 10'' maps.

The total neutral hydrogen gas column density depends on both the H i and H₂ column density. The H i column density can be taken directly from 21 cm flux maps (assuming optically thin H i), but the H₂ column density must be calculated from CO flux maps using an assumed X_CO. For these maps, we use the same X_CO used by Rahman et al. (2011) and Leroy et al. (2009) of 2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹. The total neutral hydrogen column density in a given line of sight is then

$$\begin{eqnarray}&&{N}_{{\rm{H}}}={N}_{{\rm{H}}{\rm{I}}}+2{N}_{{{\rm{H}}}_{2}}.\end{eqnarray} \tag{ 6 }$$

From Rachford et al. (2009) and Bohlin et al. (1978), for sightlines with a standard MW extinction law, one finds Equation (2). We can test if this relation applies to the galaxies in this paper by spatially correlating the neutral gas maps with the dust maps generated in Section 5. The adaptive resolution maps are useful for qualitative comparisons to the location of spiral arms and dust lanes in the galaxy; however, the radio data are also at a much lower resolution than the adaptive resolution maps. The neutral gas is therefore averaged out over a larger area, which leads to geometric effects if the resolutions are not properly matched. To compare the extinction maps to the neutral gas maps, we simply bin the stars spatially into 10'' × 10'' bins and calculate weighted averages of the extinction for each bin.

We also need to consider that, since we are using exclusively young stars, which usually lie in the midplane of the disk, this only represents half the gas column density along a line of sight. Therefore, we may need an estimate of the maximum extinction derived for individual stars along the line of sight. Since the most highly obscured stars in each bin are also likely the faintest stars, and therefore stars with the highest photometric uncertainties, a standard maximum is not necessarily representative of the maximum value of the pixel. Thus, we calculate the weighted average of the upper quartile of the stars in the bin as the maximum of that bin. This "maximum" is shown in Figure 5 for the central field of NGC 628. We also require that all pixels in the extinction map contain at least 10 stars, to ensure that we are properly sampling the line of sight of a pixel. The upper quartile therefore always consists of at least the three most-obscured stars. Both the weighted average and the "maximum" values for extinction are used to quantify the effect of estimating A_V in this way in the analysis that follows.

Zoom In Zoom Out Reset image size

The recent PHAT survey used the assumption of a log-normal distribution of reddenings along the line of sight to generate the extinction maps (Dalcanton et al. 2015) from individual red giant stars. This more accurately models a thin dust screen. The population of stars that the PHAT survey is sensitive to is much older than that of the LEGUS catalog. Thus, the potential for stars to exist significantly above or below the midplane of the disk (and therefore the majority of the dust) is much higher. Dalcanton et al. (2015) estimated the scale height of M31's stars to be more than a factor of 10 times the likely depth of the dust; thus, a thin screen with a log-normal distribution of dust is more appropriate for their sample. The vast majority of young massive stars are likely distributed with a scale height similar to that of the dust of the disk (and within the dust layer). Since our sample is limited to young massive stars, we may assume a simpler model of a smooth layer of dust within a stellar disk. We expect a symmetric histogram of extinction values around a mean for each 10'' pixel, as long as higher extinctions are not so high that we run out of highly obscured stars. On average, we may expect a typical young star to be in the middle of the disk plane and therefore probe half the extinction implied by the total gas column, in which case the expected relation would be

$$\begin{eqnarray}&&\displaystyle \frac{{A}_{V}}{{N}_{{\rm{H}}}}=2.65\times {10}^{-22}\,\mathrm{mag}\,{\mathrm{cm}}^{2}\,{{\rm{H}}}^{-1}.\end{eqnarray} \tag{ 7 }$$

We show the full-plane and half-plane relations for the MW from Equations (2) and (7) in Figure 6.

Zoom In Zoom Out Reset image size

The H i and CO data are mostly taken from the H i Nearby Galaxy Survey (THINGS), Local Irregulars That Trace Luminosity Extremes (LITTLE) THINGS, and the HERA CO-Line Extragalactic Survey (HERACLES) and supplemented where necessary by other sources. The neutral gas map data sources and beam sizes are summarized in Table 6.

Table 6.  Summary of Neutral Gas Map Data

Galaxy	H i Sourcea	H i Beam Size	CO Sourceb	CO Beam Size
NGC 628	THINGS	1188 × 93	HERACLES	11''
NGC 6503	Greisen et al. (2009)	14''	CARMA-STING	43
NGC 7793	THINGS	156 × 1085	⋯	⋯
UGC 4305 (Ho II)	LITTLE THINGS	1374 × 1257	HERACLES	11''
UGC 5139 (Ho I)	LITTLE THINGS	1466 × 1273	HERACLES	11''

Notes.

^aH i data sources. The H i Nearby Galaxy Sample (THINGS) is described in Walter et al. (2008). Local Irregulars That Trace Luminosity Extremes (LITTLE) THINGS is described in Hunter et al. (2012). ^bCO data sources. The HERA CO-Line Extragalactic Survey (HERACLES) is described in Leroy et al. (2009). The CARMA Survey Toward IR-bright Nearby Galaxies (CARMA-STING) is described in Rahman et al. (2011).

Download table as:  ASCIITypeset image

Figure 6 shows the relation between inferred extinction and gas column for 10'' × 10'' pixels. There is significant scatter in each plot, suggesting influence from geometric effects caused by variations of the location along the line of sight of stars in a pixel. On all of the plots, the vast majority of points lie to the left of the standard Galactic dust-to-gas ratio, as expected. We are unlikely to observe the full column density of dust, as we are unlikely to observe stars on the far side of the disk. In fact, the high number of points to the right of both the half- and full-MW relations in NGC 628's central field may imply a need for a higher dust-to-gas ratio in that region. NGC 628 and NGC 6503 both have most of their scatter around either the half- or full-MW relation. NGC 7793 and the two dwarfs, on the other hand, show almost vertical relations with little to no points at high extinctions for high dust column densities, implying a smaller dust-to-gas ratio. Given our limits shown in Table 4, we certainly could have observed higher extinctions for the higher gas column densities, especially for the dwarfs, but we do not.

The assumption of solar metallicity and an MW extinction law may not hold for all galaxies. The correction in Section 3 assumes both an MW reddening law and a solar metallicity set of isochrones. Individual stars may show a large variation, but these are averaged out, so the overall effect on the final maps is small. Varying the extinction law to that of the Large Magellanic Cloud (LMC) or the Small Magellanic Cloud (SMC) (Gordon et al. 2003) causes an average normalized difference of about 1%–2% in the weighted average extinction across the map, with an rms variation of about 13%. The effect from varying the metallicity to the metallicity of the young population of the LMC or SMC is smaller, about 0.5%, with an rms variation of 7%–9%. Figure 7 shows the effect of varying extinction law and metallicity on the derived A_V for extinction maps in NGC 628. Most pixels remain within about 15% for a change in extinction law and 10% for a change in isochrone metallicity. Thus, the dominant effect on Figure 6 appears to be from changes in the dust-to-gas ratio, not the local extinction law or metallicity.

Zoom In Zoom Out Reset image size

6.1. Radial Variations

We combined the two fields of NGC 628, performed a radial deprojection, and binned pixel areas from the extinction maps into five radial bins. The first four bins are 30'' in size, which corresponds to a physical size of 1.4 kpc. The outer bin encompasses everything not contained within the first four bins and is, for our two pointings, similar in size to the other bins. We took averages of the neutral gas column density and extinction and calculated average dust-to-gas ratios in each bin. As seen in Figure 8, there is a clear decrease of dust-to-gas ratio with increasing radius in the outer bins, with a turnover in the inner bins. A fit to the ratios gives a power law with a slope of −0.031 ± 0.004, slightly shallower than the observed metallicity gradient of −0.044. All of our values lie beneath the typical observed MW relationship, which is expected since we are likely not seeing the full column density of dust using our method.

Zoom In Zoom Out Reset image size

It should be noted that NGC 628 is the farthest galaxy in our sample. At 9.9 Mpc, the HST pixel scale is 1.92 pc, which is likely insufficient to resolve individual stars, especially in crowded regions such as the galactic center. Since we still find extinction values that are consistent with what might be expected from an MW relation despite this, we are likely seeing many blends in this galaxy in particular. Blends can allow us to see high extinctions for pairs or groups of massive stars, as they affect magnitudes, but blends with Q_NBVI values in the range of massive stars do still seem to enable us to derive consistent results in comparison to the other galaxies. Blends that produce redder Qs (like those produced by red supergiant stars) are excluded by our selection.

We performed a similar radial analysis with NGC 7793, with bins of the same physical size as those in NGC 628. NGC 7793 is a much smaller galaxy than NGC 628 in terms of physical size, so we were restricted to two bins instead of five. As we show in Figure 9, there is very little change in the dust-to-gas ratios between the two bins, as we might expect since we cover less than half the radial distance of NGC 628. Both bins also have significantly lower dust-to-gas ratios than NGC 628, which we might expect from NGC 7793's lower metallicity. As seen in Figure 6 and Table 5, it is possible to observe higher extinctions in NGC 7793, if they exist.

Zoom In Zoom Out Reset image size

6.2. Variations with Metallicity

Using the metallicity gradient for NGC 628 from Berg et al. (2012), we assign a metallicity for each radial bin defined in Section 6.1. We do the same for the two bins in NGC 7793 using the metallicity gradient from Stanghellini et al. (2015). For NGC 6503, UGC 4305, and UGC 5139, we calculate the dust-to-gas ratio for the whole galaxy, since either there are no published metallicity gradients (NGC 6503 and UGC 5139) or there is not enough coverage in the extinction map to effectively bin the galaxy into radial bins (UGC 4305). There are published average metallicities from Croxall et al. (2009) for UGC 4305 and UGC 5139 and from Tikhonov (2014) for NGC 6503. We plot the average dust-to-gas ratios calculated with ${N}_{{\rm{H}}}={N}_{{\rm{H}}{\rm{I}}}+2{N}_{{{\rm{H}}}_{2}}$ against metallicity in the top plot in Figure 10. The scatter in Figure 6 indicates that the weighted average extinction likely corresponds to a slightly different location in the line of sight in each individual pixel. Some pixels may inherently have more stars behind or in front of the dust layer. An average of the A_V and N_H will still be related to the average dust-to-gas ratio modulo a factor of two for those galaxies where we can probe through the dust layer (see Table 5 for estimates of maximum observable A_V). The factor might be larger than two for more dusty or distant galaxies where we do not fully sample the dust layer. We calculate the average A_V of all the pixels with at least 10 stars and an average N_H of their correlated gas column densities, then we calculate A_V/N_H from these average values.

Zoom In Zoom Out Reset image size

There is a clear trend of decreasing dust-to-gas ratio with decreasing metallicity. The lines in the plot are linear fits to the data on a logarithmic scale, implying a power-law relationship between metallicity and dust-to-gas ratio. The dashed line is the fit to the spiral galaxies only, while the solid line shows the fit to the entire sample of points. We see very little scatter in our sample given the fairly high uncertainty in the literature metallicity values. This could imply an overestimation in these metallicity uncertainties.

We do not have access to CO data of a similar resolution to the H i data for NGC 7793, so its N_H values are lower limits. The dust-to-gas ratios for this galaxy are therefore upper limits. UGC 5139 and UGC 4305 also both have negligible H₂ column densities compared to the H i. This may be more an effect of our choice of a constant X_CO factor, so it is possible that the dust-to-gas ratios for these two galaxies are also upper limits. We therefore plot the dust-to-gas ratio for each galaxy calculated using the H i column density against the metallicity in the bottom plot in Figure 10. We still observe a correlation between dust-to-gas ratio and metallicity; however, as seen in Table 7, the slope is steeper than in the full N_H case. Clearly, ${N}_{{{\rm{H}}}_{2}}$ is important for determining the relationship between dust-to-gas ratio and metallicity. If we underestimate ${N}_{{{\rm{H}}}_{2}}$ in low-metallicity galaxies by using a constant X_CO, the dust-to-gas ratios for these galaxies may in fact be upper limits, and the slope of the relation may be steeper than what we observe. This is discussed further in Section 7.1.

Table 7.  Dependence of the Dust-to-gas Ratio on Metallicity

	NH		NH i Only
	Average	"Maximum"	Average	"Maximum"
Spirals only	1.22 ± 0.16	1.30 ± 0.12	2.39 ± 0.17	2.47 ± 0.18
Full sample	1.46 ± 0.06	1.30 ± 0.04	1.82 ± 0.07	1.66 ± 0.08

Note. All values are the slopes of the linear fits from Figure 10.

Download table as:  ASCIITypeset image

Table 7 gives the slopes of the linear fits from Figure 10. Rémy-Ruyer et al. (2014) gave a slope of 1.6 ± 0.3 for a single power law and used a slope of 1 as determined by previous studies (James et al. 2002; Draine et al. 2007; Galliano et al. 2008; Leroy et al. 2011) for the high-metallicity end ($12+\mathrm{log}({\rm{O}}/{\rm{H}})\gt 7.96\pm 0.47$). All of our points lie within this metallicity regime, so it is difficult to determine if a broken power law is necessary to explain the changing dust-to-gas ratio with metallicity.

7.1. Dust-to-gas Ratio versus Metallicity

7.2. Comparison to Extinctions Derived from Spectral Energy Fitting

BEAST is an SED-fitting method used by the PHAT survey to calculate individual stellar extinctions in M31 using multiband photometry. Gordon et al. (2016) recently made this software publicly available. It serves as a useful comparison to our individual extinctions calculated from the isochrone-matching method described in Section 3 for massive stars, since both methods use the same input data and work with individual massive stars. Other methods of extinction correction (i.e., from IR emission maps) might not specifically probe the massive stars or might probe different regions of gas and dust. A comparison with BEAST allows a star-by-star comparison between methods. BEAST allows the simultaneous fitting of several different extinction parameters, including A_V and R_V. Here we compare the BEAST-derived extinctions to those of the isochrone-matched extinctions from Section 3 for the eastern field of NGC 7793. We employ identical magnitude- and Q-based cuts to both the isochrone-matched stars and the BEAST-derived stars. This is the closest spiral in our sample, so it should be least affected by crowding, and the extinction law should be fairly MW-similar (i.e., not starburst-type, as in Calzetti et al. 1994). For our comparison, we restricted the BEAST extinction law to an MW-type extinction law with R_V = 3.1, which corresponds to the reddening slope used in Section 3. The SED fitting was only performed on stars with robust PSF fits in all five bands with photometric errors ≤0.25 mag.

Figure 11 shows a comparison between the calculated A_V for individual stars for our isochrone-matching method and the BEAST SED-fitting method, color-coded by Q_NBVI, for both the entire corrected sample of BEAST-fitted stars and for BEAST-fitted stars above the absolute magnitude cut used in Section 5 to generate the extinction maps.

Zoom In Zoom Out Reset image size

BEAST tends to assign systematically lower A_V to stars compared to the isochrone method, especially at higher extinctions. There is also a large population of stars that BEAST assigns a low mass, but the photometric data show a Q_NBVI and two-color values that place the stars in the range of massive stars, as seen in Figures 11 and 12. Since there are generally many more low-mass stars than high-mass stars, it is possible that BEAST automatically assumes, given photometric equivalence, that the star is low-mass as opposed to high-mass. After applying the absolute magnitude cut used in Section 5, many of the stars in this grouping are removed. Figure 13 also shows that the stars with the greatest differences tend to be the stars with the highest photometric errors, especially in the population of stars with low assigned BEAST masses and Q_NBVI values corresponding to high masses. The systematic difference between derived extinctions using BEAST and isochrone matching is even more obvious in Figure 14, which shows the distribution of A_V(isochrone)–A_V(BEAST) for individual stars. The peak of the distribution is positive, with a large tail further to the positive.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An additional source of uncertainty is the effect of post-main-sequence stars, which can vary in color by as much as 0.1. This is within our photometric error, so it should not be much of an effect, but if BEAST puts more weight on the post-main-sequence, then it could account for at least some of the systematic difference.

The average extinction map generated using the BEAST extinctions reflects that the BEAST individual stellar extinctions tend to be systematically lower than the isochrone-matched extinctions. This can lead to differences in the extinctions of about 50% across the map. We have thus far only tested a single field of one galaxy and do not yet know if this will affect the relative changes in the dust-to-gas ratio across galaxies. More study is needed to examine the exact nature of the differences in extinction and if this behavior is consistent across other galaxies, specifically dwarf galaxies where we already find very low extinctions.

We correct the photometry of massive individual stars for extinction using the isochrone-matching method used by Kim et al. (2012) and compare these extinction corrections to those derived from SED fitting for a single pointing. We generate extinction maps using the extinction values of individual massive stars, averaged over spatial bins to account for the range of optical depths along the line of sight. This is the first time that such maps have been derived using massive stars for galaxies outside the Local Group, and this method can be expanded to other galaxies within 10 Mpc. Extinction maps such as these may be useful in probing cold and dark gas and dust that may be missed in emission (Paradis et al. 2012). We correlate the extinction maps with neutral gas maps in order to investigate the dust-to-gas mass ratio in NGC 628, NGC 7793, NGC 6503, UGC 4305, and UGC 5139. Our results can be summarized as follows.

• 1.  
  We find lower extinctions across the dwarf galaxies UGC 4305 and UGC 5139 than in the spiral galaxies NGC 628, NGC 7793, and NGC 6503.
• 2.  
  There is wide scatter in the localized dust-to-gas ratios for each galaxy. This could be explained by geometric effects caused by variations in the position of stars along the line of sight.
• 3.  
  NGC 628 shows a significant radial variation in dust-to-gas ratio with galactocentric radius, likely correlated with its metallicity gradient.
• 4.  
  We find strong evidence to support a power-law relationship between dust-to-gas ratio and metallicity consistent with other work. A metallicity-dependent CO-to-H₂ conversion factor X_CO would increase the slope compared to what is observed here.
• 5.  
  There is a systematic difference in the extinctions derived via isochrone matching and those derived using the BEAST SED-fitting algorithm. The cause of this systematic difference may be related to a preference within BEAST to assign a low mass to stars for which the Q-value would suggest a high mass.

Further investigations will provide a better understanding of the relationship between dust-to-gas ratio and metallicity, as a larger sample of LEGUS galaxies will expand and fill in our metallicity range. We will expand our comparison of the SED-fitting derived extinctions with isochrone-matching derived extinctions. We will investigate the possibility of a relationship between dust-to-gas ratio and galaxy morphology and star formation rate. If either of these factors influence the dust-to-gas ratio, it could explain the scatter seen in the dust-to-gas mass ratio seen across our extinction maps, since the uncertainties in our individual pixels can only account for ${\sigma }_{{A}_{V}}\sim 0.1$.

Our 1'' × 1'' gridded adaptive resolution maps have a variety of applications. They can be used to correct Hα and other optical images for extinction, including correcting stars in the LEGUS catalogs that cannot be corrected with isochrone matching. These stars may have been below our magnitude or Q-value cutoffs, but if they are located near stars corrected with isochrone matching, they likely have similar extinctions, allowing us to apply the extinctions from the map. Corrections to Hα maps are useful for determining accurate star formation rates. It is also possible to match these maps to maps of the diffuse ionized gas (DIG) and H ii regions in order to study differences in extinction between these two important phases of the ionized ISM. This is important for determining the fractional contribution of the DIG to the total Hα luminosity, which in turn can help determine ionization sources for the DIG (Walterbos & Braun 1994).

This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy under NASA contract NAS 5-26555. These observations are associated with program 13364 (LEGUS). Support for program 13364 was provided by NASA through a grant from the Space Telescope Science Institute. LK and RAW acknowledge support from the Space Telescope Science institute (grant Numbers HST-GO-13773.003-A). LK would additionally like to acknowledge support from the Gemini Observatory in La Serena, Chile. MF acknowledges support from the Science and Technology Facilities Council (grant Number ST/P000541/1). DAG kindly acknowledges financial support by the German Research Foundation (DFG) through program GO 1659/3-2. GA acknowledges support from the Science and Technology Facilities Council (ST/P000541/1 and ST/M503472/1). We would also like to thank the anonymous referee for useful comments that helped improve the paper.

Software: BEAST (Gordon et al. 2016), DOLPHOT photometric package (Dolphin & Kennicutt 2002), PARSEC (Bressan et al. 2012).